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.

## I. INTRODUCTION

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 propagating^{2,3} and localized surface plasmon-polaritons (SPPs) and utilizing their unique polarization features to control light propagation^{4} at various metal–dielectric interfaces.^{5} The nonlinear optics of plasmonic nanosctructures has quickly become a hot research field.^{6} Applications range from optical bistability^{7} 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 response^{32} 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.

## II. DRIVEN-DISSIPATIVE ANHARMONIC TAVIS–CUMMINGS MODEL

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 $\psi ^$ with Hermitian conjugate $\psi ^\u2020$. The SPP interacts with an ensemble of $No$ identical two-level QEs described by spin operators $\u015dn\xb1=\u015dnx\xb1i\u015dny$ and $snz$, expressed in terms of the Pauli SU(2) operators as $\u015dnj=12\sigma ^nj$ with *j* = *x*, *y*, *z* and site index $n=1,No\xaf$. The QEs are subject to incoherent energy pump, whereas the SPP interacts coherently with an incident monochromatic laser field *E*_{in}(*ω*). 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(\omega )\u223c\mu sp\psi ^(\omega )$ and $Eout(2\omega )\u223c\mu sp\psi ^(2\omega )$ oscillating with fundamental and second-harmonic frequencies, respectively. Both QEs and MNPs interact with the environment, causing dissipative processes.

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

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, *E*_{in}, is accounted for in the last term of the Hamiltonian where the coupling strength is represented by the Rabi frequency Ω_{in} = (*μ*_{sp} · *E*_{in})/*ℏ*.

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

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

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.

### A. Mean-field steady states

Starting with Eqs. (2) and (3) and the Hamiltonian (1), we derive the Heisenberg operator equations for $O^={\psi ^,\u015dn\u2212,\u015dnz}$, and adopting the mean-field approximation replaces the operators with their averages. The result is

where we introduce normalized SPP, $\phi =\u27e8\psi ^\u27e9/No$ and QE, $s\xb1=\u2211n=1No\u27e8\u015dn\xb1\u27e9/No$, coherences, and normalized QE population $sz=\u2211n=1No\u27e8\u015dnz\u27e9/No$ varying in the interval [−1/2, 1/2]. Note that the effective QE-SPP coupling, $No\lambda $, and anharmonicity, $No\alpha $, 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

referred to as the population inversion parameter. Depending on the incoherent pumping rate, *γ*_{↑}, it can be −1 < *d*_{o} < 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 < *d*_{o} < 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

where slowly varying amplitudes for the fundamental SPP, $\phi \u03031(t)$, and QE, $s\u0303\u2212(t)$, coherences as well as for the second-harmonic SPP coherence $\phi \u03032(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,

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, $\phi \u03031$. 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:

determining the steady states of $\phi \u03031$. Here, the prefactor

contains the real and imaginary part functions of $|\phi \u03031|$. The imaginary part,

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,

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,

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

Equation (14) has a trivial solution $\phi \u03031=0$ and a nontrivial solution defined by $\Sigma |\phi \u03031|=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.,

These conditions determine the lasing phase transition in the driven-dissipative Tavis–Cummings model. Roots of Eqs. (19) and (20) are coherences $|\phi \u03031|$ spontaneously appearing above critical coupling, *λ*_{c}, and the detuning parameter *δω* describing the lasing frequency shift.^{19,27,33} Below *λ*_{c}, the trivial solution, $\phi \u03031=\delta \omega =0$, holds.

### B. Polariton dispersion and SHG efficiency

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,

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 $\phi \xaf1$. Making the substitution of $\phi \u03031(t)=\phi \xaf1+\delta \phi 1(t)$ into Eq. (21) and complementing the resulting equation with a complex conjugate, we derive linearized equations of motion in the matrix form,

with the column vectors $\delta \Phi =(\delta \phi 1,\delta \phi 1*)T$ and $F(t)=e\u2212i(\omega \u2212\omega o+\delta \omega )t,\u2212ei(\omega \u2212\omega o+\delta \omega )tT$. Here, $I$ is the 2 × 2 unit matrix and

stands for the so-called stability matrix with matrix elements

depending on $\Sigma |\phi \xaf1|$ [Eqs. (15)(16)(17)(18)] and its functional derivative both functions of the steady state solutions, $\phi \xaf1$, of Eq. (14).

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

where the linear susceptibility is

and depends on the stability matrix, $M$, eigenvalues

and the upper components of orthonormal left $\delta \Phi \xb1L=(u\xb1L,v\xb1L)$ and right $\delta \Phi \xb1R=(u\xb1R,v\xb1R)T$ eigenvectors. Poles of the linear susceptibility,

determine the polariton frequency, $\omega \xb1=Re\zeta \xb1$, and dephasing rate $\gamma \xb1=Im\zeta \xb1$.

Below the critical coupling *λ* < *λ*_{c}, $\phi \u03031=\delta \omega =0$, simplifying Eq. (28) to the form

where the second term in the right-hand side is an explicit representation of the polariton gain term $G{|\phi \xaf1|=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,

where *γ*_{o} depends on the incoherent pumping rate and *d*_{o}. 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\lambda 2>\gamma sp\gamma 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 $Mo2\u2212ImMd2$. 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)], $\delta \phi 12$, as the driving force. Accordingly, the solution of Eq. (11) gives the following SPP second-harmonic polarization:

in terms of the second-order nonlinear susceptibility,

which in turn depends on *χ*^{(1)}(*ω*) [Eq. (27)].

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

explicitly indicating quadratic scaling of the SHG efficiency both with *χ*^{(1)} and the anharmonic coupling $No\alpha $.

### C. Numerical analysis

We perform the analysis by calculating the polariton dispersion [Eqs. (28) and (29)] and related absorption/emission spectra determined by $Im\chi (1)(\omega )$ [Eq. (27)] and the SHG efficiency [Eq. (34)] depending on SPP-QE coupling strength $No\lambda 2/\omega 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 *d*_{o} = 0.05. For these parameters, the critical coupling value is $No\lambda c2/\omega o2=1.9\xd710\u22123$. These parameters are similar to those adopted in the simulations discussed in Sec. III. Finally, we introduce two limiting regimes of weak, $No\alpha /\omega o=10\u22123$, and strong, $No\alpha /\omega 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.

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\alpha 2$ [see Eq. (34)].

Finally, we fix the coupling strength slightly above the critical coupling by setting $No\lambda 2/\omega o2=2.0\u2009\xd7\u200910\u22123$ and vary the incoherent pumping rate below the inversion threshold to sample the states characterized by *d*_{o} = −0.33 (*γ*_{↑} = 0.5*γ*_{↓}) and *d*_{o} = 0.05 (*γ*_{↑} = 0.9*γ*_{↓}). Comparison of the SHG efficiencies calculated below inversion with that for the inverted QEs (*d*_{o} = 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.

## III. SHG FACILITATED BY TWO-DIMENSIONAL LATTICE OF METAL NANOPILLARS

### A. Computational model

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}

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

where *P*_{o} = *ρ*_{o}*μ*_{o}*s*_{x} denotes the QE polarization with *s*_{x} = 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}

Here, the plasma frequency is $\omega p=n0e2/\epsilon 0me*$, the electron gas damping parameter is set to *γ*_{e} = 0.14 eV, and the number density of conduction electrons is *n*_{0} = 5.86 × 10^{28} 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 × 10^{7}. 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 10^{26} 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 *E*_{in} = 2 × 10^{−2} V/nm.

### B. Simulation results

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 − *T* − *R*. 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 *d*_{o} = 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 < *d*_{o} < 0.048 destroys the gain. For example, a weak inversion with *d*_{o} = 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 *d*_{o} ≪ 1 the critical parameter determining the strong QE-SPP coupling regime scales as $No\lambda c2/\u223c1/do$. The threshold quickly rises as *d*_{o} approaches zero. Accordingly, for 0 < *d*_{o} < 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, *d*_{o} = 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 $\chi xxx(2)=0$ and non-vanishing components ${\chi xyy(2),\chi xxy(2)=\chi 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 (*d*_{o} = −1), and (magenta) the array of incoherently pumped QEs above inversion (*d*_{o} = 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.

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, *d*_{o} = −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.

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.

## IV. CONCLUSIONS

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.

## AUTHORS’ CONTRIBUTIONS

All authors contributed equally to this work.

## DATA AVAILABILITY

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

## ACKNOWLEDGMENTS

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.

## REFERENCES

_{2}monolayer

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