Nonlinear atom-plasmon interactions enabled by nanostructured graphene

TL;DR

This paper proposes using nonlinear optical responses of nanostructured graphene to enable resonant coupling with quantum emitters in the near-infrared, overcoming frequency limitations of graphene plasmons.

Contribution

It introduces a novel scheme utilizing nonlinear near-fields from nanostructured graphene to coherently control quantum states in atomic systems.

Findings

Nonlinear near-fields from graphene nanodisks can excite quantum states.

Third harmonic of plasmon resonance can be tuned to atomic transitions.

The scheme avoids strong spontaneous emission enhancement.

Abstract

Electrically tunable graphene plasmons are anticipated to enable strong light-matter interactions with resonant quantum emitters. However, plasmon resonances in graphene are typically limited to infrared frequencies, below those of optical excitations in robust quantum light sources and many biologically interesting molecules. Here we propose to utilize near fields generated by the plasmon-assisted nonlinear optical response of nanostructured graphene to resonantly couple with proximal quantum emitters operating in the near-infrared. We show that the nonlinear near-field produced by a graphene nanodisk can strongly excite and coherently control quantum states in two- and three-level atomic systems when the third harmonic of its plasmon resonance is tuned to a particular electronic transition. In the present scheme, emitter and plasmon resonances are nondegenerate, circumventing strong enhancement of spontaneous emission. We envision potential applications for the proposed nonlinear plasmonic coupling scheme in sensing and temporal quantum control.

Figures

Equations

$$\frac{\partial \rho }{\partial t}=-\frac{\text{i}}{\hslash }[\mathcal{H},\rho ]+\mathcal{L}[\rho ],$$

$$\mathcal{H}=\hslash \sum _{j=1,2}{\epsilon }_j|j⟩⟨j|-{\vec{\mu }}_{12}\cdot \mathbf{E}(\mathbf{r},t)\left(|1⟩⟨2|+|2⟩⟨1|\right),$$

$$\frac{\partial {\rho }_{11}}{\partial t}$$

$$\frac{\partial {\tilde{\rho }}_{21}}{\partial t}$$

$$\Omega =$$

$$\times \frac{{\vec{\xi }}_{m^{\prime }}\cdot {\mathbf{E}}_{\mathrm{ext}}}{1-{\eta }_{\omega }^{(1)}/{\eta }_{m^{\prime }}}\frac{{\vec{\xi }}_{m^{\prime \prime }}\cdot {\mathbf{E}}_{\mathrm{ext}}}{1-{\eta }_{\omega }^{(1)}/{\eta }_{m^{\prime \prime }}}\frac{{\vec{\xi }}_{m^{\prime \prime \prime }}\cdot {\mathbf{E}}_{\mathrm{ext}}}{1-{\eta }_{\omega }^{(1)}/{\eta }_{m^{\prime \prime \prime }}}$$

$$G=\frac{\text{i}{\sigma }_{3\omega }^{(1)}}{3\hslash \omega ({ϵ}_{3\omega }^{\mathrm{a}}{)}^2{D}^4}\sum _m\frac{{\left({\vec{\mu }}_{12}\cdot {\mathbf{e}}_m\right)}^2}{1-{\eta }_{3\omega }^{(1)}/{\eta }_m}$$

$${\mathbf{e}}_m(\mathbf{r})=\int d^2\mathbf{R}{\rho }_m(\mathbf{R})\frac{\mathbf{r}-\mathbf{R}}{{\left|\mathbf{r}-\mathbf{R}\right|}^3}$$

$${\rho }_{\mathrm{ind}}^{ns}(\mathbf{R},\omega )=-\frac{\text{i}}{s\omega }{\nabla }_{\mathbf{R}}\cdot {\mathbf{j}}^{ns}(\mathbf{R},\omega ).$$

$${\mathbf{j}}^{ns}(\mathbf{R},\omega )=f_{\mathbf{R}}\left[{\mathbf{j}}_{\mathrm{NL}}^{ns}(\mathbf{R},\omega )+{\sigma }_{s\omega }^{(1)}{\mathbf{E}}^{ns}(\mathbf{R},\omega )\right],$$

$${\mathbf{E}}^{ns}(\mathbf{R},\omega )=-{\nabla }_{\mathbf{R}}$$

$$+\frac{1}{{ϵ}_{s\omega }^{\mathrm{ab}}}\int d^2{\mathbf{R}}^{\prime }\frac{{\rho }_{\mathrm{ind}}^{ns}({\mathbf{R}}^{\prime },\omega )}{\left|\mathbf{R}-{\mathbf{R}}^{\prime }\right|}],$$

$${\tilde{\mathbf{j}}}^{ns}(\vec{\theta },\omega )=$$

$$\displaystyle-\frac{\sigma^{(1)}_{s\omega}}{\epsilon^{\rm a}_{s\omega}D^{3}}\sqrt{f_{\vec{\theta}}}\,\nabla_{\vec{\theta}}\ \frac{({\vec{\theta}}-{\bf r}/D)\cdot{\bf d}_{s\omega}}{\bigl{\lvert}{\vec{\theta}}-{\bf r}/D\bigl{\rvert}^{3}}$$

$$+{\eta }_{s\omega }^{(1)}\int d^2{\vec{\theta }}^{\prime }\mathbf{M}(\vec{\theta },{\vec{\theta }}^{\prime })\cdot {\tilde{\mathbf{j}}}^{ns}({\vec{\theta }}^{\prime },\omega ),$$

$${\bf M}({\vec{\theta}},{\vec{\theta}}^{\prime})=\sqrt{f_{\vec{\theta}}f_{{\vec{\theta}}^{\prime}}}\nabla_{\vec{\theta}}\otimes\nabla_{\vec{\theta}}\,\frac{1}{\bigl{\lvert}{\vec{\theta}}-{\vec{\theta}}^{\prime}\bigl{\rvert}}.$$

$$\int d^2{\vec{\theta }}^{\prime }\mathbf{M}(\vec{\theta },{\vec{\theta }}^{\prime })\cdot {\vec{\epsilon }}_m({\vec{\theta }}^{\prime })=\frac{1}{{\eta }_m}{\vec{\epsilon }}_m(\vec{\theta })$$

$$\int d^2\vec{\theta }{\vec{\epsilon }}_m(\vec{\theta }){\vec{\epsilon }}_{m^{\prime }}(\vec{\theta })={\delta }_{m{m}^{\prime }}.$$

$${\tilde{\mathbf{j}}}^{ns}(\vec{\theta },\omega )=\sum _m{a}_m^{ns}(\omega ){\vec{\epsilon }}_m(\vec{\theta }).$$

$$a_m^{ns}(\omega )=\frac{1}{1-{\eta }_{s\omega }^{(1)}/{\eta }_m}\left[b_m^{ns}(\omega )+\frac{{\sigma }_{s\omega }^{(1)}}{{ϵ}_{s\omega }^{\mathrm{a}}{D}^3}{\mathbf{d}}_{s\omega }\cdot {\mathbf{e}}_m\right],$$

$$b_m^{ns}(\omega )=\int d^2\vec{\theta }\sqrt{f_{\vec{\theta }}}{\vec{\epsilon }}_m(\vec{\theta })\cdot {\mathbf{j}}_{\mathrm{NL}}^{ns}(\vec{\theta },\omega )$$

$${\bf e}_{m}({\bf r})=\int d^{2}{\vec{\theta}}\,\frac{{\vec{\theta}}-{\bf r}/D}{\bigl{\lvert}{\vec{\theta}}-{\bf r}/D\bigr{\rvert}^{3}}\rho_{m}({\vec{\theta}})$$

$${\mathbf{E}}_{\mathrm{ind}}^{ns}(\mathbf{r},\omega )=$$

$$\times \left[b_m^{ns}(\omega )+\frac{{\sigma }_{s\omega }^{(1)}}{{ϵ}_{s\omega }^{\mathrm{a}}{D}^3}{\mathbf{d}}_{s\omega }\cdot {\mathbf{e}}_m(\mathbf{r})\right].$$

$${\mathbf{E}}^{11}(\mathbf{R},\omega )=$$

$$-\frac{1}{{ϵ}_{\omega }^{\mathrm{ab}}}{\nabla }_{\mathbf{R}}\int d^2{\mathbf{R}}^{\prime }\frac{{\rho }_{\mathrm{ind}}^{11}({\mathbf{R}}^{\prime },\omega )}{\left|\mathbf{R}-{\mathbf{R}}^{\prime }\right|},$$

$$\vec{\epsilon }(\vec{\theta },\omega )={\vec{\epsilon }}_{\mathrm{ext}}(\vec{\theta },\omega )+{\eta }_{\omega }^{(1)}\int d^2{\vec{\theta }}^{\prime }\mathbf{M}(\vec{\theta },{\vec{\theta }}^{\prime })\cdot \vec{\epsilon }({\vec{\theta }}^{\prime },\omega ).$$

$$\vec{\epsilon }(\vec{\theta },\omega )=\sum _m\frac{c_m}{1-{\eta }_{\omega }^{(1)}/{\eta }_m}{\vec{\epsilon }}_m(\vec{\theta }),$$

$$c_m=\int d^2\vec{\theta }{\vec{\epsilon }}_m(\vec{\theta })\cdot {\vec{\epsilon }}_{\mathrm{ext}}(\vec{\theta },\omega ).$$

$${\mathbf{E}}^{11}(\vec{\theta },\omega )=\sum _m\frac{-{\vec{\xi }}_m\cdot {\mathbf{E}}_{\mathrm{ext}}}{1-{\eta }_{\omega }^{(1)}/{\eta }_m}\frac{\vec{\epsilon }(\vec{\theta })}{\sqrt{f_{\vec{\theta }}}},$$

$${\sigma }_{\omega }^{(1)}=$$

$$+\frac{e^2}{4\hslash }\left[\Theta (\hslash \omega -2E_{\mathrm{F}})+\frac{\text{i}}{\pi }\log \left|\frac{\hslash \omega -2E_{\mathrm{F}}}{\hslash \omega +2E_{\mathrm{F}}}\right|\right],$$

$${\mathbf{j}}_{\mathrm{NL}}^{33}={\sigma }_{3\omega }^{(3)}\left({\mathbf{E}}^{11}\cdot {\mathbf{E}}^{11}\right){\mathbf{E}}^{11},$$

$$b_m^{33}(\omega )=$$

Full text

1. I Introduction

2. II Theoretical model

3. III Results and discussion

4. IV Conclusions

5. A Basic concepts and approximations

6. B Electrostatic description of nonlinear plasmonic near-fields in 2D materials

7. C Nonlinear 2D plasmon–quantum emitter interaction

8. C.1 Two-level atom

9. C.2 Three-level V-type atom: Electromagnetically induced transparency

10. C.3 Three-level $\Lambda$-type atom: Temporal quantum control

I Introduction

Precise control over light-matter interactions on nanometric length scales presents opportunities in technologies such as light energy harvesting Atwater and Polman (2010), optical biosensing Fabrizio et al. (2016), and quantum information Lodahl et al. (2015). Plasmons, consisting of electromagnetic fields hybridized with collective charge-carrier oscillations, are ideal for this purpose, as they couple strongly with impinging light and can focus it into nanoscale volumes Maier (2007); Schuller et al. (2010). The interaction of enhanced local fields produced by resonantly illuminated plasmonic nanostructures with quantum emitters (QEs), such as atoms or molecules (natural or artificial), is of particular interest to the quantum optics and optical sensing communities Achermann (2010); Brolo (2012); Koenderink (2017). Typically, optimal coupling is achieved by engineering the geometry of a noble metal nanostructure so that a plasmon resonates with an electronic transition in a proximal dipole emitter (e.g., an exciton in a quantum dot). However, because noble metal plasmons cannot easily be modified using external stimuli Brown et al. (2015), resonant light-matter interactions in noble metal-QE nanocomposites typically lack the active tunability required by many photonic device functionalities, and are further hindered by large intrinsic Ohmic losses Khurgin (2015).

In its pristine form, graphene, the atomically thin carbon layer, is a zero-gap semiconductor that presents broadband 2.3% optical absorption from transitions between linear valence and conduction bands Castro Neto et al. (2009). When doped to a Fermi energy ${E_{\rm F}}$, Pauli blocking prohibits vertical interband absorption for photon energies $\hbar\omega\leq 2{E_{\rm F}}$, and graphene plasmons, consisting of coherently-coupled virtual intraband transitions, emerge within this gap. These resonances exhibit stronger confinement and longer lifetimes than their noble metal counterparts, and are readily tuned in an active manner via electrostatic gating Grigorenko et al. (2012). For these reasons, the interaction of plasmonic near fields from highly doped graphene with proximal QEs has been predicted to enable observable vacuum Rabi splittings Koppens et al. (2011), large Purcell enhancement factors Nikitin et al. (2011); Sloan et al. (2018), and electrical control of quantum states Manjavacas et al. (2012). Unfortunately, even at high doping levels, plasmon energies $\hbar\omega_{\rm p}$ in graphene nanostructures appear in the infrared or terahertz regimes, well-below the frequencies associated with long-lived electronic excitations in robust quantum light sources. Indeed, graphene structures with a characteristic size $D$ exhibit localized plasmons at energies scaling like $\hbar\omega_{\rm p}\propto\sqrt{{E_{\rm F}}/D}$, which for achievable values of ${E_{\rm F}}<1$ eV and $D>10$ nm have so far been demonstrated at mid-infrared and lower frequencies. Additionally, coupling to optical phonons persists as a strong source of damping for plasmons at energies above $\sim 0.2$ eV Yan et al. (2013), as do interband transitions unless $\hbar\omega_{\rm p}\lesssim{E_{\rm F}}$ García de Abajo (2014).

The conical electronic dispersion of graphene imposes anharmonic intraband charge-carrier motion, resulting in an intrinsically nonlinear response to external electromagnetic fields, including efficient generation of odd-ordered harmonics Mikhailov (2007); Ishikawa (2010); Cheng et al. (2015); Mikhailov (2016); Jiang et al. (2018); Soavil et al. (2018); Hafez et al. (2018). The synergistic combination of a large optical nonlinearity and strong plasmonic near-field enhancement, both arising from intraband transitions in doped graphene, is currently motivating intensive research efforts in the emerging field of nonlinear graphene plasmonics Cox and García de Abajo (2014); Jablan and Chang (2015); Christensen et al. (2015); Constant et al. (2016); Jadidi et al. (2016); Cox et al. (2017a, b); Kundys et al. (2018). In particular, 2D plasmons in graphene nanostructures have been predicted to generate harmonics with efficiency above that of similarly sized noble metal nanostructures Cox and García de Abajo (2014).

In this work, we propose to harness plasmon-enhanced harmonic generation in nanostructured graphene to couple the 2D material with a nearby QE, effectively bridging the energy mismatch between the electrically tunable infrared graphene plasmon and a near-infrared excitation in the emitter. We show that resonant excitation of a localized graphene plasmon can reduce the impinging light intensity required to generate a substantial third-harmonic (TH) near field, which can be actively tuned to the desired emitter resonance frequency by exploiting the unique electro-optic response of graphene. For the realistic parameters considered in this work, our simulations predict strong nonlinear plasmonic interactions that can drive fluorescence in two-level atoms and electromagnetically induced transparency (EIT) or coherent population control in three-level QEs. The present nonlinear coupling strategy can be used generally to probe nonlinear plasmonic near fields Metzger et al. (2017); Yang et al. (2017), while the electrical tunability of graphene plasmons can be exploited to control single-photon emission and actively manipulate quantum states.

II Theoretical model

A prototypical hybrid molecule consists of a doped graphene nanodisk with diameter $D=40$ nm and a two-level QE located a distance $d$ directly above its center, the latter characterized by the dipole moment $\vec{\mu}_{12}$ associated with the transition between states $\ket{1}$ and $\ket{2}$, which we consider to be oriented parallel to the graphene plane with magnitude $1\,e\times$ nm ($\sim 50$ Debye, commensurate with quantum dot excitons Stievater et al. (2001)). We consider local, isotropic dielectric media above and below the graphene layer characterized by the permittivities $\epsilon^{\rm a}_{\omega}$ and $\epsilon^{\rm b}_{\omega}$, respectively. As illustrated schematically in Fig. 1(a), an impinging light electric field ${\bf E}_{\rm ext}$ with frequency $\omega\approx\omega_{\rm p}$ (red beam) generates a plasmon-enhanced TH near-field, ${\bf E}^{33}_{\rm ind}({\bf r},\omega)\text{e}^{-\text{i}3\omega t}+{\rm c.c.}$, around the nanodisk (blue field), which we describe semianalytically by adopting a classical electrostatic eigenmode decomposition (see Ref. Cox et al. (2017b) and Appendix). Electron dynamics in the QE is governed by the density matrix master equation,

[TABLE]

where $\mathcal{H}$ is the system Hamiltonian and $\mathcal{L}[\rho]$ denotes the Lindblad superoperator accounting for decoherence. For a two-level system we write the Hamiltonian

[TABLE]

with $\hbar\varepsilon_{j}$ denoting the energy of the QE state $\ket{j}$. Defining the slowly-varying coherence elements $\tilde{\rho}_{12}=\rho_{12}\text{e}^{\text{i}3\omega t}$ and $\tilde{\rho}_{21}=\rho_{21}\text{e}^{-\text{i}3\omega t}$, we write the average QE transition dipole moment as ${\bf d}(t)=\vec{\mu}_{12}(\tilde{\rho}_{21}\text{e}^{\text{i}3\omega t}+\tilde{\rho}_{12}\text{e}^{-\text{i}3\omega t})$. Then, inserting Eq. (2) into Eq. (1) and using the expressions for $\mathcal{L}[\rho]$ and ${\bf E}^{33}_{\rm ind}$ provided in the Appendix, the density matrix equations of motion in the rotating-wave approximation are found to be

[TABLE]

where we adopt $\Gamma_{0}=1$ ns*-1* as a phenomenological relaxation rate from the excited state $\ket{2}$ to the ground state $\ket{1}$, $\Delta\equiv 3\omega-\varepsilon_{12}$ is the detuning of the TH near-field from the QE resonance energy $\hbar\varepsilon_{12}\equiv\hbar\varepsilon_{2}-\hbar\varepsilon_{1}$,

[TABLE]

is the Rabi frequency quantifying the direct coupling of the QE with the TH near-field produced by the 2D nanostructure, and

[TABLE]

characterizes the self-interaction strength of the induced QE dipole at frequency $s\omega$ enabled by graphene. In the above expressions,

[TABLE]

denotes the normalized electric field produced at ${\bf r}=(0,0,d)$ by integrating the so-called plasmon wave function (see Ref. Yu et al. (2017a) and Appendix) $\rho_{m}({\bf R})$ associated with mode $m$ over the 2D position vector ${\bf R}$ within the graphene nanostructure. The dimensionless parameters $\eta^{(1)}_{s\omega}\equiv 2\text{i}\sigma^{(1)}_{s\omega}/{s\omega}D(\epsilon^{\rm a}_{s\omega}+\epsilon^{\rm b}_{s\omega})$ in Eqs. (4) and (5), containing the dependence on the linear conductivities of extended graphene $\sigma^{(1)}_{s\omega}$ at frequency ${s\omega}$ (calculated in the local limit of the random-phase approximation García de Abajo (2014); Wunsch et al. (2006)), yield resonant spectral features in $\Omega$ and $G$ at frequencies ${s\omega}$ satisfying ${\rm Re}\{\eta^{(1)}_{s\omega}/\eta_{m}\}=1$, where $\eta_{m}$ are modal eigenvalues, with strengths determined by the corresponding dipolar coupling parameters $\vec{\xi}_{m}=\int d^{2}{\bf R}{\bf R}\rho_{m}({\bf R})/D^{3}$. TH generation is characterized by the local third-order conductivity of extended graphene $\sigma^{(3)}_{3\omega}$, for which we adopt the analytical result reported in Ref. Mikhailov (2016), obtained quantum-mechanically in the Dirac cone approximation, while $\zeta^{(3)}$ is a third-order coupling parameter (see Appendix). The conductivities employed here, accounting for both interband and intraband electronic transitions in extended graphene at zero temperature, combined with tabulated values of $\eta_{m}$, $\vec{\xi}_{m}$, and $\zeta^{(3)}$, can faithfully reproduce the linear and nonlinear polarizabilities of graphene nanostructures with sizes $\gtrsim 10$ nm predicted in fully-atomistic simulations Cox et al. (2017b).

III Results and discussion

Nonlinear plasmonic near-field coupling between the nanodisk and QE occurs when the transition energy $\hbar\varepsilon_{12}$ and the plasmon energy $\hbar\omega_{\rm p}$ satisfy $\varepsilon_{12}\approx 3\omega_{\rm p}$, as indicated in Fig. 1(b). The QE undergoes Rabi oscillations at a frequency $\Omega$, while the real and imaginary parts of $G$ contribute to a transition energy renormalization and decay rate enhancement, respectively Zhang et al. (2006); Malyshev and Malyshev (2011); Antón et al. (2012). Using plasmon wave functions and eigenvalues tabulated in Ref. Yu et al. (2017b) for a disk geometry, we plot $\Omega$ and $G$ in Figs. 1(c) and 1(d), respectively, for a linearly polarized field parallel to $\vec{\mu}_{12}$ impinging on the QE-nanodisk system in vacuum ($\epsilon^{\rm a}_{\omega}=\epsilon^{\rm b}_{\omega}=1$). With realistic graphene parameters (Fermi energy ${E_{\rm F}}\approx 0.3$–0.7 eV and $66$ fs plasmon lifetime, used throughout this work), a modest separation $d=20$ nm yields Rabi frequencies $\Omega\gtrsim\Gamma_{0}$, indicating strong driving of the QE via plasmon-assisted TH generation. Figures 1(c) and 1(d) show that the atom-plasmon interaction at the TH can be electrically tuned over a larger range of higher frequencies otherwise obtained though linear coupling at the fundamental frequency, while strong transition energy shifts (Re{G}) and decay rate enhancement (Im{G}) associated with $\omega_{\rm p}$ (e.g., in a linear coupling scheme) are avoided at the TH, where only weakly interacting, higher-order plasmon modes play a role.

For the system considered in Fig. 1 under continuous-wave (cw) illumination, we simulate the QE power absorption $P_{\rm abs}=\Gamma_{0}\hbar\varepsilon_{12}\rho_{22}$ from the graphene plasmon TH near-field by examining the steady-state density matrix obtained upon numerical integration of Eqs. (3). In Fig. 2(a) we plot the TH power absorption as a function of the impinging light frequency for different QE-graphene separation distances $d$ (upper panel) and light intensities $I_{\rm ext}$ (lower panel), maintaining ${E_{\rm F}}=0.5$ eV in graphene such that $\hbar\omega_{\rm p}=0.2694$ eV [see Fig. 1(c)] and $\varepsilon_{12}=3\omega_{\rm p}$. Under these conditions, the QE resonance does not interact directly with a plasmon mode of the graphene nanodisk, although a persisting self-interaction at the QE transition frequency, quantified by the parameter $G$, leads to enhanced decay rates and renormalized energies Zhang et al. (2006): this interaction is proportional to the square of the dipole coupling $\vec{\mu}_{12}\cdot{\bf e}_{m}$ [see Eq. (5)] and yields a blue-shifted and broadened power absorption peak at small separations (e.g., $d=10$ nm) but rapidly disappears with increasing separation distance, while the direct TH near-field interaction, which depends linearly on the dipole coupling, diminishes more slowly. In the present nonlinear plasmonic coupling scheme, the TH near-field enhancement of the QE excitation rate ($\Omega$) is independent of the self-interaction strength ($G$), unlike in the degenerate plasmon-exciton coupling scheme, where these competing phenomena are both maximized at the plasmon resonance frequency Doeleman et al. (2016); Koenderink (2017). Under strong driving fields, the TH near-field produces a clear signature of power broadening in the QE absorption spectrum Boyd (2008), leading towards optical bistability and hysteresis fueled by the QE self-interaction Malyshev and Malyshev (2011). In Fig. 2(b) we consider a situation where the illumination frequency is fixed at $\omega=\omega_{\rm p}=\varepsilon_{12}/3$ and the TH power absorption is actively controlled by tuning the Fermi energy ${E_{\rm F}}$ (doping charge density) in the graphene nanodisk.

The effective absorption cross-section of the hybrid system, defined as $P_{\rm abs}/I_{\rm ext}$, can reach values of approximately $\sim 0.06$ nm2 (for $I_{\rm ext}=1$ GW/m2 at $d=15$ nm), which is several orders of magnitude larger than the optimal cross-section of the dipole emitter in free space, $8\pi\mu_{12}^{2}\varepsilon_{12}/\hbar\Gamma_{0}c\approx 2.5\times 10^{-5}$ nm2, but significantly less than the geometrical area of the graphene nanodisk. We further remark that although these results are obtained for intense cw illumination, approaching the damage threshold for graphene, qualitatively similar behavior is expected for pulses of similar peak intensity and durations of several nanoseconds (i.e., with energy densities falling well-below the damage threshold Krauss et al. (2009); Currie et al. (2011); Roberts et al. (2011); Kiisk et al. (2013)).

The near-field produced by plasmon-assisted TH generation can be probed in an alternative manner by QEs with more complicated energy-level structures. In a three-level V-type atom, which we study in Fig. 3, EIT Fleischhauer et al. (2005) produces a dip in the absorption associated with one atomic transition when the other is coupled strongly to the TH near-field. This process is illustrated schematically in Fig. 3(a), where two external fields are applied to the hybrid system and the ground state $\ket{1}$ is coupled to two higher-energy states $\ket{2}$ and $\ket{3}$: we consider $\varepsilon_{12}=3\omega_{\rm p}$, while the energy associated with the transition $\ket{1}\leftrightarrow\ket{3}$ is arbitrarily far away from any prominent features in the nanodisk spectrum [see Fig. 3(b)]. By probing the $\ket{1}\leftrightarrow\ket{3}$ transition with a weak field ($\Omega_{\rm probe}=0.01\Gamma_{0}$) while resonantly coupling the $\ket{1}\leftrightarrow\ket{2}$ transition with the TH near-field (i.e., $\omega=\omega_{\rm p}=\varepsilon_{12}/3$), an absorption dip at resonance is observed when the QE and graphene nanodisk are in close proximity and/or the graphene plasmon is strongly driven [see Fig. 3(c)]. The active tunability of the graphene plasmon resonance can be exploited to actively switch EIT on and off, as we demonstrate in Fig. 3(d).

In Fig. 4 we seek to achieve electrical control of atomic state populations using the graphene nanodisk TH near-field. In this scenario, a QE with a lambda-type energy level configuration is located above the nanodisk and has orthogonal dipole moments associated with nearly-degenerate $\ket{1}\leftrightarrow\ket{2}$ and $\ket{2}\leftrightarrow\ket{3}$ transitions [we assume $\Delta=1$ meV; see Fig. 4(a)], which can be simultaneously driven by an optical pulse polarized $45^{\circ}$ from each of them Antón et al. (2012). We consider a Fourier-transform-limited Gaussian pulse of FWHM duration $\sigma=130$ fs, peak intensity $I_{\rm max}$, and carrier frequency $\omega$, such that $3\omega$ is $\pm\Delta$ away from $\varepsilon_{12}$ and $\varepsilon_{13}$, respectively. The effective maximum TH Rabi frequency $\mathinner{\!\left\lvert\Omega\right\rvert}$ for cw light at the peak pulse intensity is plotted in Fig. 4(b) for various Fermi energies, with the energies $\hbar\varepsilon_{12}/3$ and $\hbar\varepsilon_{13}/3$ indicated by vertical lines. Clearly, an asymmetry in the Rabi frequency driving the two nearly degenerate transitions can be introduced by modifying the Fermi energy, resulting in different levels of population transfer from the state $\ket{2}$ to state $\ket{3}$ [Fig. 4(c,d)].

IV Conclusions

In conclusion, we have demonstrated that plasmon-assisted up-conversion in graphene nanostructures can be used to resonantly excite a proximal QE, enabling energy transfer, EIT, or quantum state population control that would otherwise be prevented by the energy mismatch between infrared graphene plasmons and higher-energy electronic transitions in QEs. This scheme can also be used to probe the nonlinear optical response associated with graphene plasmons, potentially alleviating the sparsity of experimental work in the field of nonlinear graphene plasmonics. While the combination of an intrinsically-large nonlinear response and intense in-plane electric field enhancement associated with localized electrically tunable plasmons renders graphene an ideal material platform for nonlinear near-field coupling, this principle could be straightforwardly applied to other plasmonic materials. We anticipate that the strategy presented here can be used to couple electrically tunable graphene plasmons with higher-energy transitions in biologically interesting molecules and high-fidelity single photon sources, opening a wide range of applications in sensing and quantum nano-optics.

Acknowledgements.

The authors thank Renwen Yu for fruitful discussions and for providing the plasmon wave functions associated with a disk geometry. This work has been supported in part by the Spanish MINECO (MAT2017-88492-R and SEV2015- 0522), ERC (Advanced Grant 789104-eNANO), the European Commission (Graphene Flagship 696656), the Catalan CERCA Program, and Fundació Privada Cellex. \appendixpage\addappheadtotoc

We provide details on the semi-analytical plasmon wave function formalism used to describe the nonlinear plasmonic near-field produced by an arbitrary 2D material and its interaction with a proximal dipole emitter, summarize the theory of plasmon-assisted third-harmonic generation in a graphene nanostructure, and derive the Bloch equations governing electron dynamics in the two- and three-level atoms considered in the main text.

Appendix A Basic concepts and approximations

In the main text we consider a hybrid molecule consisting of a graphene nanostructure that occupies a finite region in the ${\bf R}=(x,y)$ plane and a proximal quantum emitter (QE) located at ${\bf r}=(x,y,z)$. The system is illuminated by an electromagnetic field ${\bf E}_{\rm ext}$ polarized in the graphene plane and oscillating at frequency $\omega\sim\omega_{\rm p}$, where $\omega_{\rm p}$ is the resonance frequency associated with a localized plasmon mode supported by the graphene structure. Through an $n^{\rm th}$-order nonlinear optical process, the graphene nanostructure generates a field ${\bf E}^{ns}_{\rm ind}({\bf r})$ oscillating at a harmonic frequency ${s\omega}$ ($\mathinner{\!\left\lvert s\right\rvert}\leq n$) that interacts with the QE. In what follows, we describe the nonlinear plasmonic near field generated by the graphene nanostructure and its interaction with the QE in a semi-analytical fashion by making the following simplifying assumptions:

• •

The 2D graphene nanostructure under consideration is much smaller than the light wavelength associated with $\omega_{\rm p}$, so that we can safely describe its response in the electrostatic limit.

• •

The QE is characterized by a point dipole placed at ${\bf r}=(x,y,z)$ with moment ${\bf d}_{s\omega}$ and resonance frequency $\sim{s\omega}$ associated with the transition between two discrete electronic states.

• •

The graphene optical response for each order $n$ and frequency ${s\omega}$ is characterized by an isotropic conductivity $\sigma^{(n)}_{s\omega}$ in the local limit, which is a reasonable approximation when the in-plane light momenta are much smaller than the involved electron momenta. This approximation is well justified for normally-impinging light interacting with localized plasmon modes in deeply-subwavelength 2D structures.

• •

The impinging light is sufficiently weak so that the optical response of graphene is well described by perturbation theory, maintaining zero electron temperature.

• •

The direct interaction of the external field with the QE is negligible compared to the near field produced by the nonlinear response of the graphene nanostructure.

• •

Inelastic scattering processes in the graphene nanostructure are frequency-independent, and occur at a phenomenological rate $\hbar\tau^{-1}=10$ meV ($\tau\approx 66$ fs) corresponding to a conservative electron mobility of $\sim 1300$ cm2 V*-1* s*-1* at ${E_{\rm F}}\approx 0.5$ eV.

Appendix B Electrostatic description of nonlinear plasmonic near-fields in 2D materials

In this section, we summarize the formalism presented in Ref. Cox et al. (2017b), which is used in this work to describe the nonlinear response of nanostructured graphene. To $n^{\rm th}$-order, an isotropic 2D material responds to an in-plane external field ${\bf E}_{\rm ext}\text{e}^{-\text{i}\omega t}+{\rm c.c.}$ by generating an induced charge $\rho^{ns}_{\rm ind}$ oscillating at a harmonic frequency ${s\omega}$ ($\mathinner{\!\left\lvert s\right\rvert}\leq n$) associated with the nonlinear 2D current ${\bf j}^{ns}$ through the continuity equation,

[TABLE]

We express ${\bf j}^{ns}$ self-consistently as the sum of currents produced by the direct nonlinear response to the total electric field at $\omega$ (${\bf j}^{ns}_{\rm NL}$) and the linear response to the induced electric field at the generated frequency ${s\omega}$ (${\bf E}^{ns}$), according to

[TABLE]

where $f_{\bf R}=1$ when ${\bf R}$ lies within the 2D material and is 0 everywhere else, effectively defining the structure morphology. For a QE with a resonance frequency $\sim{s\omega}$ positioned at a point ${\bf r}$ above the 2D material, the total electric field oscillating at this frequency within the 2D material becomes

[TABLE]

with the first and second terms arising from the QE transition dipole moment ${\bf d}_{s\omega}$ and the induced charge, respectively. Here, $\epsilon^{\rm ab}_{\omega}=(\epsilon^{\rm a}_{\omega}+\epsilon^{\rm b}_{\omega})/2$ denotes the average of the dielectric functions for media above ($\epsilon^{\rm a}_{\omega}$) and below ($\epsilon^{\rm b}_{\omega}$) the 2D layer.

Now, introducing the normalized 2D coordinate vector ${\vec{\theta}}\equiv{\bf R}/D$, where $D$ is the characteristic distance in the 2D geometry (e.g., the diameter of a disk or the side length of an equilateral triangle), we perform algebraic manipulations of Eqs. (A7-A9) to write the self-consistent current $\tilde{{\bf j}}^{ns}\equiv{\bf j}^{ns}/\sqrt{f_{\vec{\theta}}}$ as

[TABLE]

where $\eta^{(1)}_{s\omega}=\text{i}\sigma^{(1)}_{s\omega}/{s\omega}\epsilon^{\rm ab}_{s\omega}D$ is a dimensionless parameter and

[TABLE]

We identify ${\bf M}({\vec{\theta}},{\vec{\theta}}^{\prime})$ as a real and symmetric operator associated with the 2D geometry defined by $f_{\vec{\theta}}$ that admits real eigenmodes $\vec{\varepsilon}_{m}({\vec{\theta}})$ and eigenvalues $1/\eta_{m}$ satisfying

[TABLE]

and forming an orthonormal basis, that is,

[TABLE]

Following the procedure of Ref. Cox et al. (2017b), we expand the solution to Eq. (A10) in a sum of the eigenmodes of ${\bf M}({\vec{\theta}},{\vec{\theta}}^{\prime})$ according to

[TABLE]

Inserting Eq. (A14) into Eq. (A10), multiplying the resulting expression by $\vec{\varepsilon}_{m^{\prime}}({\vec{\theta}})$, and integrating over the 2D coordinates ${\vec{\theta}}$, we make use of Eqs. (A12) and (A13) to isolate the expansion coefficients

[TABLE]

where we have introduced the nonlinear coupling coefficient

[TABLE]

and the normalized electric field generated by the $m^{\rm th}$ eigenmode of the 2D geometry

[TABLE]

in terms of the so-called plasmon wave function (PWF) associated with it, $\rho_{m}({\vec{\theta}})\equiv\nabla_{\vec{\theta}}\cdot\sqrt{f_{\vec{\theta}}}\,\vec{\varepsilon}_{m}({\vec{\theta}})$ Yu et al. (2017a). In terms of normalized coordinates, Eq. (A7) provides the induced charge density $\rho^{ns}_{\rm ind}({\vec{\theta}},\omega)=-(\text{i}/{s\omega}D)\sum_{m}a^{ns}_{m}\rho_{m}({\vec{\theta}})$, from which we use Eq. (A15) to express the electric field acting on the QE at ${\bf r}$ as

[TABLE]

The first term in the above expression describes the near field produced by the 2D material oscillating at frequency ${s\omega}$ and generated directly by the external field at $\omega$, while the second term accounts for the field produced by the QE dipole oscillating at ${s\omega}$. In practice, the PWFs $\rho_{m}({\vec{\theta}})$, eigenmodes $\vec{\varepsilon}_{m}({\vec{\theta}})$, and eigenvalues $\eta_{m}$ of a given structure morphology are numerically computed once and for all using an electrostatic solver. In the main text, the PWFs and eigenvalues of a circular nanodisk are obtained from tabulated data reported in Ref. Yu et al. (2017b). From Eq. (A18) we conclude that the coupling of the QE dipole ${\bf d}_{s\omega}$ with the TH near-field produced by the $m^{\rm th}$ mode of a 2D geometry is proportional to the scalar product ${\bf d}_{s\omega}\cdot{\bf e}_{m}({\bf r})$, which is strongly dependent on the position of the dipole. In Fig. 5, we study this coupling for a dipole placed at a fixed distance away from the plane occupied by a circular disk and oriented either parallel or perpendicular to the polarization of the incident electric field. Evidently, maximal coupling should occur for a dipole aligned with the external field polarization and centered on the disk.

The 2D nanostructure response to a monochromatic external field ${\bf E}_{\rm ext}({\bf R},\omega)$ is dominated by the linear self-consistent field (i.e., $n=s=1$),

[TABLE]

where $\rho^{11}_{\rm ind}$ is the induced charge given by Eq. (A7). Following Ref. Yu et al. (2017a), we write the linear current as ${\bf j}^{11}=\sigma^{(1)}_{\omega}{\bf E}^{11}$, move to ${\vec{\theta}}$ coordinates, and define the normalized field $\vec{\varepsilon}\,({\vec{\theta}},\omega)\equiv D\sqrt{f_{\vec{\theta}}}\,{\bf E}^{11}({\vec{\theta}},\omega)$ to write

[TABLE]

Expanding the solution to the above expression in the eigenmodes of ${\bf M}({\vec{\theta}},{\vec{\theta}}^{\prime})$ [see Eqs. (A11-A13)], we obtain the normalized field

[TABLE]

where the expansion coefficients are

[TABLE]

For plane-wave illumination normal to the graphene plane, ${\bf E}_{\rm ext}({\bf R},\omega)$ is independent of ${\bf R}$ and Eq. (A21) yields the electric field

[TABLE]

where the dimensionless parameter $\vec{\xi}_{m}=\int d^{2}{\vec{\theta}}\,{\vec{\theta}}\,\rho_{m}({\vec{\theta}})$ characterizes the strength of the induced dipole Yu et al. (2017a). The in-plane near-fields given by Eq. (A23) depend on the linear conductivity $\sigma^{(1)}_{\omega}$ of the 2D material entering the parameter $\eta^{(1)}_{\omega}$. For graphene doped to a Fermi energy ${E_{\rm F}}$ and an inelastic electron scattering time $\tau$, we adopt the local conductivity obtained in the random-phase approximation at zero temperature,

[TABLE]

with the first and second terms accounting for intra- and inter-band optical transitions, respectively Wunsch et al. (2006).

The plasmon-enhanced in-plane field at frequency $\omega$, given to linear order by Eq. (A21), drives the nonlinear current ${\bf j}^{ns}_{\rm NL}$ entering Eq. (A16), producing the near-field described by Eq. (A18). For third-harmonic generation (i.e., $n=s=3$), the nonlinear source current is

[TABLE]

where $\sigma^{(3)}_{3\omega}$ is the third-order conductivity characterizing third-harmonic generation in the infinitely-extended 2D material. Here, we adopt the analytical result reported in Ref. Mikhailov (2016), which is obtained from a local quantum-mechanical description of graphene in the Dirac cone approximation at zero temperature. Inserting Eq. (A23) into Eq. (A25), we thus obtain from Eq. (A16)

[TABLE]

where

[TABLE]

is a unitless parameter. Mode dipoles $\vec{\xi}_{m}$ and eigenvalues $\eta_{m}$ of a given structure morphology are numerically computed once and for all using an electrostatic solver. Throughout the main text we consider a circular graphene nanodisk, for which the linear $\vec{\xi}_{m}$ and nonlinear $\zeta^{(3)}_{mm^{\prime}m^{\prime\prime}m^{\prime\prime\prime}}$ dipolar coupling parameters for modes up to order $m\leq 3$ have been tabulated in Ref. Cox et al. (2017b).

Appendix C Nonlinear 2D plasmon–quantum emitter interaction

Electron dynamics in the two- and three-level quantum emitters (QEs) considered in the main text is governed by the Lindblad master equation for the time-evolution of the density matrix,

[TABLE]

where $\mathcal{H}$ is the system Hamiltonian and $\mathcal{L}[\rho]$ denotes the Lindblad superoperator accounting for decoherence. In what follows, we consider that the QE states $\ket{i}$ with energies $\hbar\varepsilon_{i}$ form a complete and orthonormal set. This allows us to express the density matrix as $\rho=\sum_{ij}\rho_{ij}\sigma_{ij}$, where $\sigma_{ij}\equiv\ket{i}\bra{j}$, noting that the condition $\rho_{ij}=\rho^{*}_{ji}$ is guaranteed by the Hermitian property $\rho^{\dagger}=\rho$. Energies and dipole moments associated with transitions $\ket{i}\leftrightarrow\ket{j}$ are denoted by $\varepsilon_{ij}\equiv\varepsilon_{j}-\varepsilon_{i}$ and $\vec{\mu}_{ij}$, respectively. For simplicity we assume that all incoherent decay processes occur at the same phenomenological rate $\Gamma_{0}$.

C.1 Two-level atom

For a two-level QE with a ground state $\ket{1}$ and excited state $\ket{2}$ interacting with the classical monochromatic field ${\bf E}({\bf r},t)={\bf E}^{ns}_{\rm ind}({\bf r},\omega)\text{e}^{-\text{i}{s\omega}t}+{\rm c.c.}$, where ${s\omega}\sim\varepsilon_{12}$, we insert the Hamiltonian

[TABLE]

and the Lindblad operator

[TABLE]

into Eq. (A28). Making use of Eq. (A18) and writing the QE dipole moment as ${\bf d}(t)=\vec{\mu}_{12}(\tilde{\rho}_{21}\text{e}^{-\text{i}{s\omega}t}+\tilde{\rho}_{12}\text{e}^{\text{i}{s\omega}t})$, where $\tilde{\rho}_{21}=\tilde{\rho}^{*}_{12}=\rho_{21}\text{e}^{\text{i}{s\omega}t}$, the density matrix equations of motion in the rotating-wave approximation are obtained as

[TABLE]

where $\Delta\equiv{s\omega}-\varepsilon_{12}$ is the detuning parameter,

[TABLE]

is the Rabi frequency, which describes the direct coupling between the QE and the nonlinear near-field produced in the 2D material at frequency ${s\omega}$, and

[TABLE]

quantifies the self-interaction of the induced QE dipole at $s\omega$ enabled by the 2D nanostructure.

C.2 Three-level V-type atom: Electromagnetically induced transparency

For a V-type atom, we consider a QE with ground state $\ket{1}$ and excited states $\ket{2}$ and $\ket{3}$ interacting with the classical light field

[TABLE]

where ${s\omega}\sim\varepsilon_{12}$ and $\omega_{\rm probe}\sim\varepsilon_{13}$, but $\mathinner{\!\left\lvert\varepsilon_{2}-\varepsilon_{3}\right\rvert}\gg\Gamma_{0}$ (i.e., levels $\ket{2}$ and $\ket{3}$ are spaced far away from one another in energy so that only ${\bf E}^{ns}_{\rm ind}$ couples with $\vec{\mu}_{12}$ and only ${\bf E}_{\rm probe}$ couples with $\vec{\mu}_{13}$). In the rotating-wave approximation, the Hamiltonian describing this system is

[TABLE]

and the Lindblad operator is

[TABLE]

where $\Gamma_{j}$ is the decay rate of state $\ket{j}$ to $\ket{1}$. Inserting the above expressions into Eq. (A28), we obtain the equations of motion for the density matrix elements:

[TABLE]

where $\Omega_{2}=\Omega$ [see Eq. (A32)], $\Omega_{3}=\vec{\mu}_{13}\cdot{\bf E}_{\rm probe}/\hbar$, $\Delta_{2}\equiv{s\omega}-\varepsilon_{12}$, $\Delta_{3}\equiv\omega_{\rm probe}-\varepsilon_{13}$, and we have transformed the coherence elements according to

[TABLE]

In practical calculations, we consider a resonant driving nonlinear near-field (i.e., $\Delta_{2}=0$) with a weakly-driving probe field, such that $\Omega_{3}=0.01\Gamma_{0}$, while setting $\Gamma_{2}=\Gamma_{3}=\Gamma_{0}$ for simplicity. Electromagnetically-induced transparency is observed in the probe field absorption, which is proportional to ${\rm Im}\{\tilde{\rho}_{31}\}$, as the detuning $\Delta_{3}$ is varied.

C.3 Three-level $\Lambda$-type atom: Temporal quantum control

A $\Lambda$-type atom is characterized by two low-energy, nearly-degenerate states $\ket{1}$ and $\ket{2}$ that are each optically coupled to a higher-energy excited state $\ket{3}$. We study the interaction of the atom with a classical pulsed laser field

[TABLE]

where $t_{\rm p}$ is the temporal center of the impinging pulse and $\tau_{\rm p}$ quantifies its duration, which is considered to be much shorter than the QE excited state lifetime $\Gamma_{0}^{-1}$. We define the energy half-splitting between states $\ket{1}$ and $\ket{2}$ as $\Delta\equiv(\varepsilon_{2}-\varepsilon_{1})/2$ and assume a scenario in which the pulse carrier frequency is centered between the transitions $\ket{1}\leftrightarrow\ket{3}$ and $\ket{2}\leftrightarrow\ket{3}$, such that ${s\omega}=\varepsilon_{13}-\Delta=\varepsilon_{23}+\Delta$. The system Hamiltonian is then written in the rotating-wave approximation as

[TABLE]

while we write the Lindblad operator as

[TABLE]

where $\Gamma_{j}$ denotes the decay rate of the excited state $\ket{3}$ to state $\ket{j}$. Assuming that the pulse bandwidth is much larger than the energy splitting $2\Delta$, we use the above expressions in Eq. (A28) to write the density matrix equations of motion as

[TABLE]

In the above expressions, we have defined

[TABLE]

and

[TABLE]

where $s\omega_{j}=\varepsilon_{j3}$, and we have transformed the coherence elements according to $\tilde{\rho}_{j3}=\rho_{j3}\text{e}^{-\text{i}s\omega t}$ ($j\neq 3$). For the calculations presented in the main text, we take the transition dipole moments $\vec{\mu}_{13}$ and $\vec{\mu}_{23}$ to be orthogonal and each oriented $45^{\circ}$ degrees from the impinging field polarization, so that in the absence of material dispersion (e.g., without plasmonic enhancement) each of the transitions would be driven equally. For simplicity we set all decay rates equal to $\Gamma_{0}$.