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

Graphene has attracted tremendous interest due to its ability to support tightly confined, electrostatically tunable SPs [17–24]. More recently, its nonlinear properties have gained attention [25–29]. For example, four-wave mixing produced by single-pass transmission through a single graphene layer has been observed [26], while a second-order response at oblique incidence angles has been predicted [27], and intrinsic second-order nonlinearities have been used to excite graphene plasmons from free-space beams via difference frequency generation [28]. It has also been proposed that graphene nanostructures could enable quantum third-order nonlinearities [29].

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 [25, 29–31], which describes the evolution of the carrier distribution function ${f}_{{\bf{k}}}({\bf{r}},t)$ at position ${\bf{r}}$ and momentum ${\bf{k}}$. Within this theory ${\bf{k}}$ and ${\bf{r}}$ obey the classical equations of motion: $\dot{{\bf{r}}}={{\bf{v}}}_{{\bf{k}}}=(1/{\hslash })\partial {\epsilon }_{{\bf{k}}}/\partial {\bf{k}}$, and ${\hslash }\dot{{\bf{k}}}=-e{\bf{E}}$. We are interested in energy scales $\lesssim 1$ eV, so it is possible to linearize the graphene dispersion relation around the Dirac points, ${\epsilon }_{{\bf{k}}}=\pm {\hslash }{v}_{F}| {\bf{k}}|$, where +(−) denotes doping to positive (negative) Fermi energies E_F. The single-band approximation holds provided that the optical frequency is less than ∼$2{E}_{F}$, such that absorption arising from interband electron–hole transitions is suppressed [20]. The carrier dynamics are then described by the equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}{f}_{{\bf{k}}}({\bf{r}},t)\pm {v}_{F}\hat{{\bf{k}}}\cdot {\nabla }_{{\bf{r}}}{f}_{{\bf{k}}}({\bf{r}},t)=\displaystyle \frac{e}{{\hslash }}{\bf{E}}({\bf{r}},t)\cdot {\nabla }_{{\bf{k}}}{f}_{{\bf{k}}}({\bf{r}},t),\end{eqnarray} \tag{ 1 }$$

where ${\bf{E}}$ is the sum of the external field ${{\bf{E}}}^{\mathrm{ext}}$ and the induced field ${{\bf{E}}}^{\mathrm{ind}}$ generated by the carrier distribution.

The macroscopic quantities such as the density of charge and the surface current can be related to the microscopic dynamics of the carriers. For instance, the surface current depends on the microscopic carrier velocities as

$$\begin{eqnarray}&&{\bf{J}}({\bf{r}},t)=-{{eg}}_{v}{g}_{s}\int \displaystyle \frac{{{\rm{d}}}^{2}{\bf{k}}}{{(2\pi )}^{2}}{{\bf{v}}}_{{\bf{k}}}{f}_{{\bf{k}}}({\bf{r}},t),\end{eqnarray} \tag{ 2 }$$

where g_s = g_v = 2 are the spin and valley degeneracies of graphene. The nonlinear set of equations (1) and (2) can be solved perturbatively to give the relation between the surface current and the electric field (i.e. the conductivity). At lowest order, one assumes that ${f}_{{\bf{k}}}$ is slightly displaced from its equilibrium (zero temperature) Fermi distribution, ${f}_{{\bf{k}}}^{(0)}({\bf{r}},t)=\theta ({k}_{F}-k)$. Thus, one can substitute ${f}_{{\bf{k}}}^{(0)}$ into the term ${\nabla }_{{\bf{k}}}{f}_{{\bf{k}}}$ of equation (1), yielding a linear relationship between a perturbed distribution function ${f}^{(1)}$ and ${\bf{E}}$. Solving in the Fourier domain, the perturbed distribution function ${f}^{(1)}$ can be inserted into equation (2) to find the resulting current. This yields a linear relation between the current and field, ${\bf{J}}(k,\omega )={\sigma }^{(1)}(\omega ){\bf{E}}(k,\omega )$, where the proportionality constant is the familiar Drude conductivity [18, 19]

$$\begin{eqnarray}&&{\sigma }^{(1)}(\omega )=\displaystyle \frac{{\rm{i}}{e}^{2}| {E}_{F}| }{\pi {{\hslash }}^{2}\omega }.\end{eqnarray} \tag{ 3 }$$

Finally, the procedure can be iterated by inserting the perturbed distribution function (e.g., ${f}^{(1)}$) into equation (1) to yield higher order conductivity functions, as derived in greater detail in appendix A.

Graphene is a centro-symmetric material, which is typically associated with a vanishing second-order nonlinearity [32]. Indeed, if the nonlinear response is spatially local, ${J}^{(2)}(2\omega ,{\bf{r}})={\sigma }^{(2)}(\omega )E{(\omega ,{\bf{r}})}^{2}$, spatial inversion symmetry implies that $-J={\sigma }^{(2)}{(-E)}^{2}$, which enforces that ${\sigma }^{(2)}=0$. This argument breaks down if the conductivity is nonlocal [33], for example if $\sigma (\omega ,q)\propto q$, such that the current depends on the electric field gradient, ${J}^{(2)}={\sigma }^{(2)}(\omega )E{\partial }_{{\bf{r}}}E$.

In principle, nonlocal effects are present in any material. For a given electric field strength, the size of this nonlinear effect depends on a dimensionless parameter $k{/k}_{{nl}}$ [34]. Here k is the wavevector of the light that dictates how rapidly the field changes in space, and ${k}_{{nl}}^{-1}$ is a characteristic length scale over which carriers in the material become sensitive to field gradients. In materials where the charges are tightly bound to their atoms, the relevant length scale ${k}_{{nl}}^{-1}$ is given by the atomic size of angstroms, which is thus negligible compared to optical wavelengths. In conducting materials, the length scale is set by the typical distance between carriers, which is proportional to the inverse of the Fermi wavevector. In a typical metal like silver, the high carrier density also yields a negligible length scale of ${k}_{{nl}}^{-1}\sim {k}_{F}^{-1}\sim 1$ Å. In contrast, in graphene we can simultaneously exploit two effects to increase significantly $k{/k}_{{nl}}$. First, graphene can be electrostatically tuned to have very low carrier densities to increase ${k}_{F}^{-1}$. Second, one can use tightly confined plasmon excitations in graphene, which have been shown to yield a reduction in the wavelength (or equivalently enhancement in wavevector q_p) compared to free-space light by two orders of magnitude. Indeed, below we show specifically that $k{/k}_{{nl}}\sim {q}_{p}{/k}_{F}\lesssim 1$ emerges as the relevant quantity to characterize the strength of nonlocal nonlinearities in graphene.

After these considerations, we calculate the second-order conductivity using the procedure explained above. The second-order conductivity can be expanded in powers of q_p in the long-wavelength limit, defined by the condition ${v}_{F}{q}_{p}/{\omega }_{p}\ll 1$ (see appendix A). 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 ${\omega }_{p}$ and an induced current density at frequency $2{\omega }_{p}$

$$\begin{eqnarray}&&{J}_{i}^{2{\omega }_{p}}={\sigma }^{(1)}(2{\omega }_{p}){E}_{i}^{2{\omega }_{p}}+{\sigma }_{{ijkl}}^{(2)}(2{\omega }_{p};{\omega }_{p}){E}_{j}^{{\omega }_{p}}{\nabla }_{k}{E}_{l}^{{\omega }_{p}}.\end{eqnarray} \tag{ 4 }$$

Here ijkl denote in-plane vector indices and summation over repeated indices is implied. The nonlocal second order conductivity tensor reads

$$\begin{eqnarray}&&{\sigma }_{{ijkl}}^{(2)}(2{\omega }_{p};{\omega }_{p})=\mp \displaystyle \frac{{\rm{i}}{e}^{3}{g}_{v}{g}_{s}{v}_{F}^{2}}{32\pi {{\hslash }}^{2}{\omega }_{p}^{3}}(5{\delta }_{{ij}}{\delta }_{{kl}}-3{\delta }_{{ik}}{\delta }_{{jl}}+{\delta }_{{il}}{\delta }_{{jk}}).\end{eqnarray} \tag{ 5 }$$

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 [27].

The Drude conductivity for infinite graphene given by equation (3) provides a valid description of the carrier dynamics of graphene when ${\hslash }\omega \lesssim {E}_{F}$ [18, 19], where the interband transitions can be neglected. Like any conductor in contact with a dielectric (or vacuum, as we assume here), graphene supports SPs with a dispersion relation given by

$$\begin{eqnarray}&&\displaystyle \frac{{q}_{0}}{{q}_{p}}\approx 2\alpha \;\displaystyle \frac{{E}_{F}}{{\hslash }{\omega }_{p}},\end{eqnarray} \tag{ 6 }$$

where ${q}_{0}=\omega /c$ is the free-space wavevector at the same frequency and $\alpha \approx 1/137$ is the fine structure constant. As ${E}_{F}\gtrsim {\hslash }{\omega }_{p}$, equation (6) indicates a reduction in the plasmon wavelength compared to free space by up to two orders of magnitude, which should significantly drive up the effects of spatially nonlocal interactions.

We have seen that at fixed field strength, the nonlinear interactions between plasmons in graphene should be increased due to a large ratio of ${q}_{p}{/k}_{F}$. However, what is most important for nonlinear optics is how to maximize the interaction strength per photon (i.e. per quantized plasmon). A simple argument, made more precise below, is that because the energy of a single plasmon is fixed at ${\hslash }{\omega }_{p}$, confining it to as small volumes V as possible maximizes its intensity or electric field, ${E}_{{\rm{0}}}\sim \sqrt{{\hslash }{\omega }_{p}/{\epsilon }_{0}V}$. This motivates the study of nonlinear optical interactions between plasmons in nano-structures, which we now present in detail. As a specific example, we will focus on nanostructures that have plasmon resonances at frequencies ${\omega }_{p}$ and $2{\omega }_{p}$. This particular choice of structure is to facilitate DC or SHG.

The derivation of the quantum Hamiltonian of the system (reported in appendix B) starts from the expression of the electrostatic energy (a valid approach provided that the linear dimension of the structure D is small compared to the free-space wavelength ${\lambda }_{0}$)

$$\begin{eqnarray}&&H=\displaystyle \frac{1}{2}\displaystyle \sum _{{\omega }_{i}}{\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{{\rho }^{{\omega }_{i}}}^{*}({\bf{r}}){\phi }^{{\omega }_{i}}({\bf{r}}),\end{eqnarray} \tag{ 7 }$$

where ρ is the charge density and ϕ the electrostatic potential. The charge density can be replaced by the current density using the continuity equation, which in turn can be expressed in term of the electric field using equation (4). Expressing the potentials in terms of the electric fields we end up with an expression for the energy depending only on the electric field of the two modes, whose canonical quantization gives

$$\begin{eqnarray}&&H={\hslash }({\omega }_{p}-{\rm{i}}{\Gamma }_{a}^{\prime }/2){a}^{\dagger }a+{\hslash }(2{\omega }_{p}-{\rm{i}}{\Gamma }_{b}^{\prime }/2){b}^{\dagger }b+{\hslash }g(\;{b}^{\dagger }{a}^{2}+{\rm{h}}.{\rm{c}}.).\end{eqnarray} \tag{ 8 }$$

Here 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{\omega }_{p}$ and two plasmons with frequency ${\omega }_{p}$ [35]. Adopting a quantum jump approach we have added to the frequencies an imaginary part accounting for the total decay rates ${\Gamma }_{a}^{\prime }$ and ${\Gamma }_{b}^{\prime }$ of the two modes. The quantization associates with a single plasmon a typical electric field amplitude ${\overset{}{{\rm{E}}}}_{0}^{{\omega }_{p}}\sim {({\hslash }{\omega }_{p}{q}_{p}/{\epsilon }_{0}S)}^{1/2}$, where S is the structure area, confirming the large per-plasmon field associated with tight confinement. The quantum coupling constant g is rigorously given by the classical interaction energy between the nonlinear current at $2{\omega }_{p}$ and the fields at ${\omega }_{p}$, but with the classical field values replaced by the per-photon field strengths ${\tilde{{\bf{E}}}}^{{\omega }_{i}}({\bf{r}})$

$$\begin{eqnarray}&&{\hslash }g=| \displaystyle \frac{1}{4{\rm{i}}{\omega }_{p}}\;{\sigma }_{{ijkl}}^{(2)}(2{\omega }_{p};{\omega }_{p}){\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{\tilde{E}}_{i}^{2{\omega }_{p}}({\bf{r}}){\tilde{E}}_{j}^{{\omega }_{p}}({\bf{r}}){\nabla }_{k}{\tilde{E}}_{l}^{{\omega }_{p}}({\bf{r}})| .\end{eqnarray} \tag{ 9 }$$

Equation (9) shows that g is directly proportional to the second-order conductivity ${\sigma }_{{ijkl}}^{(2)}$ calculated in the previous section, and its dependence on the particular geometric configuration of the modes is confined to the overlap integral [36]. It should be noted that for extended graphene, the mode functions are simply propagating plane waves $E({\bf{r}})\sim {e}^{{\rm{i}}{kz}}$. Thus the integral in equation (9) produces a delta function, $g\propto \delta (2{k}_{1}-{k}_{2})$, which reflects momentum conservation. In contrast, in small structures the spatially complex modes can be thought of as a superposition of many different wavevectors, and a large interaction strength is ensured by engineering the modes such that they have good spatial overlap [37].

Using the fact that ${\overset{}{{\rm{E}}}}_{0}^{{\omega }_{p}}\sim {({\hslash }{\omega }_{p}{q}_{p}/{\epsilon }_{0}S)}^{1/2}$, that the nonlinear conductivity has an amplitude ${\sigma }^{(2)}\sim {e}^{3}{v}_{F}^{2}/{{\hslash }}^{2}{\omega }_{p}^{3}$, and that the field gradients occur over a length scale ${q}_{p}^{-1}$, one can readily verify that equation (9) predicts a general scaling of $g/{\omega }_{p}=\beta /{({k}_{F}D)}^{7/4}$. The dimensionless coefficient of proportionality, which we call β, depends only on the geometric overlap of the modes (e.g. $\beta =0$ if the modes have the wrong symmetries, or $\beta \sim 1$ for modes with good overlap). As the minimum dimension of the structure should be comparable to the plasmon wavelength, $D\sim 1{/q}_{p}$, the maximum ratio of $g{/\omega }_{p}$ scales like ${({q}_{p}{/k}_{F})}^{7/4}$, confirming the enhanced nonlinearities as q_p become comparable to k_F. Note that this relation is valid only for ${q}_{p}\lesssim {k}_{F}$, where the conductivity of graphene is Drude-like, as discussed above. In this derivation, we have assumed that a finite-size structure has the same conductivity as infinite graphene. Although this is not true for arbitrarily small structures, where quantum finite-size effects play a significant role, this approximation is already qualitatively correct for structures with $D\gtrsim 10$ nm [38].

To show that a high overlap factor of $\beta \sim 1$ can be reached in typical structures, we consider one specific example of a doped graphene isosceles triangle embedded in vacuum. This choice enables a simple optimization to obtain the desired ratio of 2 between the SP mode frequencies. Indeed, we find that an aspect ratio r = 1.3 produces plasmons at frequencies ${\omega }_{p}$ and $2{\omega }_{p}$ (see figure 1). The modes shown in figure 1 are numerically computed using a commercial finite-difference code (${\mathrm{COMSOL}}^{\circledR }$) by driving the system with a plane wave whose associated external field ${{\bf{E}}}^{\mathrm{ext}}$ is polarized along the axis of symmetry of the triangle. We model the structure as a thin slab with rounded edges and a dielectric function $\epsilon =1+4\pi {\rm{i}}{\sigma }^{(1)}/\omega t$. The thickness t is chosen to be t=0.5 nm (this value is sufficiently small that the in-plane current has converged, and the results do not depend on the specific value), and the expression of ${\sigma }^{(1)}$ is given by the equation (3). Since the characteristic length of the structure is much smaller than the free-space wavelength, the response can be determined electrostatically, where the retardation and the response to the magnetic field are neglected. Furthermore, the ratio 1 : 2 between the first and second plasmon resonances is preserved independently of the actual size of the triangle and the doping [39]. While the remaining parameters are somewhat arbitrary, as a numerical example, we consider the realistically achievable length and doping level of D = 22nm and E_F = 0.2 eV. For this choice, we observe a pronounced first harmonic mode (figures 1(a), (c)) with energy ${\hslash }{\omega }_{p}\simeq 0.20$ eV, and a second harmonic resonance (figures 1(b), (c)) twice as energetic. Once we obtain the mode profiles, their nonlinear coupling is evaluated using the equations (5) and (9). Numerical calculations for this structure yield a value of $\beta =0.34$, hence the quantum oscillation rate g reaches a remarkable 1.25% value of the dipolar frequency ${\omega }_{p}$.

Zoom In Zoom Out Reset image size

Surface plasmons in realistic graphene structures generally decay by non-radiative mechanisms, whose precise nature is still under active investigation [23, 40, 41]. We thus use a phenomenological description associating an intrinsic decay rate Γ to the modes. For our numerical calculations we will assume a mode quality factor of $Q={\omega }_{p}/\Gamma$ ranging from some tens to one hundred, close to what has been experimentally observed in nanostructures [23], although in our analytical results we will explicit keep track of the scaling with Γ.

In addition to intrinsic decay channels, graphene SPs can also be excited and detected through desirable channels, i.e. via radiative decay. We will use the notation ${\kappa }_{a,b}$ to indicate such decay rates. The total decay rate introduced in Hamiltonian (8) is thus ${\Gamma }_{a,b}^{\prime }={\Gamma }_{a,b}+{\kappa }_{a,b}$. We will also introduce the notation ${\eta }_{a,b}$ to indicate the external coupling efficiencies of the modes, defined as ${\kappa }_{a,b}/{\Gamma }_{a,b}^{\prime }$. For example, in our structure, the first and second harmonic modes radiate into free space at rates ${\kappa }_{a}\approx 2\times {10}^{-7}{\omega }_{p}$ and ${\kappa }_{b}\approx 5.4\times {10}^{-8}{\omega }_{p}$, as numerically calculated through the extinction cross sections of the incident field (see appendix C). The external coupling efficiency can be increased by using more sophisticated techniques, such as SNOM [21] or graphene nanoribbons [29].

The rate of oscillation or internal conversion between a single quantized plasmon and two lower-frequency plasmons is remarkable, particularly considering that the state-of-the-art down-conversion efficiency in conventional nonlinear crystals is $\sim {10}^{-8}$ [14, 15]. It should be pointed out that the internal conversion rate holds independently of how the plasmons are generated. Of course, for both practical observation and for technological relevance, it would be ideal if the plasmons could be efficiently excited and subsequently converted back into propagating photons (such as from free space, fiber, or other evanescent modes). Motivated by this, we now examine the coupling problem to propagating photons in more detail and investigate how their intermediate conversion and interaction as plasmons manifests itself as strong, effective nonlinearities between propagating photons.

Remarkably, the extinction cross section ${\sigma }^{\mathrm{ext}}=(3/2\pi ){\lambda }_{0}^{2}\kappa /\Gamma ^{\prime}$ of a single nano-structure (see appendix C) can exceed its physical size. However, the low values of $\kappa /\Gamma ^{\prime}$ still imply that ${\sigma }^{\mathrm{ext}}$ is much smaller than the diffraction limited area ${\lambda }_{0}^{2}$ for free-space beams, indicating that such sources cannot be used to excite plasmons efficiently. In particular, it can be shown using time-reversal symmetry that the best in-coupling (excitation) efficiency that can be achieved is the same as the out-coupling efficiency, η [42]. The situation is illustrated schematically in figure 2. This raises an important conceptual question. On one hand, graphene plasmons seem to represent the 'ultimate' quantum nonlinear optical device, capable of internal conversion at the single-photon level. However, very little incoming light enters the structure and turns into a plasmon, and vice versa, a small percentage of plasmons are radiated back into light. We now discuss various ways in which the strong quantum-level internal nonlinearities of graphene can be observed and utilized, given these limitations.

Zoom In Zoom Out Reset image size

One way of increasing the coupling to radiation, which has already been discussed in the linear optical regime, is to exploit an array of nano-structures [23, 43]. Intuitively, since the extinction cross section of a single element can exceed its physical size, having a dense array extending over an area larger than ${\lambda }_{0}^{2}$ guarantees efficient interaction with an incoming beam. We thus proceed to consider the nonlinear interaction between an incoming radiation field with frequency ${\omega }_{p}$ resonant with the fundamental mode and an array of nano-structures, as illustrated in figure 3. We expect that the efficient coupling with an array will enable the incoming photons to excite plasmons at ${\omega }_{p}$, internally convert to plasmons at $2{\omega }_{p}$, and then re-radiate into free-space as a second harmonic signal. We consider here a hexagonal lattice of nanostructures with lattice period l = 50 nm. The array is illuminated at normal incidence with a field of frequency ω, and polarized along $\hat{x}$ to maximally drive the plasmon resonance (see figure 3).

Zoom In Zoom Out Reset image size

From Hamiltonian (8) extended to include the coupling between the structures, we get the equations of motion of the operators for the first and second harmonic modes of structure i in the array are

$$\begin{eqnarray}{\dot{a}}_{i} & = & -{\rm{i}}({\omega }_{p}-{\rm{i}}{\Gamma }_{a}^{\prime }/2){a}_{i}-{\rm{i}}\displaystyle \frac{{p}_{a}}{{\hslash }}{E}_{\omega }^{\mathrm{ext}}-2{\rm{i}}g{a}_{{\rm{i}}}^{\dagger }{b}_{{\rm{i}}}+{\rm{i}}\;\displaystyle \frac{{p}_{a}^{2}}{{\hslash }}\displaystyle \sum _{j}\;{G}_{{ij}}^{{\omega }_{p}}{a}_{j},\\ {\dot{b}}_{i} & = & -{\rm{i}}(2{\omega }_{p}-{\rm{i}}{\Gamma }_{b}^{\prime }/2){b}_{i}-{\rm{i}}\displaystyle \frac{{p}_{b}}{{\hslash }}{E}_{2\omega }^{\mathrm{ext}}-{\rm{i}}g{a}_{{\rm{i}}}^{2}+{\rm{i}}\;\displaystyle \frac{{p}_{b}^{2}}{{\hslash }}\displaystyle \sum _{j}\;{G}_{{ij}}^{2{\omega }_{p}}{b}_{j},\end{eqnarray} \tag{ 10 }$$

where the last term in both equations accounts for the dipole–dipole interaction with other nanostructures j in the array, and ${G}_{{ij}}=G({r}_{i},{r}_{j})$ is the electromagnetic Green's function describing the field produced at position r_i by a dimensionless dipole oscillating at r_j, while ${p}_{a}=\sqrt{3\pi {\epsilon }_{0}{\hslash }{\kappa }_{a}{c}^{3}{/\omega }_{p}^{3}}$ is the modulus of the electric dipole moment of a single plasmon in the first mode (an equivalent expression holds for p_b at frequency $2{\omega }_{p}$). We have also included the possibility of driving either mode with classical free-space external fields, denoted by ${E}_{\omega }^{\mathrm{ext}}$ and ${E}_{2\omega }^{\mathrm{ext}}$.

Before considering the generation of a second harmonic, it is already interesting to point out that the strong internal interactions between plasmons can manifest itself in the linear optical response to an incoming laser with frequency near the second mode $2{\omega }_{p}$. We proceed by solving the coupled system of equations (10) for a weak external driving field of frequency ω around $2{\omega }_{p}$. We consider specifically an approximation where edge effects are ignored (which becomes exact in the plane-wave limit and an infinite array), which makes the sum ${\displaystyle \sum }_{j}{G}_{{ij}}$ identical for each element. As discussed in detail in appendix D the effect of the Green function is to renormalize both the resonance frequencies and the losses, so that ${\omega }_{p}\to {\tilde{\omega }}_{p},{\Gamma }_{a}^{\prime }\to {\tilde{\Gamma }}_{a}^{\prime }$, etc. We find that the linear reflection coefficient of the array is

$$\begin{eqnarray}&&{r}_{b}(\omega )=-\displaystyle \frac{{\rm{i}}{\kappa }_{b}{N}_{{\lambda }_{0}^{2}}}{2}\;\displaystyle \frac{{\tilde{\delta }}_{a}+{\rm{i}}{\tilde{\Gamma }}_{a}^{\prime }}{[{\tilde{\delta }}_{a}+{\rm{i}}{\tilde{\Gamma }}_{a}^{\prime }][{\tilde{\delta }}_{b}+{\rm{i}}{\tilde{\Gamma }}_{b}^{\prime }/2]-2{g}^{2}},\end{eqnarray} \tag{ 11 }$$

where ${\tilde{\delta }}_{a}=\omega -2\tilde{{\omega }_{p}}$ is the detuning of the input field with respect to two times the renormalized first harmonic SP frequency, and similarly for ${\tilde{\delta }}_{b}$. The quantity ${N}_{{\lambda }_{0}^{2}}=(3/2\pi ){({\lambda }_{0}/2)}^{2}/A$ is proportional to the number of structures in a diffraction limited area ${\lambda }_{0}^{2}$, as A is the area of a unit cell in the array. In figure 4, we plot ${| {r}_{b}(\omega )| }^{2}$ as a function of the detuning for different values of the ratio ${\Gamma }^{\prime }/g$. Here we have ignored the renormalized detunings, ${\tilde{\delta }}_{i}\to {\delta }_{i}$ for $i=a,b$, as the structure dimensions can be slightly altered to compensate for these shifts. We also take Q = 100 and Q = 50 for modes a and b, respectively. Note that if the nonlinear interaction between plasmons is negligible ($g\ll \Gamma ^{\prime}$), 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\lt {\Gamma }^{\prime }/2$ to the regime in which $g\gt {\Gamma }^{\prime }/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\sqrt{2}g$ does not depend on the coupling efficiency and represents a robust signature of quantum strong coupling between the SPs modes. We also emphasize that equation (11) is only obtained by solving fully the equations (10), including quantum correlations between the two plasmon modes. Solving the classical limit, in which all quantum operators are replaced with numbers, would produce a Lorentzian spectrum for any value of g, which reinforces the appearance of a mode splitting as a quantum signature.

Zoom In Zoom Out Reset image size

In a similar way, we can calculate the intensity emitted at frequency $2{\omega }_{p}$, when the system is driven at frequency ${\omega }_{p}$ by a classical external field. We find that the SHG signal intensity radiated into the far field is approximately (see appendix D for the derivation)

$$\begin{eqnarray}&&{I}_{2{\omega }_{p}}^{\mathrm{far}}\approx \displaystyle \frac{8\;{g}^{2}}{{\hslash }{\omega }_{p}{{\Gamma }_{a}^{\prime }}^{2}{\Gamma }_{b}^{\prime }}\displaystyle \frac{{[{\sigma }_{a}^{\mathrm{ext}}]}^{2}{\sigma }_{b}^{\mathrm{ext}}}{{A}^{2}}{[{I}_{{\omega }_{p}}^{\mathrm{ext}}]}^{2},\end{eqnarray} \tag{ 12 }$$

where ${\sigma }_{a,b}^{\mathrm{ext}}$ are the extinction cross sections of the two modes. This expression is valid in the undepleted pump approximation, where the converted intensity is a small fraction of the incident. Using the previously quoted parameters for the triangular nanostructure, we find that a 1% conversion efficiency can be observed for the low driving intensity of roughly 10⁸ W m⁻².

While we have presented here a semi-classical calculation, in which the input fields are treated as classical numbers, it would be interesting to find what is the conversion efficiency at the single-photon level. In particular, it would be interesting to see how graphene compares to the state-of-the-art efficiencies of $\sim {10}^{-8}$ in bulk crystals for SHG of just a two-photon input. For this purpose, in the next section we use an approach based on the S-matrix formalism.

In general, 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 [44], which specifically describes the overlap amplitude between a set of monochromatic incoming and outgoing freely propagating photons. Because monochromatic photons form a complete basis, the S-matrix thus contains all information about photon dynamics. In particular, it can be used to determine how a wave packet consisting of a superposition of monochromatic photons (i.e. decomposed into frequency components) interacts with the graphene nanostructure.

A simple example of an S-matrix element consists of the linear reflection amplitude r_b(k) of a single photon of frequency k_b, which interacts with the higher-frequency SP mode (mode b), which we have calculated in the previous section by solving the Heisenberg equations of motion. In the S-matrix language the reflection coefficient corresponds to the matrix element between an incoming photon propagating in one direction (say to the right) and a photon of the same frequency p_b = k_b scattered in the other direction (to the left). More compactly, this relation is formally written as $\langle {p}_{b}^{L}| S| {k}_{b}^{R}\rangle \equiv {r}_{b}(k)\delta (k-p)$, where $\delta (k-p)$ denotes the Dirac delta function. Such an S-matrix element can be calculated by using standard input–output techniques [44, 45], which enable one to relate the outgoing field (after interaction) to the incoming field and internal dynamics of the nanostructure (governed by the Hamiltonian of equation (8)). We assume that the incoming photon is focused at the diffraction limit, $S\sim {\lambda }_{0}^{2}$, and interacts with $N\equiv {N}_{{\lambda }_{0}^{2}}$ structures. In appendix E, we show that the resulting reflection coefficient gives a result of the form of equation (11).

Analogously, we can express the amplitude for the DC process as the S-matrix element between an incoming photon of frequency k_b near $2{\omega }_{p}$ and two outgoing photons of frequencies ${p}_{a},{q}_{a}$ near ${\omega }_{p}$. For simplicity we study the case in which the incoming photon is a superposition of a photon coming from the right and one coming from the left so that we can avoid directional labels. We thus find for an array of N structures

$$\begin{eqnarray}&&\langle {p}_{a},{q}_{a}| S| {k}_{b}\rangle =C\;{r}_{b}(k){r}_{a}(p){r}_{a}(q)\delta (k-p-q),\end{eqnarray} \tag{ 13 }$$

where r_a, r_b are respectively the reflection coefficients for photons in mode a and b, and $C=2Ng/\sqrt{2\pi {\kappa }_{a}^{2}{\kappa }_{b}}$ (see appendix E).

The S-matrix also enables one to calculate the dynamics of an incoming pulse. In particular, assuming a single-photon input wavepacket with a Fourier transform given by f(k), we find that the total DC efficiency is given by ${P}_{\mathrm{DC}}=1/2\int {\rm{d}}p\;{\rm{d}}q\;{| f(p+q){r}_{b}(p+q){r}_{a}(p){r}_{a}(q)| }^{2}$. For a near monochromatic resonant incoming photon, i.e. ${| f(k)| }^{2}\approx \delta (k-2{\omega }_{p})$, the result simplifies to

$$\begin{eqnarray}&&{P}_{\mathrm{DC}}=\displaystyle \frac{16{N}^{2}{\kappa }_{a}^{2}{\kappa }_{b}\;{g}^{2}}{{\Gamma }_{a}^{\prime }{[{\Gamma }_{a}^{\prime }{\Gamma }_{b}^{\prime }+4{g}^{2}]}^{2}},\end{eqnarray} \tag{ 14 }$$

where ${\Gamma }^{\prime }=\Gamma +N\kappa$. The value of the coupling constant that maximizes the probability of conversion is $g=\sqrt{{\Gamma }_{a}^{\prime }{\Gamma }_{b}^{\prime }}/2$, for which we have

$$\begin{eqnarray}&&{P}_{\mathrm{DC}}={N}^{2}{(\displaystyle \frac{{\kappa }_{a}}{{\Gamma }_{a}^{\prime }})}^{2}(\displaystyle \frac{{\kappa }_{b}}{{\Gamma }_{b}^{\prime }}).\end{eqnarray} \tag{ 15 }$$

In general, we expect g to exceed the plasmon linewidth in the graphene nanostructure considered, so that the condition $g=\sqrt{{\Gamma }_{a}^{\prime }{\Gamma }_{b}^{\prime }}/2$ is satisfiable, in contrast to conventional materials with weak nonlinear coefficients. For what concerns the optimal number of nanostructures, we identify two limits, one of low external coupling efficiency in which the array-enhanced external coupling does not overcome the losses i.e. $\kappa \ll \Gamma$, and the opposite case in which $\kappa \gtrsim \Gamma$. In the first limit, which is satisfied for the system parameters presented earlier, the total decay rate $\Gamma ^{\prime}$ is roughly independent of the number of structures N and ${P}_{\mathrm{DC}}^{\mathrm{max}}\approx {N}^{2}{\eta }_{a}^{2}\;{\eta }_{b}$ (we recall again that ${\eta }_{i}={\kappa }_{i}{/\Gamma }_{i}$). It is clear that in this limit the use of an array of nanostructures is an efficient way to increase the conversion (which anyway remains much smaller than 1). For our system parameters, we find that ${P}_{{DC}}^{\mathrm{max}}\approx {10}^{-7}$, which compares favorably with state-of-the-art numbers ∼ 10⁻⁸, a surprising result considering that graphene is not a bulk nonlinear crystal. In the opposite limit of good external coupling we find that ${P}_{\mathrm{DC}}^{\mathrm{max}}={N}^{-1}{\eta }_{a}^{2}\;{\eta }_{b}$. This remarkable result indicates that ultimately, there is a fundamental inequivalence between using many structures to increase the (linear) response, and working to improve the coupling to just a single structure. In particular, in the limit of efficient coupling, the strong nonlinear interaction between plasmons becomes diluted by having multiple structures. Intuitively, this ${N}^{-1}$ scaling can be understood from the complementary process of SHG (whose S-matrix is identical to DC, as shown later). Clearly, in order for two incoming photons to create a second harmonic, they must excite two plasmons in the same structure. However, with many structures, the probability that this occurs (i.e. compared to exciting single plasmons in two different structures) falls like ${N}^{-1}$. We thus argue that the development of techniques [21, 29] to efficiently couple to single structures is of fundamental importance to take maximal advantage of the strong intrinsic nonlinear interactions between graphene plasmons.

It should further be noted that the created photon pairs are frequency-entangled (see equation (13)), 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 $\Gamma ^{\prime}$. This can be seen quantitatively in figure 5(a), where Gaussian single-photon inputs with bandwidth σ are considered, i.e. $f(k)\propto {{\rm{e}}}^{-{(k-{\omega }_{p})}^{2}/4{\sigma }^{2}}$.

Zoom In Zoom Out Reset image size

In SHG two photons with frequencies centered around ${\omega }_{p}$ are (partially) converted in a single photon of frequency $2{\omega }_{p}$. By the time reversal symmetry of the scattering matrix the relation $\langle {p}_{a},{q}_{a}| S| {k}_{b}\rangle \;=\;$ $\langle {k}_{b}| S{| {p}_{a},{q}_{a}\rangle }^{*}$ holds. This implies that in principle, a maximum up-conversion efficiency of ${P}_{\mathrm{SHG}}^{\mathrm{max}}={P}_{\mathrm{DC}}^{\mathrm{max}}$ can be achieved, but only if the two-photon input itself is an entangled state. In figure 5(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 $\sigma \to 0$ and $\sigma \to \infty$. The inability to deterministically up-convert two separate photons (${P}_{\mathrm{SHG}}=1$), even for perfect coupling efficiencies, notably deviates from the semiclassical prediction that perfect conversion can be achieved [37].

We conclude showing that a single graphene nanostructure can generate nonclassical light when irradiated with weak classical light at the lower frequency. We have seen above that in the strong quantum coupling regime, $g\gt {\Gamma }^{\prime }/2$, a mode splitting at the second resonance appears. Physically, this splitting arises because the nonlinear interaction given in the Hamiltonian of equation (8) strongly mixes a single photon $| 0;1\rangle$ in mode $2{\omega }_{p}$ with two photons $| 2;0\rangle$ in mode ${\omega }_{p}$, as shown in figure 6(b). The resulting eigenstates of the Hamiltonian are symmetric and antisymmetric combinations $| 0;1\rangle \pm | 2;0\rangle$ with frequencies $2{\omega }_{p}\pm \sqrt{2}g$. The mode splitting creates an effective nonlinearity: once a single plasmon of frequency ${\omega }_{p}$ enters the system, the absence of a resonant state at $2{\omega }_{p}$ prevents a second plasmon from entering, creating a blockade effect [35]. This is a complementary signature of strong coupling observable in the lower mode. It can be quantified by considering the second-order correlation function of back-scattered photons (for instance left-propagating photons when the system is driven by right-propagating laser light) ${g}^{(2)}(t)\;=\;$ $\langle {a}_{L,\mathrm{out}}^{\dagger }(\tau ){a}_{L,\mathrm{out}}^{\dagger }(\tau +t){a}_{L,\mathrm{out}}(\tau +t){a}_{L,\mathrm{out}}(\tau )\rangle /{\langle {a}_{L,\mathrm{out}}^{\dagger }(\tau ){a}_{L,\mathrm{out}}(\tau )\rangle }^{2}$. The output field itself is related to the input field and plasmon mode by the equation ${a}_{L,\mathrm{out}}={a}_{L,\mathrm{in}}+\sqrt{{\kappa }_{a}/2}a$. However, as the left-going input field is in the vacuum state, the corresponding input operator has no effect. Thus the second-order correlation function can be written directly in terms of the plasmon mode a, ${g}^{(2)}(t)\;=\;$ $\langle {a}^{\dagger }(\tau ){a}^{\dagger }(\tau +t)a(\tau +t)a(\tau )\rangle /{\langle {a}^{\dagger }(\tau )a(\tau )\rangle }^{2}$. For t = 0 this function indicates the relative probability to detect two photons at the same time. Values of ${g}^{(2)}(0)\lt 1$ indicate the presence of nonclassical light. In the limit of weak driving amplitude we find that

$$\begin{eqnarray}&&{g}_{a}^{(2)}(0)=\displaystyle \frac{{{\Gamma }^{\prime }}^{2}(16{g}^{2}+3{{\Gamma }^{\prime }}^{2})}{3{(4{g}^{2}+{{\Gamma }^{\prime }}^{2})}^{2}}.\end{eqnarray} \tag{ 16 }$$

For g = 0 it acquires a value of ${g}_{a}^{(2)}(0)=1$, reflecting the coherent state statistics of the laser, while exhibiting strong anti-bunching $({g}_{a}^{(2)}(0)\lt 1)$ when $g\gtrsim {\Gamma }^{\prime }/2$. It is particularly important that ${g}_{a}^{(2)}(0)$ is independent of the external coupling efficiency $\kappa /\Gamma ^{\prime}$, thus making this effect a robust signature of strong quantum coupling between plasmon modes.

Zoom In Zoom Out Reset image size

We have shown that second-order nonlinear optical interactions between plasmons in graphene nanostructures can be remarkably strong. Signatures of such nonlinearities should be immediately observable in experiments involving arrays of nanostructures, where incident free-space light can undergo frequency mixing at very low input powers via interaction with plasmons.

We further show that single nanostructures should exhibit the capability to generate non-classical states of light, observable even with low coupling efficiencies, which opens up a novel route to quantum optics as compared to the conventional approach of using atom-like emitters. With improved coupling efficiencies to the modes of these nanostructures, it would become possible to realize efficient second-harmonic generation or down-conversion at the level of a few quanta, which would exceed the capabilities of current systems by several orders of magnitude. While we focused on one concrete example consisting of a graphene nanotriangle, our conclusions are quite adaptable. Thus, it would be interesting to explore further the potential of this unique 'nonlinear crystal' in a wide variety of classical and quantum nonlinear optical devices. It would also be interesting to investigate the nonlinear optical response of even smaller structures [46, 47], which is expected to deviate significantly from large-scale graphene due to quantum finite-size effects. Finally, we anticipate that our work will open up the intriguing possibility of a search for new materials that are capable of attaining the quantum nonlinear regime.

The authors thank M Jablan, J D Cox and F H L Koppens for useful discussions. DEC acknowledges support from Fundacio Privada Cellex Barcelona. MTM acknowledges the international PhD-fellowship program 'la Caixa-Severo Ochoa'. All authors acknowledge support from European Project GRASP.

We calculate in this section of the appendix the second-order conductivity of extended graphene. In the main text, we apply this result to a finite-size structure. In principle, it should be noted that the response of a finite-size structure is not spatially homogeneous, and applying the conductivity functions for extended graphene is an approximation which likely gives an error of $\sim {\lambda }_{F}/{\lambda }_{p}$. An exact treatment for finite structures requires the calculation of the response function, which involves intensive numerical computations [46].

As explained in the main text, we use a semiclassical single-band approach describing the dynamics of carriers in graphene, which is nominally a zero-gap semiconductor. However, by techniques such as electrostatic gating [16], one band can become partly filled with carriers. We characterize the carriers within this band by the distribution function ${f}_{{\bf{k}}}({\bf{r}},t)$ which is defined so that

$$\begin{eqnarray}&&{\rm{d}}N={f}_{{\bf{k}}}({\bf{r}},t){{\rm{d}}}^{2}{\bf{k}}\;{{\rm{d}}}^{2}{\bf{r}}\end{eqnarray} \tag{ A.1 }$$

is the number of carriers with positions lying within a surface element ${{\rm{d}}}^{2}{\bf{r}}$ about ${\bf{r}}$ and momenta lying within a momentum space element ${{\rm{d}}}^{2}{\bf{k}}$ about ${\bf{k}}$, at time t. When collisions between carriers are neglected, a conservation equation is satisfied by the function ${f}_{{\bf{k}}}({\bf{r}},t)$ (equation (1) of the main text [30]). In Fourier space it can be written in the form

$$\begin{eqnarray}&&{f}_{{\bf{k}}}({\bf{q}},\omega )=\displaystyle \frac{{\rm{i}}e}{{\hslash }(\omega \mp {v}_{F}\hat{{\bf{k}}}\cdot {\bf{q}})}\int \displaystyle \frac{{{\rm{d}}}^{2}{\bf{p}}}{{(2\pi )}^{2}}{\displaystyle \int }_{-\infty }^{\infty }\displaystyle \frac{{\rm{d}}\nu }{2\pi }{\bf{E}}({\bf{q}}-{\bf{p}},\omega -\nu )\cdot \displaystyle \frac{\partial {f}_{{\bf{k}}}({\bf{p}},\nu )}{\partial {\bf{k}}},\end{eqnarray} \tag{ A.2 }$$

which exhibits a nonlinear character. We assume that the electric field perturbs the equilibrium distribution weakly so that we can solve the equation (1) iteratively, obtaining a perturbation series in ${\bf{E}}$ for the distribution function.

At zero order, ignoring finite-size effects as noted above, we simply replace the distribution function on the RHS with the (zero temperature) Fermi distribution ${f}_{{\bf{k}}}^{(0)}({\bf{r}},t)=\theta ({k}_{F}-k)$, obtaining as solution the first order contribution to the conductivity, which is linear in the electric field:

$$\begin{eqnarray}&&{f}_{{\bf{k}}}^{(1)}({\bf{q}},\omega )=-\displaystyle \frac{{\rm{i}}e\hat{{\bf{k}}}\cdot {\bf{E}}({\bf{q}},\omega )}{{\hslash }(\omega \mp {v}_{F}\hat{{\bf{k}}}\cdot {\bf{q}})}\delta (k-{k}_{F}).\end{eqnarray} \tag{ A.3 }$$

In turn, inserting equation (A.3) into equation (1) of the main text, we get the second-order contribution

$$\begin{eqnarray}{f}_{{\bf{k}}}^{(2)}({\bf{q}},\omega ) & = & \displaystyle \frac{{e}^{2}}{{\hslash }(\omega \mp {v}_{F}\hat{{\bf{k}}}\cdot {\bf{q}})}\\ & & \times \displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\bf{p}}}{{(2\pi )}^{2}}{\displaystyle \int }_{-\infty }^{\infty }\displaystyle \frac{{\rm{d}}\nu }{2\pi }\;{\bf{E}}({\bf{q}}-{\bf{p}},\omega -\nu )\cdot \displaystyle \frac{\partial }{\partial {\bf{k}}}\{\displaystyle \frac{\hat{{\bf{k}}}\cdot {\bf{E}}({\bf{p}},\nu )}{{\hslash }(\nu \mp {v}_{F}\hat{{\bf{k}}}\cdot {\bf{p}})}\delta (k-{k}_{F})\}.\end{eqnarray} \tag{ A.4 }$$

Moreover equation (2) of the main text

$$\begin{eqnarray}&&{\bf{J}}({\bf{r}},t)=-{{eg}}_{v}{g}_{s}\int \displaystyle \frac{{{\rm{d}}}^{2}{\bf{k}}}{{(2\pi )}^{2}}{{\bf{v}}}_{{\bf{k}}}{f}_{{\bf{k}}}({\bf{r}},t),\end{eqnarray} \tag{ A.5 }$$

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 nth-order conductivity, we obtain that

$$\begin{eqnarray}&&{\sigma }_{{ij}}^{(1)}({\bf{q}},\omega )=\displaystyle \frac{{\rm{i}}{e}^{2}{g}_{v}{g}_{s}{v}_{F}}{{\hslash }}\int \displaystyle \frac{{{\rm{d}}}^{2}{\bf{k}}}{{(2\pi )}^{2}}\displaystyle \frac{{k}_{i}{k}_{j}}{{k}^{2}(\omega \mp {v}_{F}\hat{{\bf{k}}}\cdot {\bf{q}})}\delta (k-{k}_{F}),\end{eqnarray} \tag{ A.6 }$$

and

$$\begin{eqnarray}&&{\sigma }_{{ilm}}^{(2)}({\bf{q}},\omega ;{\bf{q}}-{\bf{p}},\omega -\nu ,{\bf{p}},\nu )=\displaystyle \frac{{e}^{3}{g}_{v}{g}_{s}{v}_{F}}{{{\hslash }}^{2}}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\bf{k}}}{{(2\pi )}^{2}}\;\displaystyle \frac{1}{{k}^{2}(\omega \mp {v}_{F}\hat{{\bf{k}}}\cdot {\bf{q}})(\nu \mp {v}_{F}\hat{{\bf{k}}}\cdot {\bf{p}})}\\ &&\times [{\delta }_{{il}}{k}_{m}-\displaystyle \frac{\omega }{{k}^{2}(\omega \mp {v}_{F}\hat{{\bf{k}}}\cdot {\bf{q}})}{k}_{i}{k}_{l}{k}_{m}\pm \displaystyle \frac{{v}_{F}}{k(\omega \mp {v}_{F}\hat{{\bf{k}}}\cdot {\bf{q}})}{k}_{i}{q}_{l}{k}_{m}]\delta (k-{k}_{F}),\end{eqnarray} \tag{ A.7 }$$

where the integration over ${\bf{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/\omega \ll 1$), by expanding the denominators in q and p. The dominant contribution to the linear conductivity is the zero-order one, giving rise to the well-known local Drude conductivity of graphene displayed in the equation (3) of the main text.

For the second-order conductivity, it can be easily proven that the zero-order expansion in q and p of equation (A.7) gives a vanishing contribution when the integral is performed, so that the dominant term is the first-order expansion in q and p, which corresponds to a nonlocal contribution.

The results are given by

$$\begin{eqnarray}{\sigma }_{{ilm}}^{(2)}({\bf{q}},\omega ;{\bf{q}}-{\bf{p}},\omega -\nu ,{\bf{p}},\nu ) & = & \pm \displaystyle \frac{{e}^{3}{g}_{v}{g}_{s}{v}_{F}^{2}}{16\pi {{\hslash }}^{2}}[\displaystyle \frac{2}{{\omega }^{2}\nu }({\delta }_{{il}}{q}_{m}-{\delta }_{{lm}}{q}_{i}+{\delta }_{{im}}{q}_{l})\\ & & +\displaystyle \frac{1}{\omega {\nu }^{2}}(3{\delta }_{{il}}{p}_{m}-{\delta }_{{lm}}{p}_{i}-{\delta }_{{im}}{p}_{l})].\end{eqnarray} \tag{ A.8 }$$

With the formula for the conductivity above, the expression

$$\begin{eqnarray}{J}_{i}^{(2)}({\bf{q}},\omega ) & = & \displaystyle \int \displaystyle \frac{{{\rm{d}}}^{2}{\bf{p}}}{{(2\pi )}^{2}}{\displaystyle \int }_{-\infty }^{\infty }\displaystyle \frac{{\rm{d}}\nu }{2\pi }\;{\sigma }_{{ilm}}^{(2)}({\bf{q}},\omega ;{\bf{q}}-{\bf{p}},\omega -\nu ,{\bf{p}},\nu )\\ & & \times {E}_{l}({\bf{q}}-{\bf{p}},\omega -\nu ){E}_{m}({\bf{p}},\nu ),\end{eqnarray} \tag{ A.9 }$$

can be Fourier transformed back to the real space, resulting finally in equations (4) and (5) of the main text.

In the limit $D{/\lambda }_{0}\ll 1$, where D is the linear dimension of the graphene structure and ${\lambda }_{0}$ is the incident radiation wavelength, the electrostatic approximation can be assumed. Thus, the total energy present in the structure is given by

$$\begin{eqnarray}&&H=\displaystyle \frac{1}{2}\displaystyle \sum _{{\omega }_{i}}{\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{{\rho }^{{\omega }_{i}}}^{*}({\bf{r}}){\phi }^{{\omega }_{i}}({\bf{r}}),\end{eqnarray} \tag{ B.1 }$$

where ρ is the charge density and ϕ the electrostatic potential. Implementing the continuity equation, the relation between the potential and the electric field, and writing explicitly the two frequencies contributions, we get that the total energy is given by

$$\begin{eqnarray}&&H=\displaystyle \frac{1}{2{\rm{i}}{\omega }_{p}}{\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{{J}_{i}^{{\omega }_{p}}}^{*}({\bf{r}}){E}_{i}^{{\omega }_{p}}({\bf{r}})+\displaystyle \frac{1}{4{\rm{i}}{\omega }_{p}}{\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{{J}_{i}^{2{\omega }_{p}}}^{*}({\bf{r}}){E}_{i}^{2{\omega }_{p}}({\bf{r}}).\end{eqnarray} \tag{ B.2 }$$

The currents can be expressed in terms of the electric fields: ${J}_{i}^{{\omega }_{p}}={\sigma }^{(1)}({\omega }_{p}){E}_{i}^{{\omega }_{p}}$, while ${J}_{i}^{2{\omega }_{p}}$ is given by equation (4) of the main text.

At this point, we can impose the quantization condition to the first mode:

$$\begin{eqnarray}&&\displaystyle \frac{{\sigma }^{(1)}({\omega }_{p})}{2{\rm{i}}{\omega }_{p}}{\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{| {E}_{i}^{{\omega }_{p}}({\bf{r}})| }^{2}={\hslash }{\omega }_{p}\;{a}^{\dagger }a,\end{eqnarray} \tag{ B.3 }$$

which can be enforced with the substitution ${E}_{i}^{{\omega }_{p}}({\bf{r}})\to {\tilde{E}}_{i}^{{\omega }_{p}}({\bf{r}})a={E}_{0}^{{\omega }_{p}}\;{f}_{i}^{{\omega }_{p}}({\bf{r}})a$. Here, ${f}^{{\omega }_{p}}({\bf{r}})$ is a vectorial function which describes the geometry of the mode and normalized such that $\mathrm{max}| {f}^{{\omega }_{p}}({\bf{r}})| =1$, ${E}_{0}^{{\omega }_{p}}={({\hslash }{\omega }_{p}{q}_{p}/{\epsilon }_{0}S\mu )}^{1/2}$ is the maximum single-photon electric field amplitude, and $\mu ={S}_{\mathrm{eff}}^{{\omega }_{p}}/S$, where ${S}_{\mathrm{eff}}^{{\omega }_{p}}={\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{| {f}_{i}^{{\omega }_{p}}({\bf{r}})| }^{2}$ is the ratio between the effective mode area and the physical area of the structure.

Now we consider the mode at frequency $2{\omega }_{p}$. We substitute the current shown in equation (4) into equation (B.2), and similarly define a single-photon electric field amplitude and annihilation operator b for this mode. This procedure yields a non-interacting term in the Hamiltonian, $2{\hslash }{\omega }_{p}{b}^{\dagger }b$, and an interacting term

$$\begin{eqnarray}&&{H}_{\mathrm{int}}=\displaystyle \frac{1}{4{\rm{i}}{\omega }_{p}}\;{\sigma }_{{ijkl}}^{(2)}(2{\omega }_{p};{\omega }_{p}){\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{[{\tilde{E}}_{i}^{2{\omega }_{p}}({\bf{r}})]}^{*}{\tilde{E}}_{j}^{{\omega }_{p}}({\bf{r}}){\nabla }_{k}{\tilde{E}}_{l}^{{\omega }_{p}}({\bf{r}}){a}^{2}{b}^{\dagger }+{\rm{h}}.{\rm{c}}.\end{eqnarray} \tag{ B.4 }$$

Note that using the freedom in the definition of the phases of ${\tilde{E}}_{i}^{{\omega }_{p}}$ and ${\tilde{E}}_{i}^{2{\omega }_{p}}$ we can make the expression in front of ${a}^{2}{b}^{\dagger }$, ${\hslash }g$, to be real, obtaining in this way equation (9) of the main text.

If we express ${\tilde{E}}_{l}^{{\omega }_{p}}({\bf{r}})$ as ${E}_{0}^{{\omega }_{p}}\;{f}_{i}^{{\omega }_{p}}({\bf{r}})$, in order to separate the geometric part of the problem, we get that the ratio between g and ${\omega }_{p}$ is equal to

$$\begin{eqnarray}\displaystyle \frac{g}{{\omega }_{p}} & = & \displaystyle \frac{1}{4{\rm{i}}{\hslash }{\omega }_{p}^{2}}\;{\sigma }_{{ijkl}}^{(2)}(2{\omega }_{p};{\omega }_{p}){[{E}_{0}^{{\omega }_{p}}]}^{2}\;{E}_{0}^{2{\omega }_{p}}{\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{{f}_{i}^{2{\omega }_{p}}}^{*}({\bf{r}}){f}_{j}^{{\omega }_{p}}({\bf{r}}){\nabla }_{k}{f}_{l}^{{\omega }_{p}}({\bf{r}})\\ & = & \mp \displaystyle \frac{1}{4{\rm{i}}{\hslash }{\omega }_{p}^{2}}\;\displaystyle \frac{{\rm{i}}{e}^{3}{v}_{F}^{2}}{8\pi {{\hslash }}^{2}{\omega }_{p}^{3}}\;\displaystyle \frac{2\pi {{\hslash }}^{3}{\omega }_{p}^{3}}{{e}^{2}| {E}_{F}| {S}_{\mathrm{eff}}^{{\omega }_{p}}}{(\displaystyle \frac{16\pi {{\hslash }}^{3}{\omega }_{p}^{3}}{{e}^{2}| {E}_{F}| {S}_{\mathrm{eff}}^{2{\omega }_{p}}})}^{1/2}\;{\Theta }_{{ijkl}}{\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{{f}_{i}^{2{\omega }_{p}}}^{*}({\bf{r}}){f}_{j}^{{\omega }_{p}}({\bf{r}}){\nabla }_{k}{f}_{l}^{{\omega }_{p}}({\bf{r}}),\end{eqnarray} \tag{ B.5 }$$

where we have called the tensorial part of the second-order conductivity ${\Theta }_{{ijkl}}=(5{\delta }_{{ij}}{\delta }_{{kl}}-3{\delta }_{{ik}}{\delta }_{{jl}}+{\delta }_{{il}}{\delta }_{{jk}})$. Now we make the lengths in the integral dimensionless, introducing $\hat{{\bf{r}}}={\bf{r}}/D$ and $\hat{\nabla }=D\nabla$. Furthermore, we introduce the dimensionless parameter ${\xi }_{1}={D}^{3}{/S}_{\mathrm{eff}}^{{\omega }_{p}}{[{S}_{\mathrm{eff}}^{2{\omega }_{p}}]}^{1/2}$. From equation (B.5), we obtain the final expression

$$\begin{eqnarray}&&\displaystyle \frac{g}{{\omega }_{p}}=\mp \displaystyle \frac{{v}_{F}^{2}{{\hslash }}^{2}}{4| {E}_{F}| {D}^{2}}{(\displaystyle \frac{\pi }{{\hslash }{\omega }_{p}| {E}_{F}| })}^{1/2}\;{\xi }_{1}\;{\Theta }_{{ijkl}}{\displaystyle \int }_{S}{{\rm{d}}}^{2}\hat{{\bf{r}}}\;{f}_{i}^{2{\omega }_{p}}(\hat{{\bf{r}}}){f}_{j}^{{\omega }_{p}}(\hat{{\bf{r}}}){\hat{\nabla }}_{k}{f}_{l}^{{\omega }_{p}}(\tilde{{\bf{r}}}).\end{eqnarray} \tag{ B.6 }$$

Note that now the integral in the last expression is fully dimensionless and depends only on the geometry of the two modes.

In a general way, the plasmon frequency ${\omega }_{p}$ is related to the dimension of the structure and the doping by the relation [39]

$$\begin{eqnarray}&&{\omega }_{p}={\xi }_{2}{(\displaystyle \frac{2\alpha c}{{\hslash }}\displaystyle \frac{{E}_{F}}{D})}^{1/2},\end{eqnarray} \tag{ B.7 }$$

where ${\xi }_{2}$ is a factor of order one that depends only on the shape of the structure considered, α is the fine-structure constant, and c is the speed of light. Using this relation in equation (B.6), and considering that $| {E}_{F}| ={\hslash }{v}_{F}{k}_{F}$, we finally find

$$\begin{eqnarray}&&\displaystyle \frac{g}{{\omega }_{p}}=\displaystyle \frac{\beta }{{({k}_{F}D)}^{7/4}},\end{eqnarray} \tag{ B.8 }$$

where we have collected in β all the geometric dimensionless factors

$$\begin{eqnarray}&&\beta =\displaystyle \frac{{\xi }_{1}}{4\sqrt{{\xi }_{2}}}{(\displaystyle \frac{{\pi }^{2}{v}_{F}}{2\alpha c})}^{1/4}\;{\Theta }_{{ijkl}}{\displaystyle \int }_{S}{{\rm{d}}}^{2}\hat{{\bf{r}}}\;{f}_{i}^{2{\omega }_{p}}(\hat{{\bf{r}}}){f}_{j}^{{\omega }_{p}}(\hat{{\bf{r}}}){\hat{\nabla }}_{k}{f}_{l}^{{\omega }_{p}}(\tilde{{\bf{r}}}).\end{eqnarray} \tag{ B.9 }$$

Note that, while the scaling of the nonlinearity as ${({k}_{F}D)}^{-7/4}$ is a completely general result, the individual terms appearing in β are defined in a somewhat arbitrary way. For instance, choosing D to be the length of the short side of a triangle instead of that of the long one, as done for the numerical example in the paper, affects the numerical values of ${\xi }_{1}$ and ${\xi }_{2}$.

There are two ways of obtaining the radiative decay rate of the graphene triangle, one through the extinction cross section, and the other directly through the dipole moment.

When using the first method, we assume that the polarizability of the graphene triangle can be expressed as a Lorentzian line shape [48]

$$\begin{eqnarray}&&\alpha (\omega )=\displaystyle \frac{6\pi {\epsilon }_{0}{c}^{3}{\kappa }_{a}}{{\omega }_{p}^{2}}\displaystyle \frac{1}{{\omega }_{p}^{2}-{\omega }^{2}-{\rm{i}}{\Gamma }^{\prime }{\omega }^{3}{/\omega }_{p}^{2}},\end{eqnarray} \tag{ C.1 }$$

with ${\omega }_{p}$ being the plasmon frequency, $\Gamma ^{\prime}$ the total decay rate, and ${\kappa }_{a}$ the radiative contribution to $\Gamma ^{\prime}$. On the other hand, the extinction cross section is directly related to the polarizability by ${\sigma }^{\mathrm{ext}}(\omega )=(\omega /c)\mathrm{Im}\{\alpha (\omega )\}{/\epsilon }_{0}$ [49]. The combination of the two former equations at the plasmon frequency results in the final expression for the radiative decay

$$\begin{eqnarray}&&{\hslash }{\kappa }_{a}({\omega }_{p})=\displaystyle \frac{{\omega }_{p}^{2}{\sigma }^{\mathrm{ext}}({\omega }_{p})}{6\pi {c}^{2}}{\hslash }{\Gamma }^{\prime }.\end{eqnarray} \tag{ C.2 }$$

Thus, we can get the value of the radiative decay by obtaining numerically the extinction cross section at the plasmon frequency.

In the second method, the radiative decay is calculated directly from the knowledge of the induced dipole moment ${\bf{p}}$ of the plasmon modes of the triangle, which are numerically computed from ${\mathrm{COMSOL}}^{\circledR }$. Indeed, the radiative decay of a dipole is given by [50]

$$\begin{eqnarray}&&{\hslash }{\kappa }_{a}({\omega }_{p})=\displaystyle \frac{{\omega }_{p}^{3}}{3\pi {\varepsilon }_{0}{c}^{3}}{| {{\bf{p}}}_{a}| }^{2},\end{eqnarray} \tag{ C.3 }$$

where the dipole moment ${{\bf{p}}}_{a}={\displaystyle \int }_{S}{{\rm{d}}}^{2}{\bf{r}}\;{\rho }_{a}\;{\bf{r}}$ of the mode can be related to the single-plasmon electric field by the continuity equation ${\rho }_{a}=(-{\rm{i}}{\sigma }^{(1)}/\omega ){\nabla }_{\parallel }\cdot \tilde{{\bf{E}}}$.

Assuming that the lattice is large enough that edge effects can be ignored, the dipoles will respond identically, so that in equation (10) of the main paper we get ${a}_{j}=a$ and ${b}_{j}=b$, and the interaction between them reduces to finding the sum $G={\displaystyle \sum }_{j}{G}_{{ij}}$. For a hexagonal lattice this sum consists of a real part that can be approximated as $\mathrm{Re}[{G}^{\omega }]\approx (1/4\pi {\epsilon }_{0})5.2{/l}^{3}$ and an imaginary contribution whose exact expression is $\mathrm{Im}[{G}^{\omega }]=(1/4\pi {\epsilon }_{0})[2\pi \omega /{cA}-2{(\omega /c)}^{3}/3]$ [43], where $A={l}^{2}\sqrt{3}/2$ is the area of the unit cell. We now focus on the case of SHG, where the fundamental mode is driven by a strong external input field with frequency around ${\omega }_{p}$, while the higher mode is undriven (${E}_{2\omega }^{\mathrm{ext}}=0$). In the strong field limit, we can replace the operators with numbers. We can solve for the dipole moments in the frequency domain,

$$\begin{eqnarray}b & = & {[2\omega -2{\omega }_{p}+{\rm{i}}{\Gamma }_{b}/2+\displaystyle \frac{{p}_{b}^{2}}{{\hslash }}\;{G}^{2{\omega }_{p}}]}^{-1}\;{{ga}}^{2}\\ & = & -\displaystyle \frac{{\hslash }}{{p}_{b}^{2}}\;{\tilde{\alpha }}_{b}(2\omega ){{ga}}^{2}\end{eqnarray} \tag{ D.1 }$$

Ignoring the depletion of a due to b (an approximation valid when $a\gg b$):

$$\begin{eqnarray}a & = & {[\omega -{\omega }_{p}+{\rm{i}}{\Gamma }_{a}/2+\displaystyle \frac{{p}_{a}^{2}}{{\hslash }}{G}^{{\omega }_{p}}]}^{-1}\;\displaystyle \frac{{p}_{a}}{{\hslash }}\;{E}_{\omega }^{\mathrm{ext}}\\ & = & -\displaystyle \frac{{E}_{\omega }^{\mathrm{ext}}}{{p}_{a}}\;{\tilde{\alpha }}_{a}(\omega ).\end{eqnarray} \tag{ D.2 }$$

In the previous two equations we have defined ${\tilde{\alpha }}_{a}(\omega )$ as

$$\begin{eqnarray}&&{\tilde{\alpha }}_{a}(\omega )=\displaystyle \frac{-{p}_{a}^{2}}{{\hslash }(\omega -{\omega }_{p}+{\rm{i}}{\Gamma }_{a}/2+\displaystyle \frac{{p}_{a}^{2}}{{\hslash }}{G}^{{\omega }_{p}})}=\displaystyle \frac{-{p}_{a}^{2}}{{\hslash }({\tilde{\delta }}_{a}(\omega )+{\rm{i}}{\tilde{\Gamma }}_{a}/2)},\end{eqnarray} \tag{ D.3 }$$

and similarly for ${\tilde{\alpha }}_{b}(\omega )$. They correspond to the polarizabilities in the proximity of the resonances, modified by the interaction between the structures of the array. Here, ${\tilde{\delta }}_{a}(\omega )$ and ${\tilde{\Gamma }}_{a}$ are the detuning and the linewidth of the plasmonic mode renormalized respectively by the real and the imaginary part of the Green function. Inserting equation (D.2) in (D.1) and multiplying by the dipole moment of a single plasmon in the second mode, we find

$$\begin{eqnarray}&&{p}_{2{\omega }_{p}}(\omega )=-\displaystyle \frac{{\hslash }g}{{p}_{a}^{2}{p}_{b}}\;{\tilde{\alpha }}_{a}^{2}(\omega ){\tilde{\alpha }}_{b}(\omega ){[{E}_{\omega }^{\mathrm{ext}}]}^{2}.\end{eqnarray} \tag{ D.4 }$$

For the lattice constant of interest, the far-field is emitted only in the perpendicular direction [43], and the far-field intensity at frequency $2{\omega }_{p}$ under driving at frequency ${\omega }_{p}$ can be directly calculated as

$$\begin{eqnarray}{I}_{2{\omega }_{p}}^{\mathrm{far}} & = & {| \displaystyle \frac{{\omega }_{p}}{c{\epsilon }_{0}A}\;{p}_{2{\omega }_{p}}({\omega }_{p})| }^{2}\displaystyle \frac{{\epsilon }_{0}c}{2}\\ & = & \displaystyle \frac{2{{\hslash }}^{2}{\omega }_{p}^{2}\;{g}^{2}}{{c}^{3}{\epsilon }_{0}^{3}{A}^{2}{p}_{a}^{4}{p}_{b}^{2}}{[{I}_{{\omega }_{p}}^{\mathrm{ext}}]}^{2}\;{| {\tilde{\alpha }}_{a}(\omega )| }^{4}{| {\tilde{\alpha }}_{b}(\omega )| }^{2}\\ & = & \displaystyle \frac{{g}^{2}\;{{\Gamma }_{a}^{\prime }}^{2}\;{\Gamma }_{b}^{\prime }}{8{\hslash }{\omega }_{p}}\displaystyle \frac{{[{\sigma }_{a}^{\mathrm{ext}}]}^{2}{\sigma }_{b}^{\mathrm{ext}}}{{A}^{2}}{[{I}_{{\omega }_{p}}^{\mathrm{ext}}]}^{2}{[{\tilde{\Gamma }}_{a}^{2}/4+{\tilde{\delta }}_{a}^{2}]}^{-2}{[{\tilde{\Gamma }}_{b}^{2}/4+{\tilde{\delta }}_{b}^{2}]}^{-1},\end{eqnarray} \tag{ D.5 }$$

where ${I}_{{\omega }_{p}}^{\mathrm{ext}}=c{\epsilon }_{0}{[{E}_{{\omega }_{p}}^{\mathrm{ext}}]}^{2}/2$.

Now we can calculate the efficiency of conversion using the values obtained for the triangular nanoisland with quality factor Q = 100, i.e. ${\kappa }_{a}\approx 2\times {10}^{-7}{\omega }_{p}$, ${\kappa }_{b}\approx 5.4\times {10}^{-8}{\omega }_{p}$, and ${\Gamma }_{a}^{\prime }={\Gamma }_{b}^{\prime }={10}^{-2}{\omega }_{p}$, $g=1.25\times {10}^{-2}{\omega }_{p}$. We assume that the lattice period is l = 50 nm. The wave vector of the driving light at frequency ${\omega }_{p}$ is $k\approx 1$ $\mu {{\rm{m}}}^{-1}$. We first calculate the effect of the array on the frequency and the linewidth. The frequency of the first harmonic is redshifted by $6\times {10}^{-3}{\omega }_{p}$, while for the second one the redshift is $2\times {10}^{-4}{\omega }_{p}$. Furthermore, the linewidths ${\Gamma }_{a}^{\prime }$ and ${\Gamma }_{b}^{\prime }$ are increased by around 17% and 1% respectively, thus we can neglect the effects on the second harmonic. We finally find

$$\begin{eqnarray}&&{I}_{2{\omega }_{p}}^{\mathrm{far}}\sim \mathrm{}{10}^{-10\;}{{\rm{W}}}^{-1\;}{{\rm{m}}}^{2}{[{I}_{{\omega }_{p}}^{\mathrm{ext}}]}^{2},\end{eqnarray} \tag{ D.6 }$$

an expression valid only when ${I}_{2{\omega }_{p}}^{\mathrm{far}}\ll {I}_{{\omega }_{p}}^{\mathrm{ext}}$. This calculation indicates that one can observe 1% of intensity conversion when the driving field has an intensity of about 10⁸ ${\rm{W}}\;{{\rm{m}}}^{-2}$.

In this section, we present the combined S-matrix and input–output formalism of [44], which is very helpful to study few-photon scattering amplitudes. We show the equations for a single structure, explaining how to generalize them to the case of the array at the end of the section. We model the incident radiation as a one-dimensional bidirectional continuum, with the two SP modes coupled to it at rates ${\kappa }_{a}$ and ${\kappa }_{b}$. The Heisenberg equations of motion for the internal mode operators, written in terms of the input operators, are given by [45]

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}a}{{\rm{d}}t} & = & -{\rm{i}}({\omega }_{p}-{\rm{i}}{\Gamma }_{a}^{\prime }/2)a-2{\rm{i}}{{gba}}^{\dagger }+{F}_{\mathrm{in}}^{a}+{F}_{\mathrm{loss}}^{a},\\ \displaystyle \frac{{\rm{d}}b}{{\rm{d}}t} & = & -{\rm{i}}(2{\omega }_{p}-{\rm{i}}{\Gamma }_{b}^{\prime }/2)b-{\rm{i}}{{ga}}^{2}+{F}_{\mathrm{in}}^{b}+{F}_{\mathrm{loss}}^{b}.\end{eqnarray} \tag{ E.1 }$$

Here ${F}_{\mathrm{in}}^{a}=-{\rm{i}}\sqrt{{\kappa }_{a}/2}{\displaystyle \sum }_{j=L,R}{a}_{\mathrm{in}}^{j}$ describes coupling from input channels, ${F}_{\mathrm{loss}}^{a}=-{\rm{i}}\sqrt{\gamma }{a}_{\mathrm{loss}}$ describes coupling to loss channels. 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

$$\begin{eqnarray}&&{a}_{\mathrm{out}}^{j}(t)={a}_{\mathrm{in}}^{j}(t)-{\rm{i}}\sqrt{\displaystyle \frac{{\kappa }_{a}}{2}}a(t),\end{eqnarray} \tag{ E.2 }$$

$$\begin{eqnarray}&&{b}_{\mathrm{out}}^{j}(t)={b}_{\mathrm{in}}^{j}(t)-{\rm{i}}\sqrt{\displaystyle \frac{{\kappa }_{b}}{2}}b(t),\end{eqnarray} \tag{ E.3 }$$

with $j=L,R$. These 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 expressed in term of the S-matrix elements

$$\begin{eqnarray}&&\langle {p}_{b}^{L}| S| {k}_{b}^{R}\rangle =\langle 0| {b}_{\mathrm{out}}^{L}(p){{b}_{\mathrm{in}}^{R}}^{\dagger }(k)| 0\rangle .\end{eqnarray} \tag{ E.4 }$$

Using the Fourier transform of equation (E.3)

$$\begin{eqnarray}&&{b}_{\mathrm{out}}^{L}(p)={b}_{\mathrm{in}}^{L}(p)-{\rm{i}}\sqrt{\displaystyle \frac{{\kappa }_{b}}{2}}b(p)\end{eqnarray} \tag{ E.5 }$$

to express the output operator in equation (E.4), we get that

$$\begin{eqnarray}&&\langle {p}_{b}^{L}| S| {k}_{b}^{R}\rangle =-{\rm{i}}\sqrt{\displaystyle \frac{{\kappa }_{b}}{2}}\langle 0| b(p)| {k}_{b}^{R}\rangle =-{\rm{i}}\sqrt{\displaystyle \frac{{\kappa }_{b}}{4\pi }}\int {\rm{d}}p\langle 0| b(t)| {k}_{b}^{R}\rangle \;{{\rm{e}}}^{{\rm{i}}{pt}}.\end{eqnarray} \tag{ E.6 }$$

To determine this matrix element, we use the equation of motion (E.2) for the operator b. Using the commutation relations between the different input operators, $[{b}_{\mathrm{in}}^{R}(k),{{b}_{\mathrm{in}}^{R}}^{\dagger }(k^{\prime} )]=\delta (k-k^{\prime} )$, $[{b}_{\mathrm{in}}^{R}(k),{{b}_{\mathrm{in}}^{L}}^{\dagger }(k^{\prime} )]=0$, etc, we see that only the term with ${b}_{\mathrm{in}}^{R}$ remains:

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle 0| b(t)| {k}_{b}^{R}\rangle & = & -{\rm{i}}(2{\omega }_{p}-{\rm{i}}{\Gamma }_{b}^{\prime }/2)\langle 0| b(t)| {k}_{b}^{R}\rangle \\ & & -{\rm{i}}g\langle 0| {a}^{2}(t)| {k}_{b}^{R}\rangle -{\rm{i}}\sqrt{\displaystyle \frac{{\kappa }_{b}}{2}}\;\langle 0| {b}_{\mathrm{in}}^{R}(t){{b}_{\mathrm{in}}^{R}}^{\dagger }(k)| 0\rangle ,\end{eqnarray} \tag{ E.7 }$$

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², one finds the equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\langle 0| {a}^{2}(t)| {k}_{b}^{R}\rangle =-2{\rm{i}}({\omega }_{p}-{\rm{i}}{\Gamma }_{a}^{\prime })\langle 0| {a}^{2}(t)| {k}_{b}^{R}\rangle -2{\rm{i}}g\langle 0| b(t)| {k}_{b}^{R}\rangle .\end{eqnarray} \tag{ E.8 }$$

The last two equations constitute a closed system of differential equations in $\langle 0| b(t)| {k}_{b}^{R}\rangle$ and $\langle 0| {a}^{2}(t)| {k}_{b}^{R}\rangle$, which can be solved by Fourier transformation. Inserting the results for the former element in equation (E.6), we finally get that

$$\begin{eqnarray}&&\langle {p}_{b}^{L}| S| {k}_{b}^{R}\rangle ={r}_{b}(k)\delta (p-k),\end{eqnarray} \tag{ E.9 }$$

with r_b(k) given by

$$\begin{eqnarray}&&{r}_{b}(k)=-\displaystyle \frac{{\rm{i}}{\kappa }_{b}}{2}\displaystyle \frac{{\delta }_{k}+{\rm{i}}{\Gamma }_{a}^{\prime }}{[{\delta }_{k}+{\rm{i}}{\Gamma }_{a}^{\prime }][{\delta }_{k}+{\rm{i}}{\Gamma }_{b}^{\prime }/2]-2{g}^{2}}.\end{eqnarray} \tag{ E.10 }$$

With a very similar procedure, the S-matrix elements can be calculated 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 impinging symmetrically from the modes $L,R$. Provided this condition, we can introduce the symmetric mode ${a}_{\mathrm{in}}=1/\sqrt{2}({a}_{\mathrm{in}}^{R}+{a}_{\mathrm{in}}^{L})$ (and similarly for b_in and the output modes) in the equations (E.1)–(E.3), and consider a vacuum input in the anti-symmetric mode. The starting point of the calculation is again the expression of the S-matrix element in term of input and output operators:

$$\begin{eqnarray}&&\langle {p}_{a},{q}_{a}| S| {k}_{b}\rangle =\langle 0| {a}_{\mathrm{out}}(p){a}_{\mathrm{out}}(q){b}_{\mathrm{in}}^{\dagger }(k)| 0\rangle .\end{eqnarray} \tag{ E.11 }$$

For what concerns the generalization to an array of N structures, one has only to replace equations (E.2)–(E.3) with

$$\begin{eqnarray}&&{a}_{\mathrm{out}}(t)={a}_{\mathrm{in}}(t)-{\rm{i}}\sqrt{N{\kappa }_{a}}A(t).\end{eqnarray} \tag{ E.12 }$$

$$\begin{eqnarray}&&{b}_{\mathrm{out}}(t)={b}_{\mathrm{in}}(t)-{\rm{i}}\sqrt{N{\kappa }_{b}}B(t),\end{eqnarray} \tag{ E.13 }$$

where $A(t)={N}^{-1/2}{\displaystyle \sum }_{i}{a}_{i}(t)$ and $B(t)={N}^{-1/2}{\displaystyle \sum }_{i}{b}_{i}(t)$ are collective modes, and solve similar differential equations.