Nonlinear optical response of a local surface plasmon coupled to a 2D material

TL;DR

This paper theoretically investigates a coupled system of a metal nanoparticle's surface plasmon and a monolayer transition-metal dichalcogenide, revealing strong nonlinear and quantum optical effects with potential applications in metamaterials and quantum photonics.

Contribution

It introduces a novel model of a nonlinear oscillator combining surface plasmons and 2D materials, predicting enhanced nonlinearity and quantum features not previously demonstrated.

Findings

Seven orders of magnitude increase in nonlinearity compared to silicon photonic cavities.

Observation of quantum features such as antibunching and non-Gaussianity.

Potential for creating nanoscale nonlinear optical metamaterials.

Abstract

We present a theoretical study of the optical response of a nonlinear oscillator formed by coupling a metal nanoparticle local surface plasmon resonance to excitonic degrees of freedom in a monolayer transition-metal dichalcogenide. We show that the combined system should exhibit strong anharmonicity in its low-lying states, predicting for example a seven order-of-magnitude increase in nonlinearity relative to a silicon photonic crystal cavity. Then, we demonstrate that such system exhibits strong quantum features such as antibunching and non-Gaussianity. Arrays of such nanoscale nonlinear oscillators could be used to realize novel optical metamaterials; alternatively, an individual nanoparticle-monolayer construct could be coupled to an optical resonator to mediate efficient input-output coupling to propagating fields.

Equations

$$H_a$$

$$H_{\mathrm{drive}}$$

$$\chi +i{\chi }^{\prime }=\frac{27{ϵ}_0\hslash {\omega }_a^2}{4}\int d^{3}r{\chi }^{(3)}(\mathbf{r})|{\mathbf{f}}_{\parallel }(\mathbf{r}){|}^4,$$

$$\dot{Q}=-\frac{i}{\hslash }[Q,H]+\sum _j\left(L_j^{†}Q{L}_j-\frac{1}{2}{L}_j^{†}{L}_{j}Q-\frac{1}{2}Q{L}_j^{†}{L}_j\right).$$

$$\dot{a}=-i{\Delta }_{a}a-i\chi (a^{†}a)a-\frac{\gamma }{2}a-\frac{{\chi }^{\prime }}{2}(a^{†}a)a-\sqrt{\kappa }\epsilon ,$$

$${\rho }_{22,\mathrm{ss}}=\frac{4\kappa {\epsilon }^2/(\gamma +{\chi }^{\prime }{)}^2}{1+(4{\chi }^2+4\kappa {\epsilon }^2)/(\gamma +{\chi }^{\prime }{)}^2}{\rho }_{11,\mathrm{ss}}.$$

$$\delta [\rho ]=S(\tau )-S(\rho )=H(\sqrt{\det \sigma })-S(\rho ),$$

$$H$$

$$H_b$$

$$H_c$$

$$H_I$$

$$H_{\mathrm{drive}}$$

$$\dot{b}$$

$$\dot{c}$$

$$\nabla \cdot \mathbf{D}=0,\nabla \cdot \mathbf{B}=0,$$

$$\nabla \times \mathbf{E}=-\frac{\partial \mathbf{B}}{\partial t},\nabla \times \mathbf{H}=\frac{\partial \mathbf{D}}{\partial t},$$

$$P_j={ϵ}_0\sum _{k=1}^3{\chi }_{jk}^{(1)}{E}_k,$$

$$\mathbf{P}={ϵ}_0\sum _{n>0}{\mathbf{\chi }}^{(n)}:{\mathbf{E}}^{\otimes n},$$

$$P_j={ϵ}_0\left[\sum _k{\chi }_{jk}^{(1)}{E}_k+\sum _{k,l}{\chi }_{jkl}^{(2)}{E}_k{E}_l+\sum _{k,l,m}{\chi }_{jklm}^{(3)}{E}_k{E}_l{E}_m+\cdots \right].$$

$$\mathbf{D}=\sum _{n>0}{\mathbf{ϵ}}^{(n)}:{\mathbf{E}}^{\otimes n},$$

$$\mathbf{E}=\sum _{n>0}{\mathbf{\eta }}^{(n)}:{\mathbf{D}}^{\otimes n},$$

$$\mathbf{E}=\sum _{n>0}{\mathbf{\eta }}^{(n)}:{\left[\sum _{m>0}{\mathbf{ϵ}}^{(m)}:{\mathbf{E}}^{\otimes m}\right]}^{\otimes n}.$$

$$D_{+}$$

$$D_{-}$$

$${ϵ}_{++}^{(1)}={ϵ}_{--}^{(1)}={ϵ}_0(1+{\chi }^{(1)}),{ϵ}_z^{(1)}={ϵ}_0,{ϵ}_{++++}^{(3)}={ϵ}_{----}^{(3)}=3{ϵ}_0{\chi }^{(3)},$$

$$H$$

$$\mathbf{D}=iϵ(\mathbf{r})\sqrt{\frac{\hslash \omega }{2}}(a\mathbf{f}(\mathbf{r})-a^{†}{\mathbf{f}}^{\ast }(\mathbf{r})),$$

$$\int d^{3}rϵ(\mathbf{r}){\mathbf{f}}^{\ast }(\mathbf{r})\cdot \mathbf{f}(\mathbf{r})=1.$$

$$\mathbf{E}=\frac{1}{ϵ(\mathbf{r})}\mathbf{D}=i\sqrt{\frac{\hslash \omega }{2ϵ(\mathbf{r})V}}(a{e}^{i\mathbf{k}\cdot \mathbf{z}}-a^{†}{e}^{-i\mathbf{k}\cdot \mathbf{z}}),$$

$$\mathbf{\Lambda }=\sqrt{\frac{\hslash }{2\mu \omega }}(a\mathbf{u}(\mathbf{r})-a^{†}{\mathbf{u}}^{\ast }(\mathbf{r})),$$

$$H=H_0+H_{\mathrm{nl}},$$

$$H_0$$

$$H_{\mathrm{nl}}$$

$$\mathbf{D}(\mathbf{r},t)=iϵ(\mathbf{r})\sqrt{\frac{\hslash {\omega }_a}{2}}\mathbf{f}(\mathbf{r})(a(t)-a^{†}(t)),$$

$$H_{\mathrm{nl}}=\int d^{3}r\frac{9}{4}{ϵ}_0\mathrm{Re}[{\chi }^{(3)}(\mathbf{r})]|{\mathbf{E}}_{\parallel }{|}^4,$$

Full text

• April 2019

1 Introduction

The quest for realizing a nanophotonic nonlinear oscillator has remained as a formidable challenge up until now, despite that many discussions on how Kerr optical nonlinearity can be used to realize such a nonlinear oscillator [1, 2]. The hardship is closely linked to another equally difficult quest for nanophotonic single-photon optical nonlinearity. The optical nonlinear oscillator, if realized, will be a game changer as it would operate at GHz or even at THz in a scalable platform. It is notable that many materials with a large Kerr coefficient have been reported so far, and yet, none of these succeeded to realize such a nonlinear oscillator. As we will explain in this paper, merely having a large Kerr nonlinearity is not sufficient since those materials with large Kerr nonlinearity tend to have also strong two-photon absorption, which forestalls the nonlinear response.

Atomically thin 2D materials are known to possess enormous optical nonlinearities within subatomic thickness [3, 4, 5, 6]. However, such a small thickness of the 2D materials is usually viewed as a drawback as their incorporation in conventional optical resonators provides disappointing nonlinear-optical performance because of small filling fraction in wavelength-scale volumes. For example, the excellent review papers by Sun et al. [3] and Autere et al. [4] illustrate the integrated optical modulation capability of 2D materials. However, the authors of these papers clearly pointed out that, to fully harness the potential of the strong nonlinearity of 2D materials at the atomic scale, it is necessary to circumvent the issues of inefficient light-matter interaction due to 2D material’s subnanometer thickness.

In order to accomplish a large optical nonlinearity in a small volume, various approaches have been extensively studied to date. For example, Nielsen *et al. *studied the plasmonic nanofocus driving strong four-wave mixing [7]. However, the authors acknowledged that their accomplished optical nonlinearity is far below the required level for quantum effects.

We propose a novel approach to harnessing the extreme field confinement of localized surface plasmon resonances (LSPR) for nonlinear optics, using 2D monolayer materials. This approach is a unique combination of the freedom to adjust the real and the imaginary values of the third-order Kerr nonlinearity from the excitons of 2D materials and the field enhancement from plasmonic coupling. Our result clearly demonstrates an evident quantum phenomenon – the strong anti-bunching effect. Thus, our result is a crucial step towards the long-standing quest for the strong quantum effect of nano-structured materials.

A LSPR of a metal nanoparticle (MNP) may be modeled as a quantized harmonic oscillator [8]. For a given resonance frequency, the energy states form a harmonic ladder broadened by Ohmic loss of the metal [9]. The LSPR field is confined within a volume slightly larger than the MNP, which is significantly smaller than the free-space wavelength of light at the LSPR resonant frequency [10]. We analyze a scheme for turning an LSPR into a quantum nonlinear resonator via near-field coupling to an atomically thin 2D material, such as a transition metal dichalcogenide (TMD). Two key considerations motivate this proposal: that monolayer TMDs can provide strong Kerr nonlinearity with a favorable coherence-to-decoherence ratio [6], and that monolayer TMDs are sufficiently thin to overlap substantially with the nanoscale LSPR evanescent field. Then, we proceed to demonstrate the remarkable quantum features of such systems, namely, photon antibunching and non-Gaussian state generation.

The structure of this article is as follows: section 2 builds up the concept of nonlinear oscillators, section 3 describes the prescription for highly nonlinear oscillators, section 4 models the quantum dynamics of the nonlinear oscillator and presents the simulation results, and finally, a conclusion and discussions follow.

2 Nonlinear oscillator

There are many complementary schemes for inducing nonlinearity in an optical oscillator [2, 1, 11]. In the (arguably) most widely applicable approach, a Kerr-nonlinear medium can be introduced into an otherwise passive optical cavity. If the ratio of nonlinear Kerr effect to absorption loss is sufficiently large, such a system may be useful for quantum optics and quantum information processing. Unfortunately, it is difficult to find materials with a large nonlinearity-to-loss ratio, which are amenable to incorporation within conventional or nanofabricated optical resonators. For example, although some bulk materials such as silicon have a significant Kerr optical nonlinearity at near-infrared optical frequencies, they also have a large imaginary nonlinear response (two-photon absorption) that gives rise to comparably large losses [12]. Setting aside the issue of losses, nanophotonic resonators incorporating conventional optical materials such as silicon generally do not possess a sufficiently large Kerr nonlinearity at low photon numbers for quantum information processing applications.

Monolayer MoS2 has an order-of-magnitude larger real $\chi^{(3)}$ than bulk silicon and can have a significant real-to-imaginary ratio of $\chi^{(3)}$ [6]. Its 2D structure furthermore makes it a promising candidate to induce a large coherent nonlinearity through direct coupling to the strongly confined field of an LSPR. Such an LSPR-TMD system could be efficiently coupled to optical input-output modes by incorporation within an optical cavity, in a manner analogous to the way that a Josephson-Junction nonlinear LC oscillator can be coupled to propagating fields by a microwave resonator. With regard to the potential scalability of such a device concept, we note that, while automated assembly of the type of system we envision could be challenging to realize, the material components we require may be producible in bulk by chemical MNP synthesis [13] and chemical-vapor-deposition (CVD) growth of TMDs [15, 14]. One can alternatively envision a large number of MNPs dispersed on a TMD membrane to form an array of nonlinear oscillators with complex couplings in both real and frequency spaces; such an optical metamaterial could provide rich nonlinear dynamics for reservoir computing-type architectures [16].

Let us consider an MNP driven by a classical light source. A monolayer TMD is positioned nearby the MNP so that the LSPR field has a considerably large overlap with the monolayer TMD (see Fig. 1). To understand the system, we adopt the quantum master equation [18, 17], which is particularly suitable for optical frequencies. In this framework, the unitary evolution of the system is modeled through a Hamiltonian $H=H_{a}+H_{\mathrm{drive}}$ where [21, 8, 19, 2, 20, 11, 22] (see A):

[TABLE]

where $\omega_{a}$ is the resonant frequency of the LSPR field, $\chi$ is the Kerr nonlinear coefficient, and $a,a^{\dagger}$ are the annihilation and the creation operators of the LSPR field. Here, $\kappa$ is the coupling coefficient, $\varepsilon$ is the drive, and $\omega$ is the frequency of the drive field. The nonunitary evolution of the system undergoes dissipation through the Lindblad operators $L_{1}=\sqrt{\gamma}a$ with $\gamma=\gamma_{r}+\gamma_{o}$, which includes both the radiative ($\gamma_{r}$) and nonradiative ($\gamma_{o}$) decays where the nonradiative Ohmic loss is dominant ($\gamma\simeq\gamma_{o}\gg\gamma_{r}$) [9], and $L_{2}=\sqrt{(\chi^{\prime}/2)}aa$, which is the nonlinear loss (two-photon absorption) through the coupling with the nonlinear 2D material. As a consequence, the effective energy levels appear as in Fig. 1. Without nonlinear coupling to the 2D material, the energy levels for a given frequency $\omega_{a}$ are a linear harmonic ladder with an equal level spacing of $\hbar\omega_{a}$, with a level broadening $\hbar\gamma$. The nonlinearity modifies the energy levels in two ways: both the energy level shift and the level broadening vary with LSPR photon number, resulting in an anharmonic energy ladder. Assuming a single-frequency optical driving field, it is apparent that, in order not to excite an effectively linear (trivial) response, the nonlinear energy shifts must be larger than the level broadenings.

The nonlinear energy shift and broadening depend on the coefficients $\chi$ and $\chi^{\prime}$ (both real-valued), which are calculated using the results in [21, 19, 23, 24] (see the detailed derivation in B):

[TABLE]

where $\epsilon_{0}$ is the vacuum permittivity and $\chi^{(3)}(\bm{r})$ is the position ($\bm{r}$) dependent third-order nonlinear susceptibility. Since monolayer TMDs have a negligibly small optical response to out-of-plane-polarized electric fields [25, 26], we count only the in-plane LSPR field component (with respect to the monolayer 2D material): $\bm{f}_{\parallel}(\bm{r})$ is the in-plane component of the LSPR mode function $\bm{f}(\bm{r})$, normalized as $\int d^{3}r~{}\epsilon_{0}\epsilon_{r}(\bm{r})|\bm{f}(\bm{r})|^{2}=1$ with the position dependent relative permittivity $\epsilon_{r}(\bm{r})$, Evidently, the overlap between the LSPR mode field and the 2D material volume must be large to increase $\chi$. We note that a large overlap also increases the nonlinear broadening of two-photon absorption. Hence, the the ratio $\chi/\chi^{\prime}$ must be large, which depends in general on the material and the operating frequency.

To appreciate how difficult it is to make a practical nonlinear (anharmonic) oscillator, let us consider an optical resonator filled with silicon. We make an example of the state-of-the-art smallest photonic crystal cavity to increase the nonlinear energy shift: a cavity of volume $V_{c}=0.02$ $\mu$m3 having a cavity $Q=2\times 10^{5}$ [27]. Then, using $\mathrm{Re}[\chi^{(3)}]\sim 4.8\times 10^{-20}$ m2/V2 of silicon at 1 eV photon frequency (energy) [12], Equation (3) leads to $\hbar\chi\sim 1\times 10^{-9}$ eV. However, the linear cavity-lifetime-induced level broadening (ignoring the intrinsic linear loss of silicon) is $\hbar\gamma\sim 5\times 10^{-6}$ eV. Hence, the line broadening is dominant over the nonlinear energy level shift at low photon number. Although quantum well materials such as GaAs/InGaAs have appreciable Kerr nonlinearity ($\mathrm{Re}[\chi^{(3)}]=1.3\times 10^{-19}$ m2/V2 at 1.3 $\mu$m wavelength [28]), the two-photon absorption is large ($\mathrm{Im}[\chi^{(3)}]=1.4\times 10^{-18}$ m2/V2 [28]) and the linear absorption is excessively large ($\alpha>10$  cm*-1*) [30, 29] so that the nonlinearity-to-broadening ratio is even less favourable than our silicon example. Chalcogenide glass materials are not much better due to the same reason of large linear and two-photon absorption [31].

3 Localized surface plasmon coupled to 2D materials

We will show below that our proposed system could achieve few-photon nonlinear energy shifts that dominate broadening, and thus realize a quantum nonlinear oscillator. We first discuss several design aspects for the system. We utilize the large ratio of the real and the imaginary values of $\chi^{(3)}$ in a monolayer MoS2 by detuning slightly from the two-photon resonant frequency 1.06 eV [6]. We must also design an MNP to support an LSPR resonance near 1.06 eV with a sufficiently small $\hbar\gamma$. The resonance frequency $\omega_{a}$ of the LSPR mode is well known to depend upon the geometry of the MNP [32, 33, 34], with thin MNPs exhibiting lower resonance frequencies [35]. The reason is that the portion of the LSPR field inside the metal reduces as the metal layer becomes thinner, which tends to red-shift the resonance frequencies [34].

We need to maximize the in-plane component of the LSPR field that overlaps the monolayer MoS2. A natural choice of geometry is a disk shape MNP with a large ratio of diameter to thickness. A large diameter improves the in-plane LSPR field component while a thinner disk shape pulls $\omega_{a}$ down from the nominally visible plasmonic resonances of silver. Using finite-element-method (FEM) software (COMSOL), we found that a diameter-to-thickness ratio of 22 for a silver disk should provide an LSPR resonance frequency near the target. On the other hand, the value of $\gamma$ of an LSPR is known to be independent of shape and size, since it is completely determined by the complex dielectric function of the material; once the material and operating frequency are determined, there is not much one can do to adjust $\gamma$ [34]. For this theoretical study, we take from the work of Wang and Shen a Q-factor $Q=\omega_{a}/\gamma=50$ [34]. In reality, the surface imperfections may reduce the Q-factor.

We also need to consider potential impacts on the intrinsic optical properties of a TMD monolayer when an MNP is placed nearby—an excessively small gap between the two may alter the band structure of the 2D material and impact the nonlinearity adversely. To address this problem, we consider a tight binding band Hamiltonian $H_{\mathrm{band}}=\sum_{<i,j>}J_{ij}c_{i}^{\dagger}c_{j}+\sum_{i}L_{i}c_{i}^{\dagger}c_{i}+\sum_{i}M_{i}(c_{i}^{\dagger}d+c_{i}d^{\dagger})$ where $c_{i}^{\dagger},d^{\dagger}$ are the second-quantized fermionic creation operators for the electrons in the 2D material lattice and the impurity nearby, respectively [36]. The first term is the nearest-neighbor (nn) hopping energy, the second is the on-site energy, and the third is the interaction energy between the 2D material and the impurity sites. We note that $J_{ij}$ and $M_{i}$ originate from the Coulomb potential. It is customary to ignore the next-nearest neighborhood (nnn) hopping since the matrix element for this remote potential is negligibly small compared to the nearest-neighbor hopping. Hence, to preserve the band structure intact, we require the condition $|M_{i}|\ll|J_{ij}|$, which can be fulfilled if the distance between the 2D lattice sites and the impurity sites becomes larger than the next-nearest neighbor (nnn) distance of the 2D lattice. When the hexagonal atomic arrangement of monolayer TMDs is considered, the nnn distance is 0.32 nm (MoS2), approximately 1.7 times longer than nn. Indeed, Liu *et al. *recently discussed the encapsulation of monolayer TMDs and found that using a few-layer hexagonal boron-nitride (hBN) buffer preserved the properties of monolayers successfully [37]. Man *et al. *discussed the successful protection of monolayer properties using even a single layer hBN buffer [38]. Hence, in this theoretical paper, we set the minimum distance between the monolayer TMD and the MNP as 0.3 nm, corresponding to the thickness of a single layer hBN [39]. Additionally, since the operating frequency is detuned from both the one- and two-photon excitonic resonance of the monolayer TMD, exciton generation will be minimal and we neglect metal-dielectric effects such as mirror charges.

4 Dynamical model and simulation results

To understand the dynamics of the system, we consider the adjoint master equation for a Heisenberg picture operator $Q$ as [17]:

[TABLE]

Then, we easily obtain the dynamical equation for the LSPR field using the Hamiltonian in Equation (1) and (2) and the Lindblad operators:

[TABLE]

where $\Delta_{a}=\omega_{a}-\omega$. In the limit of considering only the lowest three energy levels with a moderate drive power, the steady-state population in the second excited level can be expressed using the population in the first level as (C):

[TABLE]

Note that for a linear system ($\chi=\chi^{\prime}=0$), the population is $\rho_{22,\mathrm{ss}}=\rho_{11,\mathrm{ss}}(4\kappa\varepsilon^{2}/\gamma^{2})/(1+4\kappa\varepsilon^{2}/\gamma^{2})$, which exhibits the conventional saturation feature with respect to the drive power. The nonlinearity (anharmonicity) kicks in if the quantity $4\chi^{2}/(\gamma+\chi^{\prime})^{2}$ is sufficiently large to suppress $\rho_{22,\mathrm{ss}}$. Hence, we should adjust the system design parameters to maximize the nonlinear figure of merit (FOM) $4\chi^{2}/(\gamma+\chi^{\prime})^{2}$.

The FOM depends mainly on two design parameters: the MNP size and the drive frequency detuning. The MNP size determines the overlap between the monolayer TMD and the LSPR field, scaling $\chi$ and $\chi^{\prime}$ in the same fashion. A large MNP has a more delocalized field outside the MNP and, thus, the relative portion overlapping the thin nonlinear 2D material is small. On the other hand, detuning of the drive frequency adjusts the ratio between $\chi$ and $\chi^{\prime}$. The values $\mathrm{Im}[\chi^{(3)}]$ of a monolayer MoS2 roughly follows a typical Lorentzian lineshape while $\mathrm{Re}[\chi^{(3)}]$ is the frequency-derivative of $\mathrm{Im}[\chi^{(3)}]$ [6]. Hence, FOM is expected to be nearly zero at the two-photon resonance of the monolayer TMD exciton while the maximum occurs an appropriately detuned frequency. D presents the detailed design considerations.

To maximize the overlap, we place two monolayer MoS2 planes sandwiching the disk-shaped MNP, with hBN buffers separating the MNP from each of the monolayers. Using an FEM software package (COMSOL), we calculate the detailed field distribution of a given MNP design and subsequently determine $\chi,\chi^{\prime}$ according to Equation (3). Fig. 2 shows the calculated FOM as a function of the silver disk MNP’s diameter (while fixing the diameter-to-thickness ratio at 22) and the drive frequency detuning from the monolayer MoS2’s two-photon resonance at 1.06 eV. We note that with the fixed diameter-to-thickness ratio of 22, the MNP diameter must be larger than 4 nm so that the MNP disk thickness is larger than the atomic size of silver. Remarkably the nonlinear FOM is larger than unity over a wide range of parameters; the unity FOM value is the on-set for a system to be nonlinear. The maximum FOM of 10.7 occurs at a detuning of $-11$ meV with the smallest MNP diameter of 4 nm. The real and the imaginary values of $\chi^{(3)}$ of the monolayer MoS2 at this detuning are $4.0\times 10^{-19}$ m2/V2 and $1.0\times 10^{-19}$ m2/V2, respectively [6]. The subsequent $\hbar\chi$ and $\hbar\chi^{\prime}$ values are obtained as 53 meV and 13 meV, respectively. It is noteworthy that the improvement of FOM by placing additional stacks of monolayer MoS2 with gap-buffers is negligible due to the fast decaying LSPR mode field.

We should mention that a nanodisk with 4 nm diameter with an aspect ratio of 22 approaches the material limit of a single atom layer silver disk. In this limit, one might suspect the validity of the classical approach (COMSOL) as well as manufacturability of such thin nanoparticles. However, several experimental results on 2D single-layer silver and gold nanoparticles for bio-imaging applications were reported [40, 41, 42]. The 2D silver and gold sheets were fabricated using a self-assembly technique on a water-air interface. Such top-down fabrication approach may place the nanoparticle precisely in the desired location, but the distance among nanoparticles is challenging to be smaller than 100 nm, while other bottom-up (growth) fabrication technique may randomly place the sheet nanoparticles, but with much higher nanoparticle density [40]. Moreover, a good agreement between experimental results and theoretical analysis based on classical Maxwell equations was reported [40, 41]. This implies that the theoretical extrapolation towards such thin nanoparticle geometry based on the classical approach may still be acceptable. Moreover, the FOM has a long tail for larger diameter nanodisks. For example, a nanodisk with 10 nm diameter that has $\sim$ 0.5 nm thickness still has FOM of 2.3. We note that even larger (thicker) nanodisks still maintain FOM values larger than unity.

Nonlinear behavior of the system can be clearly seen through the population distribution among energy eigenstates when driven by an external field. For this, we solved the dynamic Equation (5) by setting $\sqrt{\kappa}\varepsilon=10$ meV driving on resonance ($\omega=\omega_{a}$). The steady-state population distributions are shown in Fig. 3. The linear system ($\hbar\chi=\hbar\chi^{\prime}=0$ meV) shows a typical coherent state population distribution, while the optimal nonlinear system (4 nm diameter, $\hbar\chi,\hbar\chi^{\prime},\hbar\gamma=53,13,21$ meV) has suppressed population in the upper levels. Even the simulation with a larger nanodisk (10 nm diameter, $\hbar\chi,\hbar\chi^{\prime},\hbar\gamma=28,16,21$ meV) shows an excellent suppression of the higher energy level populations, although, for this case, $\rho_{33}$ is nonzero. The population suppression agrees with the above analysis (Equation (6)). We adjusted the drive power slightly to normalize $\rho_{11}$ since the same drive power will lead to slightly different $\rho_{11}$ populations among the three systems. The residual population $\rho_{22}$ in the nonlinear systems is caused by the saturation of $\rho_{11}$, which could be effectively further reduced using a smaller drive power.

The quantum behavior of the nonlinear system can be clearly seen through the two-photon correlation $g^{(2)}(\tau)=\lim_{t\rightarrow\infty}\langle a^{\dagger}(t)a^{\dagger}(t+\tau)a(t+\tau)a(t)\rangle/(\langle a^{\dagger}(t)a(t)\rangle\langle a^{\dagger}(t+\tau)a(t+\tau)\rangle)$ [43]. A linear system such as an empty optical cavity will produce $g^{(2)}(\tau)=1$ for all $\tau$ when driven by a classical source, which corresponds to a coherent state. However, a strongly coupled atom–cavity system exhibits photon blockade and, thus, $g^{(2)}(0)$ close to zero [43, 44]. We obtained the two-photon correlation $g^{(2)}(\tau)$ as shown in Fig. 4 by solving the dynamic equation (Equation (5)). As expected, a linear system ($\chi=\chi^{\prime}=0$) shows a flat $g^{(2)}(\tau)=1$ for all $\tau$, regardless of the drive power. In contrast, the nonlinear system (4 nm nanodisk) shows a strong nonlinearity having $g^{(2)}(0)<1$. Obviously $g^{(2)}(0)$ depends on the population in the upper excited levels and, hence, it depends on the drive power. When weakly driven, the $g^{(2)}(0)$ value reaches as low as 0.15, predicting strong antibunching in optical fields escaping from the LSPR. As the drive power increases, the upper level starts being populated, degrading $g^{(2)}(0)$ as shown. We note that even the larger 10 nm nanodisk with less FOM also shows a strong antibunching behavior because of the FOM value above unity.

Finally, we investigate the non-Gaussianity of the states. The non-Gaussianity (nG) is crucial for a number of continuous-variable (CV) quantum information theory [45], such as the continuous-variable (CV) entanglement distillation protocol, CV quantum error correction, and cluster-state quantum computation. One of the measures for the non-Gaussianity of a state is the relative entropy measure [46, 45]

[TABLE]

where $H(x)=(x+1/2)\ln(x+1/2)-(x-1/2)\ln(x-1/2)$, $\sigma$ is the covariance matrix of the state $\rho$ using the quadrature operators $q=(1/\sqrt{2})(a+a^{\dagger})$, $p=(1/i\sqrt{2})(a-a^{\dagger})$ with the boson annihilation operator $a$, and $S(\rho)=-\mathrm{Tr}[\rho\ln\rho]$ is the von Neumann entropy. Here, $\tau$ is a Gaussian state having the same $\langle p\rangle,\langle q\rangle$ and the same $\sigma$ as the state $\rho$. While this nG measure does not immediately reveal the usefulness of the state, (which indeed requires a careful consideration of the particular CV protocol of interest,) this nG measure clearly shows how far the state is from a Gaussian state, and hence, the state with a larger nG measure $\delta[\rho]$ has a potential to be useful.

It is known that the nG measure $\delta[\rho]=0$ for any general Gaussian state $\rho=D(\alpha)S(\zeta)\nu(n)S^{\dagger}(\zeta)D^{\dagger}(\alpha)$ with arbitrary displacement operator $D(\alpha)$ and squeezing operator $S(\zeta)$ on the thermal state $\nu(n)$ with the average photon number $n$ [46]. It is also known that $\delta[\rho]$ is non-negative, and zero only for a Gaussian state $\rho$ including $\rho=\left|0\rangle\langle 0\right|$ (vacuum) [46]. As a reference, we note that $\delta[\rho]=1.39$ for $\rho=\left|1\rangle\langle 1\right|$ (single photon Fock state).

Figure 5 shows the simulated results of nG measure for the 4 nm nanodisk and 10 nm nanodisk cases, respectively, as a function of time. The initial state is vacuum and the pumping starts at time $t=0$ with the pump power of $\sqrt{\kappa}\varepsilon=27$ meV. Both 4 nm and 10 nm nanodisks show substantial nG. Particularly, the nG peaks in the first Rabi oscillation cycle, reaching 0.42 (0.24) at time $t=1.3/\gamma(t=1.6/\gamma)$ for 4 nm (10 nm) nanodisk case, respectively. The transient state with the peak nG shows a strong quantum feature of negative Wigner function values, although the purity of the state degraded. The purity at the peak nG is 0.79 (0.82) for 4 nm (10 nm) nanodisks. Such purity degradation of the state is expected in such a dissipative system. Then, one can control the pump to last only up to the point of the nG peak, which will prepare the state having strong quantum signatures. Such state with negative Wigner function values and meaningful nG values may become quite useful for CV quantum information applications.

5 Conclusion and discussions

We have demonstrated a nonlinear oscillator composed of an LSPR with added Kerr nonlinearity via coupling to monolayer MoS2. To qualify as a quantum nonlinear oscillator, the nonlinear FOM $4\chi^{2}/(\gamma+\chi^{\prime})^{2}$ must be larger than unity. The nonlinearity of a state-of-the-art smallest photonic-crystal point defect cavity filled with silicon is inadequate (FOM $\sim 10^{-6}$), whereas our proposed system has predicted FOM $\sim 10$, which is a seven order-of-magnitude improvement. The improvement stems from the strong field enhancement provided by a metal nanoparticle LSPR and direct coupling to a highly nonlinear monolayer TMD material, which is thin enough to overlap well with the tightly confined LSPR evanescent field. (Although graphene possesses an order-of-magnitude larger real $\chi^{(3)}$ at the same optical frequency than monolayer TMDs, the unfavorable real-to-imaginary $\chi^{(3)}$ ratio of graphene spoils the potential advantage [5], making graphene a less desirable material than the monolayer TMDs.) Moreover, by appropriately detuning the drive frequency, one can maximize the nonlinear FOM thanks to the discrete nonlinear susceptibility spectrum of the monolayer TMDs. It is possible that the practically achievable LSPR loss rate $\gamma$ could be larger than the theoretical value of 21 meV that we used in our calculation, for example, due to surface impurities, and in this case, the FOM would degrade accordingly. Even in the weakly nonlinear regime, however, pulse sequences analogous to those designed for suppression of population outside the qubit subspace in microwave anharmonic oscillators [47] could be used to enable manipulation of the LSPR state in an effectively nonlinear fashion.

Our proposed MNP-TMD quantum nonlinear oscillator could be further coupled to an optical cavity, facilitating its use to construct nonlinear quantum photonic devices. Recent studies demonstrated a decoherence-free operation of a system combining an LSPR mode and an emitter [19, 48]. For example, one can couple our quantum nonlinear oscillator with a cavity. In this case, one can consider the following Hamiltonian

[TABLE]

where $b$, $c$ are the annihilation operators for the nonlinear oscillator mode and the cavity mode, respectively, and $g$ is the coupling coefficient between the cavity and the nonlinear oscillator. The nature of the coupling is through the electron multipole interaction between the cavity and the nonlinear oscillator. The Lindblad operators are $L_{1}=\sqrt{\gamma_{b}}b$, the (Ohmic) dissipation of the LSPR, $L_{2}=\sqrt{\chi^{\prime}/2}b^{2}$, the two-photon absorption from the monolayer MoS2, and $L_{3}$ = $\sqrt{\gamma_{c}}c$, the cavity dissipation. Although the detailed description of this advanced system is beyond the scope of the current manuscript, one can use the adjoint master equation [17] to obtain the following Heisenberg dynamical equation for the operators:

[TABLE]

It is noteworthy that, in our preliminary study of such systems, the dissipated cavity mode $\sqrt{\gamma_{c}}c$ inherited the non-Gaussianity and the antibunching properties of the LSPR field substantially when $\gamma_{c}$ is comparably large as $\gamma_{b}$. Then, one can extract the nonlinear quantum behavior of the LSPR efficiently into a traveling wave, which is essential for a potential large-scale quantum information processing network.

Acknowledgment

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International, Inc., for the DOE’s National Nuclear Security Administration under contract DE-NA0003525.

Appendix A Field quantization in a dielectric cavity weakly coupled with a 2D nonlinear material

We follow the treatment in Drummond and Hillery [21]. We will use the displacement field $\bm{D}$ instead of the electric field $\bm{E}$ and make the quantization consistent with the Maxwell equation. The Maxwell equation for a nonmagnetic material is given as

[TABLE]

where $\bm{D}=\epsilon_{0}\bm{E}+\bm{P}$ is the displacement field, $\bm{E}$ is the electric field, $\bm{P}$ is the polarization, $\bm{B}=\mu\bm{H}$ is the magnetic field. Here, $\epsilon_{0}$ is the vacuum permittivity and $\mu$ is the magnetic permeability. The linear polarization is given by the following relation:

[TABLE]

where $j,k=x,y,z$ and $\chi^{(1)}_{jk}$ is the linear susceptibility tensor element of the medium.

The nonlinear polarization is typically defined by

[TABLE]

where $[\bm{\chi}^{(n)}]_{jklm\cdots}=\chi^{(n)}_{jklm\cdots}$. Hence, we have

[TABLE]

One can express the displacement field as

[TABLE]

where $\bm{\epsilon}$ is the material’s electric permittivity tensor.

We next introduce equivalent formalism based on the expansion of the electric field with respect to the polynomials of the displacement field:

[TABLE]

where $\bm{\eta}^{(n)}$ is the inverse permittivity tensor. The relation between $\bm{\epsilon}$ and $\bm{\eta}$ is found by equating

[TABLE]

In the case of the 2D TMD material, when we use the polarization basis such that $\hat{\varepsilon}_{+},\hat{\varepsilon}_{-},\hat{\varepsilon}_{z}$ where the first two are the clockwise and the counter-clockwise circular polarization in the x-y plane, the dielectric tensor is indeed diagonal [6]. Then, the displacement field is given by

[TABLE]

where $\chi^{(1)}=\chi^{(1)}_{+}=\chi^{(1)}_{-}$ and $\chi^{(3)}=\chi^{(3)}_{++++}=\chi^{(3)}_{----}$. Here, the factor 3 is the degeneracy such that $\chi^{(3)}=\chi^{(3)}(\omega;\omega,\omega,-\omega)=\chi^{(3)}(\omega;\omega,-\omega,\omega)=\chi^{(3)}(\omega;-\omega,\omega,\omega)$. From this, we have

[TABLE]

while all the other elements up to the third order tensors are zero.

The Hamiltonian is given as [21]

[TABLE]

For a a linear medium, the quantized displacement field is obtained as [21]:

[TABLE]

where the mode function $\bm{f}(\bm{r})$ is normalized such that

[TABLE]

For example, the electric field in a 1D Fabry-Perot linear isotropic dielectric cavity is quantized through

[TABLE]

where $V$ is the cavity volume.

Next, we consider a case where a small piece of nonlinear material is coupled to the cavity field. When the nonlinear material is substantially occupying the cavity volume, one must follow the procedure in Drummond and Hillery [21] where a dual potential $\bm{\Lambda}$ such that $\bm{D}=\bm{\nabla}\times\bm{\Lambda}$ is introduced and the quantized field is not directly $\bm{D}$, but this dual potential such that

[TABLE]

where $\bm{u}(\bm{r})$ is the mode function of the dual potential and $a,a^{\dagger}$ are the annihilation and the creation operators of the quantized field. Then, $\bm{D}$ is obtained by taking the curl of this dual potential. However, when the overlap between the cavity field and the nonlinear material is sufficiently small, one can avoid such cumbersome treatment by taking the leading order perturbing term of the Hamiltonian such that

[TABLE]

where

[TABLE]

where $m$ is the leading order of the perturbation. Then, the first term is converted to $H_{0}=\hbar\omega a^{\dagger}a$ using the quantized field $\bm{D}$ in Equation (21) and the field annihilation and creation operators.

We note that the same approach was adopted in many earlier papers regarding the Kerr-nonlinear optical cavity [2, 20, 11, 22]. In this approximate approach, $H_{0}$ represents the unperturbed Hamiltonian and $H_{\mathrm{nl}}$ is the perturbing Hamiltonian. This approach is valid when $E_{\mathrm{nl}}\ll E_{\mathrm{0}}$ where $E_{i}=\langle H_{i}\rangle$ with $i=\mathrm{0},\mathrm{nl}$. This condition is usually met when the nonlinear 2D material’s overlap with the cavity field is sufficiently small. The advantage of this approach is that the unperturbed Hamiltonian is easily quantized while the perturbing nonlinear term can be expanded using the quantized field $\bm{D}$ from the unperturbed Hamiltonian, which is given by Equation (21).

Appendix B Nonlinear coefficient $\chi,\chi^{\prime}$ and unperturbed Hamiltonian

The nonlinear coefficient $\chi_{b}$ appearing in the Hamiltonian depends on the details such as the shape, the size, and the material properties of the MNP and the nonlinear material. Once they are known, the value of $\chi_{b}$ for a particular LSPR mode can be accurately calculated. We solve $\chi_{b}$ for an arbitrary quantized field of an MNP.

The second-quantized displacement field of the LSPR mode in an arbitrary MNP is [8, 19]:

[TABLE]

where $\omega_{a}$ is the resonance frequency of the LSPR, $a$ and $a^{\dagger}$ are the annihilation and the creation operators of the LSPR with the usual bosonic commutator $[a,a^{\dagger}]=1$, and the mode function $\bm{f}(\bm{r})$ is normalized according to Equation (22). This form is consistent with the quantized $\bm{D}$ field appearing in Equation (21).

The nonlinear optical material, namely MoS2 in our case, is embedded into a background material such as the hexagonal boron nitride that fills the space. The refractive index of MoS2 at 1.06 eV is calculated through DFT method [23] to be 2.2. The refractive index of the hexagonal boron nitride is also found to be 2.2 [24]. Therefore, for simplicity, we assume that the linear dielectric function of the external space of the MNP is a constant with respect to the position $\bm{r}$. We adopt the approximate perturbative approach in the previous section to consider an LSPR mode with a weakly coupled nonlinear monolayer TMD material. The perturbation is expected to be small due to the small spatial overlap of the LSPR field with the atomically thin monolayer.

Let us now consider the nonlinear perturbation $H_{\mathrm{nl}}$ in Equation (26). The leading order is the third-order nonlinearity [6] and, therefore, we obtain

[TABLE]

where $\bm{E}_{\parallel}$ is the electric field component parallel to the 2D surface. We also explicitly noted that $\chi^{(3)}$ is a function of space so that only the in-plane electric field component in the nonlinear 2D material contributes to the third-order optical nonlinearity. The quantized in-plane component of the electric field of the unperturbed Hamiltonian is given from the quantized unperturbed displacement field $\bm{D}$ as

[TABLE]

where $\bm{f}_{\parallel}(\bm{r})$ is the in-plane component of $\bm{f}(\bm{r})$ in the 2D material. Applying the rotating wave approximation, we calculate

[TABLE]

We then obtain the nonlinear part

[TABLE]

where the first term is an uninsteresting constant and the second term is the Lamb shift from the Kerr process. Here, $\bm{f}_{\parallel}(\bm{r})$ is the in-plane component of the LSPR mode function. Therefore, we obtain

[TABLE]

Since $\chi$ and $\chi^{\prime}$ are related through the Kramers-Kronig relation, the above leads to the result in the main text (Eq. (3)). This is indeed a general result regardless of the details of the LSPR field distribution $\bm{f}(\bm{r})$. It also applies to any optical cavity.

It is straightforward to derive the equation (1) of the main text as follows. The quantized electric field in a 1D Fabry-Perot cavity is given from Equation (21) as

[TABLE]

where $\omega_{c}$ is the cavity’s mode’s frequency and $c$ is the cavity’s annihilation operator. Hence, we obtain

[TABLE]

Hence, one can calculate $\chi$ using this through the Equation (32).

Finally, we discuss the unperturbed Hamiltonian of the combined system. The nonlinear term is an addition to the linear harmonic oscillator Hamiltonian $\hbar a^{\dagger}a$ with the annihilation operator $a$ of the quantized LSPR field [8], so that the unperturbed Hamiltonian is

[TABLE]

which coincides with Eq. (1) in the main text.

Appendix C Population of the LSPR levels

To study the population of the anharmonic ladder of the LSPR, we consider the anharmonic energy eigenstates of the LSPR mode: particularly three lowest levels $\left|0\right\rangle_{b},\left|1\right\rangle_{b},\left|2\right\rangle_{b}$ whose energy eigenvalues are $0,\hbar\omega_{b},\hbar(2\omega_{b}+\chi_{b})$. An external field drives the transition among the levels. We are particularly interested in the ratio of the population at the second level to that at the first level. The dynamics of this anharmonic system can be modeled using a Hamiltonian $H^{\prime}=H^{\prime}_{a}+H^{\prime}_{\mathrm{drive}}$ with

[TABLE]

and the Lindblad dissipation operators:

[TABLE]

The master equations for the density matrix elements are obtained by [17]:

[TABLE]

For simplicity, we consider the case where $\omega=\omega_{a}$ (resonant drive) and $\varepsilon^{*}=\varepsilon$. When we separately treat the dynamics only between $\left|1\right\rangle$ and $\left|2\right\rangle$ levels and transform onto a rotating frame, the equations of motion are

[TABLE]

where we used the fact $\rho_{20}=0$. Here, we denoted $\rho_{ij}=\left\langle i\right|\rho\left|j\right\rangle$. One can solve the above equation using the Fourier transformation technique:

[TABLE]

From the second and the third equations, we obtain

[TABLE]

From this, we obtain the steady-state solution by setting $\omega^{\prime}\rightarrow 0$:

[TABLE]

When $\kappa\varepsilon^{2}$ is large, the population in $\left|2\right\rangle$ becomes saturated. For a sufficiently small drive $\varepsilon$ such that $4\kappa\varepsilon^{2}\ll{\gamma}^{2}$, we obtain the steady-state population

[TABLE]

Appendix D An example design of a metal nanoparticle

In this section, we illustrate one possible design of the metal nanoparticle and the associated LSPR mode that provides a sufficiently large Kerr coefficient in the Hamiltonian at our desired driving optical frequency. We consider a configuration in the figure 1 of the main text where an LSPR field is coupled to the two 2D material planes in the left and the right. One drives the LSPR resonance using a classical light source (lasers) as shown.

The mode field function and the resonance frequency for the spherical MNPs are well known (for example, see [9]). Unfortunately, however, the spherical MNP does not provide a large overlap between the in-plane LSPR field component and the sandwiching 2D material planes: To excite the in-plane field component, the electrical field component of the driving field must be parallel to the 2D planes. Then, most of the in-plane field component of spherical LSPR is concentrated in the middle of the sphere [8] (also see the FEM solution in the figure 6). Therefore, we must find another geometry of the MNP that provides a large overlap. It is also well known that the detailed LSPR field function (distribution) and the resonance depend largely on the shape and size of the MNP [32, 33]. A reasonable path forward is to utilize the finite element method to solve the Maxwell equations in details. Hence, we adopt the commercial FEM software package (COMSOL) to design a suitable MNP that provides the resonance peak at the desired optical frequency at 1.06 eV and maximizes the nonlinear overlap.

A natural selection of the MNP geometry that may improve the in-plane component overlapping the nonlinear 2D material is a disk shape; the large bottom and top surfaces induce large in-plane fields, while the thin layer of the disk may help make the resonant frequency lower. For example, we solved the field distribution of an LSPR mode from a silver disk MNP with diameter 6 nm and thickness 0.27 nm at 1140 nm wavelength, which is the resonant frequency of this MNP design. For this, we applied a $z$-polarized external field, traveling in $y-$ direction with 1 V/m amplitude to the system and obtained the induced electric field (subtracting the background). Figure 7 shows the result. The two monolayers of MoS2 overlap well with the in-plane component of the LSPR induced field. The E-field vectors in the monolayer planes are significant, compared to the case of the spherical MNP.