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

Daniel B. S. Soh1,2, Ryotatsu Yanagimoto1, Eric Chatterjee1, and Hideo
Mabuchi1 1E. L. Ginzton Laboratory, Stanford University, Stanford, California
94305, USA.  2Sandia National Laboratories, Livermore, California 94550, USA.
<EMAIL_ADDRESS>
###### 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.

  * 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.

Figure 1: Left: schematic of a metal nanoparticle (MNP) and 2D material.
Right: energy levels of the coupled system.

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):

$\displaystyle H_{a}$
$\displaystyle=\hbar\omega_{a}a^{\dagger}a+\frac{\hbar\chi}{2}a^{\dagger}a^{\dagger}aa,$
(1) $\displaystyle H_{\mathrm{drive}}$
$\displaystyle=-i\hbar\sqrt{\kappa}(a^{\dagger}\varepsilon e^{-i\omega
t}-a\varepsilon^{*}e^{i\omega t}),$ (2)

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):

$\chi+i\chi^{\prime}=\frac{27\epsilon_{0}\hbar\omega_{a}^{2}}{4}\int
d^{3}r\;\chi^{(3)}(\bm{r})|\bm{f}_{\parallel}(\bm{r})|^{4},$ (3)

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]:

$\dot{Q}=-\frac{i}{\hbar}[Q,H]+\sum_{j}\left(L_{j}^{\dagger}QL_{j}-\frac{1}{2}L_{j}^{\dagger}L_{j}Q-\frac{1}{2}QL_{j}^{\dagger}L_{j}\right).$
(4)

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

$\dot{a}=-i\Delta_{a}a-i\chi(a^{\dagger}a)a-\frac{\gamma}{2}a-\frac{\chi^{\prime}}{2}(a^{\dagger}a)a-\sqrt{\kappa}\varepsilon,$
(5)

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):

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

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.

Figure 2: Nonlinear FOM dependence on the diameter of a silver disk MNP and
the drive frequency detuning. Solid red lines represent the points with unity
FOM (onset of nonlinearity). Inset: geometry of MNP-TMD nonlinear system (not
in scale, only the lower layer of monolayer MoS2 is shown for clarity).

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.

(a) Linear system

(b) Nonlinear system ($4$ nm)

(c) Nonlinear system ($10$ nm)

Figure 3: Comparison of the steady-state excited level populations when driven
by an external field at resonance between a linear and a nonlinear system. The
middle nonlinear system is from a 4 nm diameter nanodisk ($\hbar\chi=53$ meV,
$\hbar\chi^{\prime}=13$ meV) while the right is from a 10 nm diameter nanodisk
($\hbar\chi=28$ meV, $\hbar\chi^{\prime}=16$ meV).

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.

(a) 4 nm nanodisk

(b) 10 nm nanodisk

Figure 4: Two-photon correlation $g^{(2)}(\tau)$ of linear (L) and nonlinear
systems (NL) with various drive powers. Left: a 4 nm nanodisk ($\hbar\chi=53$
meV, $\hbar\chi^{\prime}=13$ meV). Right: a 10 nm nanodisk ($\hbar\chi=28$
meV, $\hbar\chi^{\prime}=16$ meV).

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]

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

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).

(a) 4 nm nanodisk

(b) 10 nm nanodisk

(c) 4 nm nanodisk

(d) 10 nm nanodisk

Figure 5: Upper: transient evolution of the nG measure $\delta[\rho(t)]$ and
the purity $\mathrm{Tr}[\rho(t)]$, lower: Wigner function of the transient
states with the peak nG values, for a 4 nm nanodisk ($\hbar\chi=53$ meV,
$\hbar\chi^{\prime}=13$ meV) and a 10 nm nanodisk ($\hbar\chi=28$ meV,
$\hbar\chi^{\prime}=16$ meV).

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

$\displaystyle H$ $\displaystyle=H_{b}+H_{c}+H_{I}+H_{\mathrm{drive}},$
$\displaystyle H_{b}$
$\displaystyle=\hbar\omega_{b}b^{\dagger}b+\hbar(\chi_{b}/2)(b^{\dagger})^{2}b^{2},$
$\displaystyle H_{c}$ $\displaystyle=\hbar\omega_{c}c^{\dagger}c,$
$\displaystyle H_{I}$ $\displaystyle=-g(b^{\dagger}c+bc^{\dagger}),$
$\displaystyle H_{\mathrm{drive}}$
$\displaystyle=-i\hbar\sqrt{\kappa}(b^{\dagger}\varepsilon e^{-i\omega
t}-b\varepsilon^{*}e^{i\omega t}),$ (8)

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:

$\displaystyle\dot{b}$
$\displaystyle=-i\Delta_{b}b-i\chi_{b}(b^{\dagger}b)b+igc-\frac{\gamma_{b}}{2}b-\frac{\chi^{\prime}_{b}}{2}(b^{\dagger}b)b+\sqrt{\kappa}\varepsilon,$
$\displaystyle\dot{c}$
$\displaystyle=-i\Delta_{c}c+igb-\frac{\gamma_{c}}{2}c.$ (9)

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

$\displaystyle\bm{\nabla}\cdot\bm{D}=0,\quad\bm{\nabla}\cdot\bm{B}=0,$ (10)
$\displaystyle\bm{\nabla}\times\bm{E}=-\frac{\partial\bm{B}}{\partial
t},\quad\bm{\nabla}\times\bm{H}=\frac{\partial\bm{D}}{\partial t},$ (11)

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:

$P_{j}=\epsilon_{0}\sum_{k=1}^{3}\chi_{jk}^{(1)}E_{k},$ (12)

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

$\bm{P}=\epsilon_{0}\sum_{n>0}\bm{\chi}^{(n)}:\bm{E}^{\otimes n},$ (13)

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

$P_{j}=\epsilon_{0}\left[\sum_{k}\chi^{(1)}_{jk}E_{k}+\sum_{k,l}\chi^{(2)}_{jkl}E_{k}E_{l}+\sum_{k,l,m}\chi^{(3)}_{jklm}E_{k}E_{l}E_{m}+\cdots\right].$
(14)

One can express the displacement field as

$\bm{D}=\sum_{n>0}\bm{\epsilon}^{(n)}:\bm{E}^{\otimes n},$ (15)

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:

$\bm{E}=\sum_{n>0}\bm{\eta}^{(n)}:\bm{D}^{\otimes n},$ (16)

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

$\bm{E}=\sum_{n>0}\bm{\eta}^{(n)}:\left[\sum_{m>0}\bm{\epsilon}^{(m)}:\bm{E}^{\otimes
m}\right]^{\otimes n}.$ (17)

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

$\displaystyle D_{+}$
$\displaystyle=\epsilon_{0}\left((1+\chi^{(1)})E_{+}+3\chi^{(3)}E_{+}^{3}\right),$
$\displaystyle D_{-}$
$\displaystyle=\epsilon_{0}\left((1+\chi^{(1)})E_{-}+3\chi^{(3)}E_{-}^{3}\right),$
(18)

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

$\displaystyle\epsilon^{(1)}_{++}=\epsilon^{(1)}_{--}=\epsilon_{0}(1+\chi^{(1)}),\quad\epsilon^{(1)}_{z}=\epsilon_{0},\quad\epsilon^{(3)}_{++++}=\epsilon^{(3)}_{----}=3\epsilon_{0}\chi^{(3)},$
(19)

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

The Hamiltonian is given as [21]

$\displaystyle H$ $\displaystyle=\int
d^{3}r\left[\frac{1}{2\mu}|\bm{B}|^{2}+\sum_{n\geq
1}\frac{1}{n+1}\bm{D}\cdot(\bm{\eta}^{(n)}:\bm{D}^{\otimes n})\right]$ (20)
$\displaystyle=\int d^{3}r\left[\frac{1}{2\mu}|\bm{B}|^{2}+\sum_{n\geq
1}\frac{n}{n+1}\bm{E}\cdot(\bm{\epsilon}^{(n)}:\bm{E}^{\otimes n})\right].$

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

$\bm{D}=i\epsilon(\bm{r})\sqrt{\frac{\hbar\omega}{2}}(a\bm{f}(\bm{r})-a^{\dagger}\bm{f}^{*}(\bm{r})),$
(21)

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

$\int d^{3}r~{}\epsilon(\bm{r})\bm{f}^{*}(\bm{r})\cdot\bm{f}(\bm{r})=1.$ (22)

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

$\bm{E}=\frac{1}{\epsilon(\bm{r})}\bm{D}=i\sqrt{\frac{\hbar\omega}{2\epsilon(\bm{r})V}}(ae^{i\bm{k}\cdot\bm{z}}-a^{\dagger}e^{-i\bm{k}\cdot\bm{z}}),$
(23)

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

$\bm{\Lambda}=\sqrt{\frac{\hbar}{2\mu\omega}}(a\bm{u}(\bm{r})-a^{\dagger}\bm{u}^{*}(\bm{r})),$
(24)

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

$H=H_{0}+H_{\mathrm{nl}},$ (25)

where

$\displaystyle H_{0}$ $\displaystyle=\int
d^{3}r\left[\frac{1}{2\mu}|\bm{B}|^{2}+\frac{1}{2}\bm{D}\cdot(\bm{\eta}^{(1)}:\bm{D})\right]$
$\displaystyle=\int
d^{3}r\left[\frac{1}{2\mu}|\bm{B}|^{2}+\frac{1}{2}\bm{E}\cdot(\bm{\epsilon}^{(1)}:\bm{E})\right],$
$\displaystyle H_{\mathrm{nl}}$ $\displaystyle=\int
d^{3}r\frac{m}{m+1}\bm{E}\cdot(\bm{\epsilon}^{(m)}:\bm{E}^{\otimes m}),$ (26)

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]:

$\bm{D}(\bm{r},t)=i\epsilon(\bm{r})\sqrt{\frac{\hbar\omega_{a}}{2}}\bm{f}(\bm{r})(a(t)-a^{\dagger}(t)),$
(27)

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

$H_{\mathrm{nl}}=\int
d^{3}r\frac{9}{4}\epsilon_{0}\mathrm{Re}[\chi^{(3)}(\bm{r})]|\bm{E}_{\parallel}|^{4},$
(28)

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

$\bm{E}_{\parallel}(\bm{r})=i\sqrt{\frac{\hbar\omega_{a}}{2}}\bm{f}_{\parallel}(\bm{r})(a-a^{\dagger}),$
(29)

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

$(a-a^{\dagger})^{4}=6(a^{\dagger})^{2}a^{2}+12a^{\dagger}a+3.$ (30)

We then obtain the nonlinear part

$\displaystyle H_{\mathrm{mnp,Kerr}}$
$\displaystyle=\frac{27\epsilon_{0}\hbar^{2}\omega_{a}^{2}}{16}(2(a^{\dagger})^{2}a^{2}+4a^{\dagger}a+1)\int
d^{3}r\mathrm{Re}[\chi^{(3)}(\bm{r})]|\bm{f}_{\parallel}(\bm{r})|^{4}$ (31)
$\displaystyle=E_{0}+\hbar\Delta_{k}^{\prime}a^{\dagger}a+\hbar\frac{\chi}{2}a^{\dagger}a^{\dagger}aa,$

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

$\chi=\frac{27\epsilon_{0}\hbar\omega_{a}^{2}}{4}\int
d^{3}r\mathrm{Re}[\chi^{(3)}(\bm{r})]|\bm{f}_{\parallel}(\bm{r})|^{4}.$ (32)

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

$\bm{E}(z,t)=i\sqrt{\frac{\hbar\omega_{c}}{2\epsilon
V_{c}}}\sin(n(\omega_{c}/c_{0})z)\hat{x}\left(ce^{-i\omega_{c}t}-c^{\dagger}e^{i\omega_{c}t}\right),$
(33)

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

$\bm{f}(\bm{r})=\frac{1}{\sqrt{\epsilon
V_{c}}}\sin(n(\omega_{c}/c_{0})z)\hat{x}.$ (34)

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

$H_{a}=\hbar\omega_{a}a^{\dagger}a+\frac{\hbar\chi}{2}a^{\dagger}a^{\dagger}aa,$
(35)

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

$\displaystyle H^{\prime}_{a}$
$\displaystyle=\hbar\omega_{a}\left|1\rangle\langle
1\right|+\hbar(2\omega_{a}+\chi)\left|2\rangle\langle 2\right|,$
$\displaystyle H^{\prime}_{\mathrm{drive}}$
$\displaystyle=-i\hbar\sqrt{\kappa}\left[(\left|1\rangle\langle
0\right|+\left|2\rangle\langle 1\right|)\varepsilon e^{-i\omega
t}-(\left|0\rangle\langle 1\right|+\left|1\rangle\langle
2\right|)\varepsilon^{*}e^{i\omega t}\right],$ (36)

and the Lindblad dissipation operators:

$L_{1}=\sqrt{\gamma}\left|0\rangle\langle 1\right|,\quad
L_{2}=\sqrt{\gamma}\left|1\rangle\langle 2\right|,\quad
L_{3}=\sqrt{\chi^{\prime}}|0\rangle\langle 2|.$ (37)

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

$\dot{\rho}=\frac{i}{\hbar}[\rho,H]+\sum_{j}\left(L_{j}\rho
L_{j}^{\dagger}-\frac{1}{2}L_{j}^{\dagger}L_{j}\rho-\frac{1}{2}\rho
L_{j}^{\dagger}L_{j}\right).$ (38)

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

$\displaystyle\dot{\rho}_{22}$
$\displaystyle=-\sqrt{\kappa}\varepsilon(\rho_{21}+\rho_{12})-(\gamma+\chi^{\prime})\rho_{22},$
$\displaystyle\dot{\rho}_{21}$
$\displaystyle=-i\chi\rho_{21}+\sqrt{\kappa}\varepsilon(\rho_{22}-\rho_{11})-\frac{\gamma+\chi^{\prime}}{2}\rho_{21},$
$\displaystyle\dot{\rho}_{12}$
$\displaystyle=+i\chi\rho_{12}+\sqrt{\kappa}\varepsilon(\rho_{22}-\rho_{11})-\frac{\gamma+\chi^{\prime}}{2}\rho_{12},$
(39)

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:

$\displaystyle-i\omega^{\prime}\tilde{\rho}_{22}(\omega^{\prime})$
$\displaystyle=-\sqrt{\kappa}\varepsilon(\tilde{\rho}_{21}(\omega^{\prime})+\tilde{\rho}_{12}(\omega^{\prime}))-(\gamma+\chi^{\prime})\tilde{\rho}_{22},$
$\displaystyle-i\omega^{\prime}\tilde{\rho}_{21}(\omega^{\prime})$
$\displaystyle=-i\chi\tilde{\rho}_{21}(\omega^{\prime})+\sqrt{\kappa}\varepsilon(\tilde{\rho}_{22}(\omega^{\prime})-\tilde{\rho}_{11}(\omega^{\prime}))-\frac{\gamma+\chi^{\prime}}{2}\tilde{\rho}_{21},$
$\displaystyle-i\omega^{\prime}\tilde{\rho}_{12}(\omega^{\prime})$
$\displaystyle=+i\chi\tilde{\rho}_{12}(\omega^{\prime})+\sqrt{\kappa}\varepsilon(\tilde{\rho}_{22}(\omega^{\prime})-\tilde{\rho}_{11}(\omega^{\prime}))-\frac{\gamma+\chi^{\prime}}{2}\tilde{\rho}_{12},$
(40)

From the second and the third equations, we obtain

$\displaystyle\tilde{\rho}_{21}(\omega^{\prime})$
$\displaystyle=\frac{\sqrt{\kappa}\varepsilon(\tilde{\rho}_{22}(\omega^{\prime})-\tilde{\rho}_{11}(\omega^{\prime}))}{i(\chi-\omega^{\prime})+(\gamma+\chi^{\prime})/2},$
$\displaystyle\tilde{\rho}_{12}(\omega^{\prime})$
$\displaystyle=\frac{\sqrt{\kappa}\varepsilon(\tilde{\rho}_{22}(\omega^{\prime})-\tilde{\rho}_{11}(\omega^{\prime}))}{-i(\chi+\omega^{\prime})+(\gamma+\chi^{\prime})/2}.$
(41)

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

$\displaystyle\tilde{\rho}_{22}(\omega^{\prime}\rightarrow
0)=\frac{4\kappa\varepsilon^{2}}{(\gamma+\chi^{\prime})^{2}+4\chi^{2}+4\kappa\varepsilon^{2}}\tilde{\rho}_{11}.$
(42)

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

$\rho_{\mathrm{22,ss}}\simeq\frac{4\kappa\varepsilon^{2}}{(\gamma+\chi^{\prime})^{2}(1+4\chi^{2}/(\gamma+\chi^{\prime})^{2})}\rho_{\mathrm{11,ss}}.$
(43)

##  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.

Figure 6: In-plane component of the LSPR field in a spherical MNP. Also shown
are the two monolayer MoS2 planes in the left and the right. (Solved in COMSOL
FEM software package)

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) Side

(b) Front

(c) E-field vectors

Figure 7: E-field distribution of a silver disk MNP driven by an external
field. Solution obtained using a finite-element-method software package
(COMSOL).

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.

## References

## References

  * [1] Mabuchi H 2012 Physical Review A 85 015806
  * [2] Imamoglu A, Schmidt H, Woods G and Deutsch M 1997 Physical Review Letters 79 1467
  * [3] Sun Z, Martinez A and Wang F 2016 Nature Photonics 10 227
  * [4] Autere A, Jussila H, Dai Y, Wang Y, Lipsanen H and Sun Z 2018 Advanced Materials 30 1705963
  * [5] Soh D B, Hamerly R and Mabuchi H 2016 Physical Review A 94 023845
  * [6] Soh D B, Rogers C, Gray D J, Chatterjee E and Mabuchi H 2018 Physical Review B 97 165111
  * [7] Nielsen M P, Shi X, Dichtl P, Maier S A and Oulton R F 2017 Science 358 1179–1181
  * [8] Waks E and Sridharan D 2010 Physical Review A 82 043845
  * [9] Shvets G and Tsukerman I 2012 Plasmonics and Plasmonic Metamaterials: Analysis and Applications vol 4 (World Scientific)
  * [10] Jensen T R, Malinsky M D, Haynes C L and Van Duyne R P 2000 The Journal of Physical Chemistry B 104 10549–10556
  * [11] Miranowicz A, Paprzycka M, Liu Y x, Bajer J and Nori F 2013 Physical Review A 87 023809
  * [12] Lin Q, Zhang J, Piredda G, Boyd R W, Fauchet P M and Agrawal G P 2007 Applied physics letters 91 021111
  * [13] Jana N R, Gearheart L and Murphy C J 2001 The Journal of Physical Chemistry B 105 4065–4067
  * [14] Rogers C, Gray D, Bogdanowicz N and Mabuchi H 2018 Physical Review Materials 2 094003
  * [15] Novoselov K, Mishchenko A, Carvalho A and Neto A C 2016 Science 353 aac9439
  * [16] Lukoševičius M and Jaeger H 2009 Computer Science Review 3 127–149
  * [17] Breuer H P, Petruccione F et al 2002 The theory of open quantum systems (Oxford University Press on Demand)
  * [18] Carmichael H J 2009 Statistical methods in quantum optics 2: Non-classical fields (Springer Science & Business Media)
  * [19] Peng P, Liu Y C, Xu D, Cao Q T, Lu G, Gong Q, Xiao Y F et al 2017 Physical review letters 119 233901
  * [20] Ferretti S and Gerace D 2012 Physical Review B 85 033303
  * [21] Drummond P D and Hillery M 2014 The quantum theory of nonlinear optics (Cambridge University Press)
  * [22] Li A, Zhou Y and Wang X B 2017 Scientific reports 7 7309
  * [23] Ullah M S, Yousuf A H B, Es-Sakhi A D and Chowdhury M H 2018 Analysis of optical and electronic properties of mos2 for optoelectronics and fet applications AIP Conference Proceedings vol 1957 (AIP Publishing) p 020001
  * [24] Cappellini G, Satta G, Palummo M and Onida G 2001 Physical Review B 64 035104
  * [25] Wang G, Robert C, Glazov M, Cadiz F, Courtade E, Amand T, Lagarde D, Taniguchi T, Watanabe K, Urbaszek B et al 2017 Physical review letters 119 047401
  * [26] Echeverry J, Urbaszek B, Amand T, Marie X and Gerber I 2016 Physical Review B 93 121107
  * [27] Kuramochi E, Kim J K, Taniyama H, Shinya A, Kita S and Notomi M 2017 Ultrahigh-q/v single point-defect photonic crystal nanocavity with embedded sub-wavelength air-slot CLEO: Science and Innovations (Optical Society of America) pp JTh3M–5
  * [28] Hales J M, Chi S H, Allen T, Benis S, Munera N, Perry J W, McMorrow D, Hagan D J and Van Stryland E W 2018 Third-order nonlinear optical coefficients of si and gaas in the near-infrared spectral region CLEO: Applications and Technology (Optical Society of America) pp JTu2A–59
  * [29] Yoshida H, Mozume T, Nishimura T and Wada O 1998 Electronics Letters 34 913–915
  * [30] Akiyama T, Georgiev N, Mozume T, Yoshida H, Gopal A V and Wada O 2001 Electronics Letters 37 129–130
  * [31] Zakery A and Elliott S 2003 Journal of Non-Crystalline Solids 330 1–12
  * [32] Link S and El-Sayed M A 2000 International reviews in physical chemistry 19 409–453
  * [33] Kelly K L, Coronado E, Zhao L L and Schatz G C 2003 The optical properties of metal nanoparticles: the influence of size, shape, and dielectric environment
  * [34] Wang F and Shen Y R 2006 Physical review letters 97 206806
  * [35] Kirakosyan A S, Stockman M I and Shahbazyan T V 2016 Physical Review B 94 155429
  * [36] Marder M P 2010 Condensed matter physics (John Wiley & Sons)
  * [37] Liu Y, Weiss N O, Duan X, Cheng H C, Huang Y and Duan X 2016 Nature Reviews Materials 1 16042
  * [38] Man M K, Deckoff-Jones S, Winchester A, Shi G, Gupta G, Mohite A D, Kar S, Kioupakis E, Talapatra S and Dani K M 2016 Scientific reports 6 20890
  * [39] Golla D, Chattrakun K, Watanabe K, Taniguchi T, LeRoy B J and Sandhu A 2013 Applied Physics Letters 102 161906
  * [40] Usukura E, Shinohara S, Okamoto K, Lim J, Char K and Tamada K 2014 Applied physics letters 104 121906
  * [41] Tanaka D, Imazu K, Sung J, Park C, Okamoto K and Tamada K 2015 Nanoscale 7 15310–15320
  * [42] Masuda S, Yanase Y, Usukura E, Ryuzaki S, Wang P, Okamoto K, Kuboki T, Kidoaki S and Tamada K 2017 Scientific reports 7 3720
  * [43] Walls D F and Milburn G J 2007 Quantum optics (Springer Science & Business Media)
  * [44] Kimble H, Polzik E, Rempe G and Thompson R 1992 Structure and dynamics in cavity quantum electrodynamics Quantum Electronics and Laser Science Conference (Optical Society of America) p QTuA2
  * [45] Genoni M G and Paris M G 2010 Physical Review A 82 052341
  * [46] Genoni M G, Paris M G and Banaszek K 2008 Physical Review A 78 060303
  * [47] Motzoi F, Gambetta J M, Rebentrost P and Wilhelm F K 2009 Physical review letters 103 110501
  * [48] Gurlek B, Sandoghdar V and Martín-Cano D 2017 ACS Photonics 5 456–461