We report on a theoretical study of second-harmonic generation (SHG) in plasmonic nanostructures interacting with two-level quantum emitters (QEs) under incoherent energy pump. We generalize the driven-dissipative Tavis–Cummings model by introducing the anharmonic surface plasmon-polariton (SPP) mode coupled to QEs and examine physical properties of corresponding SPP-QE polariton states. Our calculations of the SHG efficiency for strong QE-SPP coupling demonstrate orders of magnitude enhancement facilitated by the polariton gain. We further discuss time-domain numerical simulations of SHG in a square lattice comprising Ag nanopillars coupled to QEs utilizing a fully vectorial nonperturbative nonlinear hydrodynamic model for conduction electrons coupled to Maxwell–Bloch equations for QEs. The simulations support the idea of gain enhanced SHG and show orders of magnitude increase in the SHG efficiency as the QEs are tuned in resonance with the lattice plasmon mode and brought above the population inversion threshold by incoherent pumping. By varying pump frequency and tuning QEs to a localized plasmon mode, we demonstrate further enhancement of the SHG efficiency facilitated by strong local electric fields. The incident light polarization dependence of the SHG is examined and related to the symmetries of participating plasmon modes.

Research in optics of nanomaterials including the nanoplasmonics has significantly advanced due to both tremendous progress in nanofabrication and further characterization employing various new ultrafast spectroscopic techniques.1 The nanoplasmonics has enjoyed substantial growth scrutinizing physical properties of propagating2,3 and localized surface plasmon-polaritons (SPPs) and utilizing their unique polarization features to control light propagation4 at various metal–dielectric interfaces.5 The nonlinear optics of plasmonic nanosctructures has quickly become a hot research field.6 Applications range from optical bistability7 through second-harmonic generation (SHG)8 to difference frequency generation.9 For instance, by manipulating geometrical parameters, one can achieve substantial increase in four-wave mixing processes.10 Both experiment and theory show that the SHG by plasmonic nanoparticles with no center of the inversion symmetry is extremely sensitive to particle shape even within few nanometers, opening new opportunities to designing functional materials for frequency conversion applications.11 

On the other hand, strong spatial localization of SPPs could be used to investigate how ensembles of quantum emitters (QEs), e.g., organic dyes or semiconductor quantum dots, behave if resonantly coupled to corresponding surface electromagnetic modes.12 When in such systems the coupling strength surpasses all damping rates, the so-called strong coupling regime is achieved, resulting in formation of QE-SPP hybridize polariton states.13 These states will be referred to as the polaritons throughout this paper. Recent studies show that in the strong coupling regime, the second-order nonlinear processes that are mostly caused by highly polarizable metal electrons are significantly altered by QEs.14 The polariton states contribute on equal footing to the SHG resulting in detectable Rabi splitting of the signal even at a single nanoparticle level.15 Numerical simulations show that the presence of a gain medium in core–shell metal–dielectric nanostructures supporting localized surface-plasmon modes leads to an enhancement of SHG.16 However, mechanisms of enhanced SHG by the polariton states contributed by QEs in the population inversion regime require theoretical examination.

From quantum plasmonics point of view, highly polarizable metal nanostructures can play a role of open optical cavities.17 Corresponding plasmonic cavity modes are typically described as perfect harmonic oscillators, while the nonlinearity arises from the interaction with QEs. In a strong coupling regime, such a nonlinear system shows unique properties of the polariton states. As recently demonstrated for plasmonic lattices, these states are responsible for nonlinear emission, lasing,18–26 non-equilibrium superradiant phase transitions,27 and Bose–Einstein condensation.22,28–31 An ability of plasmonic nanostructures to support anharmonic response32 brings an interesting notion of anharmonic quantum plasmonic cavities. In this paper, we study an interplay of the QE nonlinearity and anharmonicity effects on the polariton state structure and nonlinear optical response in such plasmonic cavities.

For this purpose, we generalize the driven dissipative Tavis–Cummings model, allowing for lasing phase transitions.19,27,33 Specifically, we introduce an ensemble of two-level QEs subject to incoherent energy pump and dissipation, which interacts with a single anharmonic SPP mode. In practice, incoherent pumping can be realized via optical excitation of high electronic states of QEs followed by the excitation relaxation to the optically active states or via direct carrier injection. If the pumping rate is high enough to create QE population inversion, a gain builds up compensating plasmonic losses. Thus, we focus on examining formation of the polariton gain and subsequent enhancement of SHG in the strong QE-SPP coupling regime for various strengths of SPP anharmonicity.

Although we quantize the anharmonic SPP mode, our results are valid in the semi-classical limit holding for a high population of SPP quanta and typical for the SHG experiment. To this extent, our minimal model is capable of providing an upper limit estimate for the SHG efficiency. To account for the effect of the complex structure of SPP modes in plasmonic lattices, we further performed time-domain numerical simulations of SHG in the 2D lattice of Ag nanopillars coupled to QEs. These calculations utilize a fully vectorial nonlinear hydrodynamic model for conduction electrons coupled to Maxwell–Bloch equations describing QEs. The simulation results can be qualitatively interpreted with the help of our minimal model and give a lower boundary for experimentally expected values of SHG efficiency. On the other hand, the simulations reveal enhanced SHG dependence on the incident light polarization, which can be experimentally verified.

This paper is organized as follows: In Sec. II, we discuss the anharmonic Tavis–Cummings model and its mean-field steady states pinpointing the physics of SHG enhancement by polariton gain. In Sec. III, we discuss the results of the SHG simulations in the Ag nanopillar lattice identifying the role of various SPP modes. Conclusions are drawn in Sec. IV.

Our model illustrated in Fig. 1 consists of an array of metal nanoparticles (MNPs) supporting the delocalized SPP mode treated as an anharmonic oscillator described by a bosonic operator ψ^ with Hermitian conjugate ψ^. The SPP interacts with an ensemble of No identical two-level QEs described by spin operators ŝn±=ŝnx±iŝny and snz, expressed in terms of the Pauli SU(2) operators as ŝnj=12σ^nj with j = x, y, z and site index n=1,No¯. The QEs are subject to incoherent energy pump, whereas the SPP interacts coherently with an incident monochromatic laser field Ein(ω). By taking into account that the SPP transition dipole, μsp, exceeds the QE transition dipole by orders of magnitude, we consider the output electric field generated by the SPP only. The field is further partitioned into two modes Eout(ω)μspψ^(ω) and Eout(2ω)μspψ^(2ω) oscillating with fundamental and second-harmonic frequencies, respectively. Both QEs and MNPs interact with the environment, causing dissipative processes.

FIG. 1.

A schematic of the adopted driven-dissipative anharmonic Tavis–Cummings model.

FIG. 1.

A schematic of the adopted driven-dissipative anharmonic Tavis–Cummings model.

Close modal

The coherent dynamics of our model is described by the Hamiltonian,

H^=ωoψ^ψ^+αψ^ψ^ψ^+ψ^ψ^ψ^+ωon=1Noŝnz+No2+λn=1Noψ^ŝn+ŝn+ψ^Ωinψ^eiωt+ψ^eiωt,
(1)

represented in units of . Here, we set the SPP and all QEs to have the same transition frequency ωo and identical coupling rate λ. The second term of the Hamiltonian reflects the SPP cubic anharmonicity with the coupling rate α. Interaction between the SPP and incident driving field, Ein, is accounted for in the last term of the Hamiltonian where the coupling strength is represented by the Rabi frequency Ωin = (μsp · Ein)/.

To include the driven-dissipative dynamics, we employ the Heisenberg equation of motion in the form

tO^=iH^,O^+2γspD^ψ^[O^]+γn=1NoD^ŝn+[O^]+γn=1NoD^ŝn[O^]+γϕn=1NoD^ŝnz[O^],
(2)

where the first term in the right-hand side is responsible for coherent dynamics due to Hamiltonian H^. The rest of the terms describe the SPP decay with the rate 2γsp, incoherent pumping (population decay) of QEs with the rate γ (γ), and QE pure dephasing with the rate γϕ, respectively. The associated Lindblad operator in the Heisenberg representation reads

D^Ĉ[O^]=12ĈO^,Ĉ+12Ĉ,O^Ĉ.
(3)

In the limit of α = Ωin = 0, Eqs. (1)(2)(3) recover the open Tavis–Cummings model, allowing for the lasing phase transition.19,27,33 Our generalization of the model includes the SPP anharmonicity and external coherent laser field. As demonstrated below, this generalization constitutes a minimal model to examine SHG and the effect of the QE induced gain to enhance this process.

Starting with Eqs. (2) and (3) and the Hamiltonian (1), we derive the Heisenberg operator equations for O^={ψ^,ŝn,ŝnz}, and adopting the mean-field approximation replaces the operators with their averages. The result is

tφ=(iωo+γsp)φiNoαφ2+2|φ|2iNoλsiΩinNoeiωt,
(4)
ts=iωo+γos+2iNoλszφ,
(5)
tsz=γtszdo2+iNoλφ*ss+φ,
(6)

where we introduce normalized SPP, φ=ψ^/No and QE, s±=n=1Noŝn±/No, coherences, and normalized QE population sz=n=1Noŝnz/No varying in the interval [−1/2, 1/2]. Note that the effective QE-SPP coupling, Noλ, and anharmonicity, Noα, rates scale with the system size as No and can be enhanced by increasing the number of QEs.

The SPP dephasing rate, γsp, entering Eq. (4) is the same as in Eq. (2). However, the QE dissipation rates in Eqs. (5) and (6) become linear combinations of those in Eq. (2). Specifically, the QE dephasing and population decay rates are γo = γ/2 + γ/2 + γϕ and γt = γ + γ, respectively. Notice their dependence on the incoherent pumping rate, γ. Furthermore, Eq. (6) contains the quantity

do=γγγ+γ,
(7)

referred to as the population inversion parameter. Depending on the incoherent pumping rate, γ, it can be −1 < do < 0, indicating that γ < γ and QEs are pumped below the population inversion threshold. Increase in the pumping to γ > γ brings QEs to the inversion regime with 0 < do < 1.

To solve Eqs. (4)(5)(6), accounting for the fundamental frequency and the second-harmonic response, we transform the coherences into rotating frames as

φ(t)=φ̃1(t)eiωLt+φ̃2(t)e2iωLt,
(8)
s(t)=s̃(t)eiωLt,
(9)

where slowly varying amplitudes for the fundamental SPP, φ̃1(t), and QE, s̃(t), coherences as well as for the second-harmonic SPP coherence φ̃2(t) are introduced. Substitution of Eqs. (8) and (9) into Eqs. (4)(5)(6) and subsequent elimination of the fast oscillating terms result in the following set of coupled equations of motion for the slowly varying amplitudes,

tφ̃1=iδω+γspφ̃12iαφ̃1*φ̃2iNoλs̃+iΩinNoei(ωωo+δω)t,
(10)
tφ̃2=i[2δωωo]+γspφ̃2iNoαφ̃12,
(11)
ts̃=iδω+γos̃+2iNoλszφ̃1,
(12)
tsz=γtszdo2+iNoλφ̃1*s̃s̃+φ̃1.
(13)

Here, δω = ωoωL is the frequency detuning parameter.

Now, let us examine the steady states of Eqs. (10)(11)(12)(13), focusing on the properties of fundamental SPP coherence, φ̃1. Keeping in mind that the external field is weak enough to have a perturbative effect on the steady states, we set Ωin = 0. By solving Eqs. (11)(12)(13) for the steady states and using these solutions to eliminate the second-harmonic coherence and the QE variables in Eq. (10), we obtain the following nonlinear algebraic equation:

Σ|φ̃1|φ̃1=0,
(14)

determining the steady states of φ̃1. Here, the prefactor

Σ|φ̃1|=iΔΩ|φ̃1|+Γ|φ̃1|G|φ̃1|
(15)

contains the real and imaginary part functions of |φ̃1|. The imaginary part,

ΔΩ|φ̃1|=δω2Noα2(2δωωo)(2δωωo)2+γsp2|φ̃1|2+doNoλ2δωδω2+γo2+4Noλ2[γt/γo]|φ̃1|2,
(16)

describes renormalization of the detuning parameter by self-energies due to the anharmonic interaction (second term) and the SPP-QE coupling (third term).

The real part of Eq. (15) describes an interplay of two dephasing terms, specifically,

Γ|φ̃1|=γsp+2Noα2γsp|φ̃1|2(2δωωo)2+γsp2
(17)

being the SPP dephasing due to the energy decay into heat and far-field radiation (first term) and due to the energy transfer to the second-harmonic (second term). The other contribution to the real part is the gain rate,

G|φ̃1|=doNoλ2γoδω2+γo2+4Noλ2[γt/γo]|φ̃1|2,
(18)

compensating the dephasing term, provided 0 < do ≤ 1, i.e., an inverted QE population. Below the population inversion threshold, −1 ≤ do < 0, the sign of Eq. (18) becomes negative, resulting in additional contribution to the losses.

Equation (14) has a trivial solution φ̃1=0 and a nontrivial solution defined by Σ|φ̃1|=0. The latter is satisfied if the real part (frequency shift) and the imaginary part (gain-loss term) of Eq. (15) both become zero, i.e.,

ΔΩ|φ̃1|=0,
(19)
Γ|φ̃1|G|φ̃1|=0.
(20)

These conditions determine the lasing phase transition in the driven-dissipative Tavis–Cummings model. Roots of Eqs. (19) and (20) are coherences |φ̃1| spontaneously appearing above critical coupling, λc, and the detuning parameter δω describing the lasing frequency shift.19,27,33 Below λc, the trivial solution, φ̃1=δω=0, holds.

Having defined the steady state solutions, we proceed to examine both linear and non-linear responses induced by the incident driving field. For the linear response, we use the time-dependent equation of motion for the SPP coherence at the fundamental frequency,

t+Σ|φ̃1|φ̃1=iΩinNoei(ωωo+δω)t,
(21)

derived in the same way as Eqs. (14)(15)(16)(17)(18) but retaining the time-dependence. We further seek the solution of this equation in terms of time-dependant perturbations δφ1(t) of the mean-field steady state φ¯1. Making the substitution of φ̃1(t)=φ¯1+δφ1(t) into Eq. (21) and complementing the resulting equation with a complex conjugate, we derive linearized equations of motion in the matrix form,

It+MδΦ=iΩinNoF(t),
(22)

with the column vectors δΦ=(δφ1,δφ1*)T and F(t)=ei(ωωo+δω)t,ei(ωωo+δω)tT. Here, I is the 2 × 2 unit matrix and

M=MdMoMo*Md*,,
(23)

stands for the so-called stability matrix with matrix elements

Md=Σ|φ¯1|+φ¯1φ¯1Σ|φ¯1|,
(24)
Mo=φ¯1φ¯1*Σ|φ¯1|,
(25)

depending on Σ|φ¯1| [Eqs. (15)(16)(17)(18)] and its functional derivative both functions of the steady state solutions, φ¯1, of Eq. (14).

Solution of Eq. (22) in the presence of an oscillating incident field defines SPP polarization at the fundamental frequency,

μspδφ1(t)=χ(1)(ω)Einei(ωωo+δω)t,
(26)

where the linear susceptibility is

χ(1)(ωL)=μsp2Nou+Lu+Rωωo+δω+iΛ++uLuRωωo+δω+iΛ,
(27)

and depends on the stability matrix, M, eigenvalues

Λ±=ReMd±Mo2ImMd2,
(28)

and the upper components of orthonormal left δΦ±L=(u±L,v±L) and right δΦ±R=(u±R,v±R)T eigenvectors. Poles of the linear susceptibility,

ζ±=ωoδωiΛ±,
(29)

determine the polariton frequency, ω±=Reζ±, and dephasing rate γ±=Imζ±.

Below the critical coupling λ < λc, φ̃1=δω=0, simplifying Eq. (28) to the form

Λ±=γspdoNoλ2γo,
(30)

where the second term in the right-hand side is an explicit representation of the polariton gain term G{|φ¯1|=0} [Eqs. (18)]. According to Eqs. (29) and (30), the polariton frequency is constant ω± = ωo. In contrast, the polariton dephasing γ± = Λ± decreases with the coupling increase until the gain rate fully compensates the losses, i.e., γ± = 0. This condition determines the critical coupling,

Noλc2=γspγodo,
(31)

where γo depends on the incoherent pumping rate and do. Equation (31) is known as the lasing threshold in the open Tavis–Cummings model,19,27 which in our generalized case does not acquire dependence on the anharmonicity rate, α. According to Eq. (31), requirement that λ > λc implies that Noλ2>γspγo, i.e., the QE-SPP coupling strength exceeds the system losses. This is a definition of the strong coupling regime.31 

In the strong coupling regime, λ > λc, the polariton splitting is determined by the square root expression in Eq. (28) containing Mo2ImMd2. Depending on the incoherent pumping rate and anharmonicity values, this difference can become either positive or negative. In the former case, the square root becomes real and the polariton dephasing γ± experiences splitting into two branches, while the frequency ω± stays degenerate. In the latter case, the square root becomes imaginary, leading to splitting of the polariton energy ω±, but the dephasing γ± stays degenerate. As shown below, the former situation represents the lasing/amplification regime, whereas the later one does not.

Assuming that the anharmonicity rate αωo, λ, we perturbatively calculate the second-harmonic response, δφ2(t), using Eq. (11) with the linear solution [Eqs. (26) and (27)], δφ12, as the driving force. Accordingly, the solution of Eq. (11) gives the following SPP second-harmonic polarization:

μspδφ2(t)=χ(2)(2ω,ω,ω)Ein2e2i(ωωo+δω)t,
(32)

in terms of the second-order nonlinear susceptibility,

χ(2)(2ω,ω,ω)=Noαχ(1)(ω)2ωo2ωiγsp,
(33)

which in turn depends on χ(1)(ω) [Eq. (27)].

By evaluating the time-averaged radiated power34 due to the SPP polarization at the second harmonic [Eq. (32)] and the fundamental frequency [Eq. (26)], we introduce their ratio P(2ω)/P(ω)=16|χ(2)Ein2|2/|χ(1)Ein|2 quantifying the SHG efficiency.35 With the help of Eq. (33), this expression simplifies to the form

P(2ω)P(ω)=16Noα2|χ(1)(ω)Ein|2(2ωωo)2+γsp2,
(34)

explicitly indicating quadratic scaling of the SHG efficiency both with χ(1) and the anharmonic coupling Noα.

We perform the analysis by calculating the polariton dispersion [Eqs. (28) and (29)] and related absorption/emission spectra determined by Imχ(1)(ω) [Eq. (27)] and the SHG efficiency [Eq. (34)] depending on SPP-QE coupling strength Noλ2/ωo2. Here and below, all the rates/frequencies are represented in units of ωo. The relaxation rates and incident field Rabi frequency are parameterized as follows: γsp/ωo = 7 × 10−3, γ/ωo = 7 × 10−3, γϕ/ωo = 5 × 10−3, and Ωin/ωo = 4 × 10−3. The incoherent pumping rate is fixed at γ = 1.1γ, corresponding to the population inversion characterized by do = 0.05. For these parameters, the critical coupling value is Noλc2/ωo2=1.9×103. These parameters are similar to those adopted in the simulations discussed in Sec. III. Finally, we introduce two limiting regimes of weak, Noα/ωo=103, and strong, Noα/ωo=0.2, anharmonicity and compare them below.

Figure 2(a) shows polariton frequency (cyan curve) and the dephasing rate (red curves) as functions of the coupling strength between QEs and SPP calculated for the weak anharmonicity. The dispersion follows the trends discussed in Sec. II B: Below the strong coupling threshold, the frequency is constant, ω± = ωo, and the dephasing rate γ± linearly falls to zero [see Eq. (30)]. Above the threshold, the frequency does not show visible changes but the polariton dephasing splits into the lower and upper branches as a result of the real square root in Eq. (28). Associated linear spectra presented in panel (b) of Fig. 2 show the rise of a single absorption peak below the threshold flipping into the narrow negative dip above the threshold, i.e., emission of the amplified incident field typical for the lasing steady state. This peak is associated with the lower-polariton dephasing branch, i.e., γ ≈ 0, shown in panel (a). The upper-polariton dephasing branch contributes to the widths of a broad absorption feature around the emission dip. Maximum amplification occurs in the vicinity of the critical coupling value. Further increase in the coupling strengths reduces the depth of the dip. Panel (c) of Fig. 2 shows more than five orders of magnitude enhancement in the SHG efficiency facilitated by the amplification for the coupling value slightly about the lasing threshold (red curve). According to Eq. (34), such an enhancement rises from quadratic scaling of the linear susceptibility compensating weak anharmonic coupling.

FIG. 2.

[(a) and (d)] Polariton frequency ω± and dephasing rate γ± dispersion entering the strong coupling regime above critical coupling at Noλc2/ωo2=1.9×103. [(b) and (e)] Linear absorption/emission spectra of the polariton states. [(c) and (f)] SHG efficiency [Eq. (34)] calculated for various QE-SPP coupling strengths Noλ2/ωo2. Panels (a)–(c) present weak anharmonic coupling, Noα/ωo=103, and panels (d)–(f) present strong anharmonic coupling, Noα/ωo=0.2. Curves in panels (b) and (e) are vertically staggered for clarity with the color coded numbers indicating corrections to the base shift values.

FIG. 2.

[(a) and (d)] Polariton frequency ω± and dephasing rate γ± dispersion entering the strong coupling regime above critical coupling at Noλc2/ωo2=1.9×103. [(b) and (e)] Linear absorption/emission spectra of the polariton states. [(c) and (f)] SHG efficiency [Eq. (34)] calculated for various QE-SPP coupling strengths Noλ2/ωo2. Panels (a)–(c) present weak anharmonic coupling, Noα/ωo=103, and panels (d)–(f) present strong anharmonic coupling, Noα/ωo=0.2. Curves in panels (b) and (e) are vertically staggered for clarity with the color coded numbers indicating corrections to the base shift values.

Close modal

For the adopted strong anharmonicity value, the polariton dispersion [Fig. 2(d)] above the strong QE-SPP coupling threshold is quite different compared to the weak anharmonicity case [panel (a)]. The increase in the anharmonicity makes the expression under the square root in Eq. (28) negative, resulting in its imaginary values and subsequently in the splitting of the polariton frequency (cyan curve). However, the dephasing rate stays degenerate as pointed out in Sec. II B. In this case, the absorption spectra shown in Fig. 2(e) acquire the well known polariton splitting showing up in the strong coupling regime. In contrast to the weak anharmonicity case, the peaks do not become negative dips, which is the signature of the incident field absorption rather than amplification. The SHG efficiency shown in panel (f) of Fig. 2 demonstrates maximum three orders of magnitude enhancement (compare peak values of red and black curves), but the absolute value of the maximum efficiency (red curve) is comparable with that in the weak anharmonicity regime [red curve in panel (c)]. In this case, the efficiency gains additional increase due to its quadratic scaling with the anharmonicity parameter Noα2 [see Eq. (34)].

Finally, we fix the coupling strength slightly above the critical coupling by setting Noλ2/ωo2=2.0×103 and vary the incoherent pumping rate below the inversion threshold to sample the states characterized by do = −0.33 (γ = 0.5γ) and do = 0.05 (γ = 0.9γ). Comparison of the SHG efficiencies calculated below inversion with that for the inverted QEs (do = 0.05) is shown in panels (a) and (b) of Fig. 3 for the weak and strong anharmonicities, respectively. As expected, the inverted system in the strong coupling regime shows tremendous enhancement of the SHG efficiency. However, comparison of Fig. 3(a) with Fig. 2(c) and Fig. 3(b) with Fig. 2(f) clearly demonstrates that bringing the QEs to the population inversion but keeping the QE-SPP strength below the critical value does not result in significant SHG enhancement. Thus, to achieve high efficiency values, the system should be in the strong coupling regime.

FIG. 3.

Comparison of the SHG efficiency for various values of the inversion parameter do in the case of (a) weak, Noα/ωo=103, and (b) strong, Noα/ωo=0.2, anharmonicity. The QE-SPP coupling strength is fixed at Noλ2/ωo2=2.0×103.

FIG. 3.

Comparison of the SHG efficiency for various values of the inversion parameter do in the case of (a) weak, Noα/ωo=103, and (b) strong, Noα/ωo=0.2, anharmonicity. The QE-SPP coupling strength is fixed at Noλ2/ωo2=2.0×103.

Close modal

Here, we consider a periodic square lattice comprised of Ag nanopillars placed on top of the semi-infinite glass substrate, as schematically illustrated in Fig. 4(a). The input side of the lattice is covered by a 200 nm thick dielectric. Such a setup has been shown to support both the localized plasmon (2.1 eV) and the surface lattice plasmon seen at ℏωo = 1.96 eV in Fig. 4(b). Combined with colloidal quantum dots, this system was experimentally demonstrated to exhibit lasing.24 

FIG. 4.

(a) Schematics of a repeat unit forming two-dimensional lattice. (b) Linear absorption spectrum of the plasmonic lattice calculated without QEs (black line), with QEs but no incoherent pump, do = −1 (γ = 0) (red line), and in the inverted regime, do = 0.048 (γ = 1.1γ) (blue line).

FIG. 4.

(a) Schematics of a repeat unit forming two-dimensional lattice. (b) Linear absorption spectrum of the plasmonic lattice calculated without QEs (black line), with QEs but no incoherent pump, do = −1 (γ = 0) (red line), and in the inverted regime, do = 0.048 (γ = 1.1γ) (blue line).

Close modal

Similar to Sec. II, our computational semi-classical model treats QEs as two-level systems whose time-evolution is described by Bloch equations in the form

t2+2γot+ωo2Po=2ωoρoμoμoEsz,
(35)
tsz+γtszdo2=2ωoρotPoE,
(36)

where Po = ρoμosx denotes the QE polarization with sx = Re[s], ρo denotes QE density, and μo stands for the QE transition dipole. The relaxation rates in Eqs. (35) and (36) have the same definition as in Eqs. (5) and (6). However, in contrast to the Tavis–Cummings model, we explicitly introduce classical local electric field E interacting with QEs. This electric field is calculated by numerically solving Maxwell’s equations governing the dynamics of electromagnetic radiation including the coherent incident field. The optical response of the plasmonic lattice is modeled using the nonlinear hydrodynamic model resulting in the following equation of motion for the metal polarization, P(r, t):14 

t2P+γetP=eme*n0eE+tP×BEP1n0etPtP+tPtP.
(37)

Here, the plasma frequency is ωp=n0e2/ε0me*, the electron gas damping parameter is set to γe = 0.14 eV, and the number density of conduction electrons is n0 = 5.86 × 1028 m−3. Note that in the linear regime, Eq. (37) reduces to the conventional Drude model.

The QEs are tuned in resonance with the surface lattice plasmon mode at 1.96 eV and uniformly distributed in the dielectric surrounding the nanopillars. Coupled Maxwell–Bloch equations are propagated in the time-domain using home-build codes utilizing the finite-difference time-domain method and the weakly coupling method.36 To speed up the calculations, we employ three-dimensional domain decomposition and split the simulation domain into 1152 sub-domains, each carried out by a single processor. With the spatial resolution of 1.5 nm and dimensions summarized in Fig. 4(a), the total number of grid points is 6.7 × 107. To ensure numerical convergence, the time step is set at 2.5 as. Typical execution times of our codes vary between 20 min (for linear simulations) and 3 h (for nonlinear simulations with the hydrodynamic model employed). Other parameters for emitters are as follows: the number density is 1026 m−3, the transition dipole moment is 10 D, the pure dephasing rate is γϕ = 10−2 eV, and the non-radiative decay rate is 6 × 10−3 eV. The system is coherently driven by a plane wave linearly polarized along the x-axis and propagating in the z-direction normal to the lattice plane. Due to limited computational resources, the excitation pulse of 500 fs duration is used. For better numerical convergence of power spectra at high pump intensities, we implement the Blackman–Harris time envelope. All nonlinear results are obtained for the incident field amplitude of Ein = 2 × 10−2 V/nm.

As a measure of the linear response of the plasmonic system, we calculate absorption as a function of incident photon energy, as shown in Fig. 4(b). To simulate linear absorption, we implement the short pulse method allowing us to evaluate a wide frequency response in a single run.13 The system is excited by a short pulse, and transmission, T, and reflection, R, are calculated. The absorption, A, is then evaluated according to the following expression: A = 1 − TR. Simulations without QEs (black line) result in three distinct peaks nearly identical to those obtained in Ref. 24. The localized plasmon resonance is seen at 2.12 eV. The two other modes peaked at 1.67 eV and 1.96 eV are due to the surface lattice mode coupled to the guiding mode inside the dielectric covering nanopillars.

In the absence of an incoherent pump (γ = 0), the spectrum (red line) exhibits a large peak near QE transition frequency, indicating the response of strongly coupled QE-SPP states. Turing on the incoherent pump and setting the pumping rate to γ = 1.1γ, we reach the inversion regime characterized by do = 0.048. In this regime, the system acquires polariton gain, resulting in the negative absorption (i.e., amplification) at the QE transition frequency (blue line). Interestingly, due to strong coupling between QEs and SPP and significant material dispersion, the gain is also showing up at the nearby plasmon resonances (1.7 eV and 2.0 eV). Noteworthy is that reduction in the incoherent pumping rate in the range of 0 < do < 0.048 destroys the gain. For example, a weak inversion with do = 0.020 gives an absorption spectrum nearly equivalent to that with no pump at all. This can be understood by looking at Eq. (31) according to which and for do ≪ 1 the critical parameter determining the strong QE-SPP coupling regime scales as Noλc2/1/do. The threshold quickly rises as do approaches zero. Accordingly, for 0 < do < 0.048, the systems get into no amplification regime. Exact determination of the pumping threshold resulting in the gain effect is complicated by the finite propagation time of the coupled Maxwell–Bloch equations. Even in the presence of a gain, there might not be enough time to develop detectable increase in the signal to see the amplification effect. Due to the complexity of the problem and thus long execution times, the total propagation time is limited. Thus, for the analysis below, we adopt the smallest inversion parameter, do = 0.048, allowing for an observable gain associated with 500 fs of the propagation time.

Due to the in-plane symmetry of the array, the second harmonic is generated primarily due to the breaking up-down symmetry in the z-direction, i.e., the pre-defined directionality of the k-vector of the pump results in SHG. It should be noted that the efficiency of SHG can, of course, be greatly enhanced by properly changing the shape of the nanoparticles (e.g., replacing nanopillars with L-shaped particles). For the adopted geometry with an in-plane mirror symmetry and incident field polarized along the x-axis, the diagonal element of the second order susceptibility tensor is χxxx(2)=0 and non-vanishing components {χxyy(2),χxxy(2)=χxyx(2)} make the y-component of lattice polarization the primary contributor to the SHG. Figure 5 examines three scenarios: (blue) the array without emitters, (red) the array with emitters and no incoherent pumping (do = −1), and (magenta) the array of incoherently pumped QEs above inversion (do = 0.048). For all these cases, the incident field is driving the system at QE transition frequency. Here, we calculate separately power spectra polarized parallel (x-axis) to the incident driving field, Fig. 5(a), and perpendicular (y-axis) to it, Fig. 5(b). In the plot, the second and third harmonics are clearly observed for both polarizations.

FIG. 5.

Nonlinear power spectra calculated for the plasmonic lattice without QE (blue), with QEs and no incoherent pump (do = −1) (red), and inverted QEs with do = 0.048 (magenta). Panel (a) shows a signal polarized along the incident field polarization (x-axis) and panel (b) shows a signal polarized perpendicular to the incident field (y-axis). Curves in both panels are normalized to the total signal value at the fundamental frequency.

FIG. 5.

Nonlinear power spectra calculated for the plasmonic lattice without QE (blue), with QEs and no incoherent pump (do = −1) (red), and inverted QEs with do = 0.048 (magenta). Panel (a) shows a signal polarized along the incident field polarization (x-axis) and panel (b) shows a signal polarized perpendicular to the incident field (y-axis). Curves in both panels are normalized to the total signal value at the fundamental frequency.

Close modal

We note that the y-polarized signal exhibits significantly stronger SHG (relative to the signal at the fundamental frequency) due to the symmetry considerations above. Once QEs are added, the SHG seen in the y-polarized signal is slightly reduced. This is primarily due to the fact that QEs do not generate second harmonic directly and act as absorbers at the fundamental frequency. However, when conditions for the polariton gain are met, the picture drastically changes. The y-polarized signal at the second harmonic exceeds that at the fundamental frequency as a result of strong local field enhancement. We note that this effect is predominantly observed for the y-polarized signal. This is again due to the in-plane symmetry and the fact that emission by the QEs is isotropic, thus leading to the SHG enhancement only in the direction perpendicular to the incident field polarization. It is important to emphasize that the case with gain is not trivial as emitters do not directly generate SHG nor do they have pre-set coherences that would lead to the emission. The coherences are induced by the local interaction with the corresponding plasmon mode of the lattice, which in turn enhances stimulated emission and this in turn leads to further local field enhancement significantly altering the overall signal.

Choosing QEs resonant to the corresponding surface lattice plasmon mode may seem like an obvious choice to enhance SHG. However, due to strong coupling between QEs and SPP and because of high complexity of the parameter landscape that could significantly alter SHG, we scan the driving field frequency and explore the SHG efficiency. Similar to Sec. II B, the latter is defined as the ratio of radiated powers at the second harmonic and the fundamental frequency. When calculating the SHG efficiency, one needs to carefully take into account the material dispersion as transmission and reflection significantly vary with the driving field frequency.

Figure 6 presents main results of this section. First, we consider no incoherent pumping of QEs, do = −1 (black and green lines). The SHG efficiency varies substantially over a wide range exhibiting strong enhancement near 2 eV and 2.12 eV (reaching the value of 6 × 10−8), which corresponds to the localized plasmon resonance [see Fig. 4(b)]. This is not surprising as the localized plasmon is characterized by a strong local field enhancement (which in turn leads to high damping since the mode is partially localized inside the metal). When QEs are resonant with the lattice plasmon mode at 1.96 eV, the SHG efficiency (blue line) displays similar complex behavior and peaks at the transition frequency. Additionally and most surprisingly, we observe another broad peak at 2.12 eV that is noticeably different from the QE resonance and at which the SHG efficiency is ten times higher compared to the value at the driving field frequency. When QE transition frequency changed to 2.12 eV (red line), the SHG efficiency is slightly increased at the resonance compared to the off-resonant case (blue line). As expected, the peak at 1.96 eV is no longer seen but the enhancement of the SHG efficiency at 2.12 eV is not as strong as it is observed for the 1.96 eV case. The latter raises the efficiency by two orders of magnitude compared to the case with no incoherent pump. When QEs are strongly coupled to the localized surface plasmon mode, our simulations demonstrate that the SHG efficiency is significantly increased. This is attributed to the local field enhancement near the surface of a metal known to be higher for localized plasmon modes compared to the Bragg plasmons and lattice modes.

FIG. 6.

SHG efficiency as a function of the incident photon energy. Black and green lines show data with QEs with no incoherent pumping (do = −1) and transition frequencies of 2.12 eV (black) and 1.96 eV (green), respectively. The blue line is for the QEs with do = 0.048 and the transition frequency of 1.96 eV. The red line shows results for the QEs with do = 0.048 and the transition frequency of 2.12 eV. Two vertical dashed lines indicate corresponding QEs’ transition frequencies.

FIG. 6.

SHG efficiency as a function of the incident photon energy. Black and green lines show data with QEs with no incoherent pumping (do = −1) and transition frequencies of 2.12 eV (black) and 1.96 eV (green), respectively. The blue line is for the QEs with do = 0.048 and the transition frequency of 1.96 eV. The red line shows results for the QEs with do = 0.048 and the transition frequency of 2.12 eV. Two vertical dashed lines indicate corresponding QEs’ transition frequencies.

Close modal

Comparing absolute values of SHG efficiency calculated using the anharmonic Tavis–Cummings model [Figs. 2(c), 2(f), and 3] with full scale three-dimensional simulations (Fig. 6), it is important to emphasize several major differences between these models. First and foremost, the numerical simulations take into account multiple channels that lead to losses, thus noticeably diminishing the SHG efficiency. Such channels include radiation to the far-field and spatially dependent local field distribution, which both significantly lower the coupling between QEs and the corresponding resonant plasmon mode. Second, the anharmonicity parameter used in the model is most likely overestimated compared to the effective second order susceptibility, which directly comes from our hydrodynamic approach. The latter depends on the geometry used, and as we pointed out, the symmetry of the square lattice along with a shape of each nanoparticle will affect the efficiency of SHG as well. Having non-centrosymmetric nanoparticles (L-shaped particles, nanocrescents, etc.) and arranging them in hexagonal lattices will further increase the anharmonicity and thus the SHG efficiency.

We have considered plasmonic nanostructures interacting with two-level QEs under incoherent energy pump and demonstrated that in the strong coupling regime, such a system exhibits noticeable enhancement of the SHG. The driven-dissipative Tavis–Cummings model accounting for the anharmonic SPP provides us with an insight into the mechanisms behind polariton gain and their interplay with the SPP anharmonicity. As demonstrated, this results in strong SHG efficiency enhancement. However, this minimal model provides an upper boundary estimate for the SHG efficiency.

To examine more realistic scenarios, we performed three-dimensional numerical simulations of the SHG by a square lattice of Ag nanopillars coupled to QEs. The numerical model employed is based on the fully vectorial nonperturbative nonlinear hydrodynamic model for conduction electrons coupled to Maxwell–Bloch equations for QEs. The simulations support the idea of the gain enhanced SHG showing orders of magnitude increase when inverted QEs are tuned in resonance with the lattice plasmon mode and brought to the strong coupling regime. By varying the driving field frequency and tuning QEs to a localized plasmon mode, we demonstrate further enhancement of the SHG efficiency facilitated by strong local electric fields. The incident light polarization dependence of the second harmonic generation is examined and related to the symmetries of participating plasmon modes.

Finally, we point out that our local anharmonic Tavis–Cummings model provides rather a qualitative insight into the results of the simulations based on the semiclassical non-perturbative and non-local Maxwell–Bloch hydrodynamic model. Establishing a quantitative relationship between these models would require extraction of properly averaged χ(1)(ω) and χ(2)(ω) from the semiclassical calculations and their subsequent fitting with Eqs. (27) and (33). This is a subject of a separate study focused on quantization of anharmonic plasmonic cavities extending beyond the scope of this report.

All authors contributed equally to this work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

This work was performed, in part, at the Center for Integrated Nanotechnologies, an Office of Science User Facility operated for the U.S. Department of Energy (DOE) Office of Science by Los Alamos National Laboratory (Contract No. 89233218CNA000001) and Sandia National Laboratories (Contract No. DE-NA-0003525). M.S. is thankful to the Air Force Office of Scientific Research under Grant No. FA9550-19-1-0009 for the financial support.

1.
M. I.
Stockman
, “
Nanoplasmonics: Past, present, and glimpse into future
,”
Opt. Express
19
,
22029
(
2011
).
2.
J.
Weiner
, “
The physics of light transmission through subwavelength apertures and aperture arrays
,”
Rep. Prog. Phys.
72
,
064401
(
2009
).
3.
F. J.
Garcia-Vidal
,
L.
Martin-Moreno
,
T. W.
Ebbesen
, and
L.
Kuipers
, “
Light passing through subwavelength apertures
,”
Rev. Mod. Phys.
82
,
729
(
2010
).
4.
M.
Sukharev
and
T.
Seideman
, “
Phase and polarization control as a route to plasmonic nanodevices
,”
Nano Lett.
6
,
715
(
2006
).
5.
T. W.
Ebbesen
,
C.
Genet
, and
S. I.
Bozhevolnyi
, “
Surface-plasmon circuitry
,”
Phys. Today
61
(
5
),
44
(
2008
).
6.
M.
Kauranen
and
A. V.
Zayats
, “
Nonlinear plasmonics
,”
Nat. Photonics
6
,
737
(
2012
).
7.
G.
Wurtz
,
R.
Pollard
, and
A. V.
Zayats
, “
Optical bistability in nonlinear surface-plasmon polaritonic crystals
,”
Phys. Rev. Lett.
97
,
057402
(
2006
).
8.
N. C.
Panoiu
,
W. E. I.
Sha
,
D. Y.
Lei
, and
G.-C.
Li
, “
Nonlinear optics in plasmonic nanostructures
,”
J. Opt.
20
,
083001
(
2018
).
9.
S.
Sideris
and
T.
Ellenbogen
, “
Terahertz generation in parallel plate waveguides activated by nonlinear metasurfaces
,”
Opt. Lett.
44
,
3590
(
2019
).
10.
Y.
Blechman
,
E.
Almeida
,
B.
Sain
, and
Y.
Prior
, “
Optimizing the nonlinear optical response of plasmonic metasurfaces
,”
Nano Lett.
19
,
261
(
2018
).
11.
H.
Maekawa
,
E.
Drobnyh
,
C. A.
Lancaster
,
N.
Large
,
G. C.
Schatz
,
J. S.
Shumaker-Parry
,
M.
Sukharev
, and
N.-H.
Ge
, “
Wavelength and polarization dependence of second-harmonic responses from gold nanocrescent arrays
,”
J. Phys. Chem. C
124
,
20424
(
2020
).
12.
P.
Törmä
and
W. L.
Barnes
, “
Strong coupling between surface plasmon polaritons and emitters: A review
,”
Rep. Prog. Phys.
78
,
013901
(
2014
).
13.
M.
Sukharev
and
A.
Nitzan
, “
Optics of exciton-plasmon nanomaterials
,”
J. Phys.: Condens. Matter
29
,
443003
(
2017
).
14.
E.
Drobnyh
and
M.
Sukharev
, “
Plasmon enhanced second harmonic generation by periodic arrays of triangular nanoholes coupled to quantum emitters
,”
J. Chem. Phys.
152
,
094706
(
2020
).
15.
C.
Li
,
X.
Lu
,
A.
Srivastava
,
S. D.
Storm
,
R.
Gelfand
,
M.
Pelton
,
M.
Sukharev
, and
H.
Harutyunyan
, “
Second harmonic generation from a single plasmonic nanorod strongly coupled to a WSe2 monolayer
,”
Nano Lett.
(published online
2020
).
16.
G.-M.
Pan
,
D.-J.
Yang
,
L.
Zhou
,
Z.-H.
Hao
, and
Q.-Q.
Wang
, “
Enhanced second harmonic generation by mode matching in gain-assisted double-plasmonic resonance nanostructure
,”
Sci. Rep.
7
,
9776
(
2017
).
17.
M. S.
Tame
,
K. R.
McEnery
,
Ş. K.
Özdemir
,
J.
Lee
,
S. A.
Maier
, and
M. S.
Kim
, “
Quantum plasmonics
,”
Nat. Phys.
9
,
329
(
2013
).
18.
M.
Ramezani
,
Q.
Le-Van
,
A.
Halpin
, and
J. G.
Rivas
, “
Nonlinear emission of molecular ensembles strongly coupled to plasmonic lattices with structural imperfections
,”
Phys. Rev. Lett.
121
,
243904
(
2018
).
19.
M.
Richter
,
M.
Gegg
,
T. S.
Theuerholz
, and
A.
Knorr
, “
Numerically exact solution of the many emitter–cavity laser problem: Application to the fully quantized spaser emission
,”
Phys. Rev. B
91
,
035306
(
2015
).
20.
Y.
Zhang
and
V.
May
, “
Theory of molecule metal nano-particle interaction: Quantum description of plasmonic lasing
,”
J. Chem. Phys.
142
,
224702
(
2015
).
21.
T. B.
Hoang
,
G. M.
Akselrod
,
A.
Yang
,
T. W.
Odom
, and
M. H.
Mikkelsen
, “
Millimeter-scale spatial coherence from a plasmon laser
,”
Nano Lett.
17
,
6690
(
2017
).
22.
M.
Ramezani
,
A.
Halpin
,
A. I.
Fernández-Domínguez
,
J.
Feist
,
S. R.-K.
Rodriguez
,
F. J.
Garcia-Vidal
, and
J. G.
Rivas
, “
Plasmon-exciton-polariton lasing
,”
Optica
4
,
31
(
2017
).
23.
M. O.
Ramírez
,
P.
Molina
,
A.
Gómez-Tornero
,
D.
Hernández-Pinilla
,
L.
Sánchez-García
,
S.
Carretero-Palacios
, and
L. E.
Bausá
, “
Hybrid plasmonic–ferroelectric architectures for lasing and SHG processes at the nanoscale
,”
Adv. Mater.
31
,
1901428
(
2019
).
24.
J.
Guan
,
L. K.
Sagar
,
R.
Li
,
D.
Wang
,
G.
Bappi
,
W.
Wang
,
N.
Watkins
,
M. R.
Bourgeois
,
L.
Levina
,
F.
Fan
 et al, “
Quantum dot-plasmon lasing with controlled polarization patterns
,”
ACS Nano
14
,
3426
(
2020
).
25.
A.
Pusch
,
S.
Wuestner
,
J. M.
Hamm
,
K. L.
Tsakmakidis
, and
O.
Hess
, “
Coherent amplification and noise in gain-enhanced nanoplasmonic metamaterials: A Maxwell-Bloch Langevin approach
,”
ACS Nano
6
,
2420
(
2012
).
26.
J. M.
Winkler
,
M. J.
Ruckriegel
,
H.
Rojo
,
R. C.
Keitel
,
E.
De Leo
,
F. T.
Rabouw
, and
D. J.
Norris
, “
Dual-wavelength lasing in quantum-dot plasmonic lattice lasers
,”
ACS Nano
14
,
5223
(
2020
).
27.
A.
Piryatinski
,
O.
Roslyak
,
H.
Li
, and
E. R.
Bittner
, “
Nonequilibrium states of a plasmonic Dicke model with coherent and dissipative surface-plasmon–quantum-emitter interactions
,”
Phys. Rev. Res.
2
,
013141
(
2020
).
28.
J.-P.
Martikainen
,
M. O. J.
Heikkinen
, and
P.
Törmä
, “
Condensation phenomena in plasmonics
,”
Phys. Rev. A
90
,
053604
(
2014
).
29.
S.
Zaster
,
E. R.
Bittner
, and
A.
Piryatinski
, “
Quantum symmetry breaking of exciton/polaritons in a metal-nanorod plasmonic array
,”
J. Phys. Chem. A
120
,
3109
(
2016
).
30.
T. K.
Hakala
,
A. J.
Moilanen
,
A. I.
Väkeväinen
,
R.
Guo
,
J.-P.
Martikainen
,
K. S.
Daskalakis
,
H. T.
Rekola
,
A.
Julku
, and
P.
Törmä
, “
Bose–Einstein condensation in a plasmonic lattice
,”
Nat. Phys.
14
,
739
(
2018
).
31.
M.
Ramezani
,
M.
Berghuis
, and
J. G.
Rivas
, “
Strong light–matter coupling and exciton-polariton condensation in lattices of plasmonic nanoparticles
,”
J. Opt. Soc. Am. B
36
,
E88
(
2019
).
32.
J.
Butet
,
P.-F.
Brevet
, and
O. J. F.
Martin
, “
Optical second harmonic generation in plasmonic nanostructures: From fundamental principles to advanced applications
,”
ACS Nano
9
,
10545
(
2015
).
33.
P.
Kirton
,
M. M.
Roses
,
J.
Keeling
, and
E. G.
Dalla Torre
, “
Introduction to the Dicke model: From equilibrium to nonequilibrium, and vice versa
,”
Adv. Quantum Technol.
2
,
1800043
(
2019
).
34.
L.
Novotny
and
B.
Hecht
,
Principles of Nano-Optics
(
Cambridge University Press
,
New York
,
2008
).
35.

Here, numerical prefactor of 16 reflects scaling of the radiated power with emission frequency as P(ω) ∼ ω4.

36.
M.
Sukharev
and
A.
Nitzan
, “
Numerical studies of the interaction of an atomic sample with the electromagnetic field in two dimensions
,”
Phys. Rev. A
84
,
043802
(
2011
).