Near-field enhancement of optical second harmonic generation in hybrid gold–lithium niobate nanostructures

Nanophotonics research has focused recently on the ability of nonlinear optical processes to mediate and transform optical signals in a myriad of novel devices, including optical modulators, transducers, color filters, photodetectors, photon sources, and ultrafast optical switches. The inherent weakness of optical nonlinearities at smaller scales has, however, hindered the realization of efficient miniaturized devices, and strategies for enhancing both device efficiencies and synthesis throughput via nanoengineering remain limited. Here, we demonstrate a novel mechanism by which second harmonic generation, a prototypical nonlinear optical phenomenon, from individual lithium niobate particles can be significantly enhanced through nonradiative coupling to the localized surface plasmon resonances of embedded gold nanoparticles. A joint experimental and theoretical investigation of single mesoporous lithium niobate particles coated with a dispersed layer of ~10 nm diameter gold nanoparticles shows that a ~32-fold enhancement of second harmonic generation can be achieved without introducing finely tailored radiative nanoantennas to mediate photon transfer to or from the nonlinear material. This work highlights the limitations of current strategies for enhancing nonlinear optical phenomena and proposes a route through which a new class of subwavelength nonlinear optical platforms can be designed to maximize nonlinear efficiencies through near-field energy exchange.


I. INTRODUCTION
Ultrafast optical frequency conversion using nonlinear optical (NLO) harmonic generation in micro-and nanoscale materials is emerging as an important phenomenon in the basic and applied study of photonics, 1 as well as for the development of novel sensing and imaging techniques in materials science, chemistry, and biology. 2 Perhaps the simplest of the ultrafast NLO processes relevant to each of these fields is second harmonic generation (SHG), in which two input fields oscillating at a fundamental frequency are coherently combined into an output field oscillating at twice the fundamental. 3,4 In general, SHG is prized for its narrowband operation, relatively high efficiency among NLO processes, and natural occurrence in a wide array of inorganic materials. 3,4 However, even in "good" NLO materials, the nonlinear wave motion underpinning SHG is generally either weak 5 or slower than the ∼1-ps timescales desired in ultrafast applications. 6 To overcome this difficulty, SHG photonics designs have until recently centered around the use of large NLO crystals within which the fundamental wavelength (FW) and SH light waves can be constructively interfered to overcome intrinsically low SHG efficiencies. 7,8 It has become clear, however, that much smaller micro-and nanoscopic SHG devices are more desirable for cutting-edge applications like fewmolecule sensing, 9 bioimaging, 10 and phase-resolved light emission 11 due to the unique degree of control over the spectrum and spatial distribution of the scattered fields that miniaturized optical systems provide. 12 Unfortunately, due to the complicated field profiles of nanostructured particles, interference and phasematching techniques are rendered ineffective. As a result, the absolute efficiency of an ultrafast NLO process at the nanoscale tends to depend on the microscopic volume of NLO material used. 13 Researchers continue to pursue fundamentally new NLO enhancement techniques and novel nanofabrication solutions to optimize or supersede this limit, exploring the NLO properties of light confined within resonant dielectric, 14 metallic, 15 and hybrid 16 micro-and nanostructures, as well as within nanopatterned waveguides. 17 Our investigation addresses both the scientific and engineering challenges posed in NLO nanophotonics by revealing a new and practical strategy for enhancing SHG radiation. Specifically, we use the surface-localized nearfields of disordered arrays of Au nanoparticles (NPs) to create hitherto unexplored energy-transfer pathways for nanolocalized upconverted light. These disordered nanostructures, comprised of ∼1-µm diameter mesoporous lithium niobate (LiNbO 3 ) microspheres coated in a dispersed layer of 10-nm diameter Au NPs, allow for a straightforward and high-throughput synthetic process devoid of the difficulties of precision nanofabrication while simultaneously producing an NLO system with a large number of regions of enhanced field intensity (socalled "hot spots") between the LiNbO 3 and Au surfaces.
In stark contrast, many past investigations [18][19][20][21][22][23] have achieved enhancement factors between 5 and 20 by employing carefully-tailored hybrid nanostructures. The manufacture of such structures is technically challenging, causing difficulties with the reproduction of large enhancement ratios between samples. 22 Other studies have demonstrated enhancement factors of ∼100-1000 more reliably using simple core-shell SHG systems. 16,24 How-ever, the field profiles and resonance structures of shells, which are often assembled as collections of large NPs, can be complex and difficult to separate from the influence of imperfections, 25 and can be susceptible to degradation under exposure to high-intensity lasers. 26 As such, we have chosen to use a more rarefied ensemble of smaller and more easily modeled NPs to isolate the NLO enhancing abilities and governing physical principles of an optically robust and reproducible proof-of-concept nanostructure design.
As will be shown below, our multiple-hot-spot geometry enhances SHG by a factor of 32. Further, this enhancement is clearly attributed through a complete make-measure-model synthesis and characterization approach to the near-field localization and resonant behaviors of metallic nanoantennas that are generally considered unsuitable for the task of enhancing far-field (i.e., radiative) emission. Although a few very recent studies 6,27 have measured moderate (e.g., ∼7×) SHG enhancements using similarly disordered NP ensembles, each was either unable to isolate single nanostructures (Ref. 27) or to thoroughly explore the governing mechanism behind the observed SHG enhancements (Ref. 6). Our results, therefore, represent both the first demonstration of NLO enhancement values of >10 from single disordered nanostructures, as well as the discovery of as-yet unknown mechanisms by which strongly-subwavelength NPs can enhance SHG.

A. Preparation of hybrid SHG microspheres
Selecting an NLO material with a large NLO susceptibility coefficient is critical to improve the efficiency of the SHG conversion process at smaller scales. Relative to many NLO materials, LiNbO 3 is a unique photonic material, often referred to as the "silicon of photonics" due to its relatively large second-order susceptibilities (e.g., 41.7 pm/V), higher optical damage resistance (e.g., 20 W/cm 2 for congruent LiNbO 3 crystals), and a relatively wide window of optical transparency (e.g., 400 to 5 000 nm). 3,28 In the first step, monodisperse mesoporous LiNbO 3 particles were prepared by modifying a previously reported method (details in Materials and Methods Section S.1). 4 Scanning transmission electron microscopy (STEM) images confirm that diameters of the LiNbO 3 particles were ∼1 000 nm (Figure 1b-d and Figure S1) and that each particle has a mesoporous structure. These mesoporous particles contained a randomly disordered network of distinct grains in its bulk and a textured surface that contains randomly distributed ∼20nm diameter pores. In addition, X-ray diffraction (XRD) patterns of the LiNbO 3 particles in a powder form were acquired and correlate well with a reported LiNbO 3 reference ( Figure S2), displaying a distinct rhombohedral phase (space group R3c) that accompanies strong SHG behavior. 3,28 Further evidence of the rhombohedral phase of the LiNbO 3 particles was provided by Raman spectroscopy of a powdered sample ( Figure S3), in which the observed bands are characteristic of rhombohedral LiNbO 3 and agree well with a commercially available LiNbO 3 standard. 28 In the second step to prepare the hybrid Au-LiNbO 3 particles, we utilized an in situ seed-mediated growth to load the surfaces of the mesoporous LiNbO 3 with Au NPs (Figure 1a). This approach used a one-pot synthetic technique that eliminated the need for precision nanoengineering techniques to position the plasmonic materials upon their NLO support. 29 Briefly, an aqueous suspension of bare, mesoporous LiNbO 3 particles in the presence of Au chloride as a precursor was heated at 70 • C for 3 h. The in situ reduction of the Au precursor using sodium citrate as a reducing agent led to the formation of the hybrid structures, wherein the Au NPs were observed both through electron microscopy imaging (Figure 1b . These Au nanoparticles were deposited on the outer surfaces, as well as within the pores of the mesoporous LiNbO 3 . Synthesis of hybrid particles prepared from porous LiNbO 3 supports revealed that the number of Au NPs loaded onto each LiNbO 3 particle was ∼900-1 000.

B. Linear and nonlinear optical spectroscopy
Analysis of the SHG emission from the hybrid nanostructures was performed using a Leica SP5 commercial two-photon microscope. The incident light was generated using a pulsed, mode-locked Ti:sapphire laser with a pulse width of ∼140 fs. The resonances of the NLO particles have periods between 0.1 and 0.3 fs such that their driven oscillatory behavior and associated SHG scattering are independent of the pulse envelope. Additionally, with a repetition rate of 80 MHz, the spacing between pulses (12.5 ns) is much longer than the lifetimes of the target resonances (2-30 fs) such that scattering data from each pulse is considered to be uncorrelated to the data obtained from prior pulses. Moreover, the use of short, well-separated pulses also enabled the laser to deliver the high-intensity field required to generate the observable SHG signal. The scattered light and SHG emission from single particles were each analyzed by collecting the light with a microscope objective lens and analyzing the signal using a spectrometer that enabled separation of the input signal (e.g., fundamental wavelength) from the output signal (e.g., SH wavelength). Finally, The laser frequency was tuned from 680 nm to 1 080 nm to analyze the wide spectral response of these NLO materials. Additional details of the experimental setup are shown in NPs with LiNbO3 as achieved using an in situ synthesis method. Assemblies of Au-LiNbO3 hybrid particles as characterized by scanning transmission electron microscopy (STEM) operating in: (b, c, d) a high-angle annular dark-field (HAADF) mode. Energy dispersive X-ray spectroscopy (EDS) analysis of the Au-LiNbO3 hybrids obtained by STEM techniques was used to create EDS maps of: (e) Au NPs; (f ) Au NPs overlaid on the HAADF image of the assemblies; and (g) overlaid Au-Nb signals within these assemblies. First, we acquired and compared the extinction spectra of the hybrid Au-LiNbO 3 particles, pristine LiNbO 3 particles, and Au NPs. Aqueous suspensions of each sample were prepared and characterized using a UV-visible spectrometer. No extinction peak was observed for the pristine LiNbO 3 particles indicating their optical transparency in the region from 400 to 800 nm (Figure 3a). We prepared ∼10-nm diameter Au NPs and acquired and characterized their extinction spectra to obtain a detailed comparison to the hybrid Au-LiNbO 3 particles. This suspension of Au NPs had an extinction peak at ∼520 nm due to their plasmonic response ( Figure S6). The hy-brid particles of Au-LiNbO 3 , however, exhibited a broad extinction band from the visible to the near-IR with a maximum centered at 530 nm ( Figure 3a). The presence of this peak in the hybrid Au-LiNbO 3 particles was due to the contributions of the Au NPs. The broadness and red shift in the resonance peak of the Au-LiNbO 3 particles could be attributed to the formation of Au NP aggregates on the surfaces of the LiNbO 3 particles. 30 For the linear scattering and nonlinear optical analyses of individual particles, we prepared separate, dilute aqueous suspensions (0.1 mg/mL) of both the LiNbO 3 and Au-LiNbO 3 particles via sonication. The resultant suspensions were drop cast onto glass coverslips and dried under a vacuum to evaporate the solvent. Microscopy analyses of the resulting substrates indicated the presence of well-dispersed, individual particles. The linear optical spectra of individual Au-LiNbO 3 and pristine LiNbO 3 particles were obtained to evaluate their scattering profiles (Figure 3b). A dark-field optical spectroscopy setup was used for these scattering measurements where individual particles were imaged using a 50× dark-field objective. The individual particles were illuminated by white light generated from a halogen lamp. The scattered light was collected by the same objective lens and detected with an imaging spectrometer through a pinhole that restricted the collection of signal from the area around the measured nanoparticle. We determined the position and the number of the scattering peaks in the linear optical spectra of each type of particle. Four distinct scattering peaks were observed for both the Au-LiNbO 3 hybrids and pristine LiNbO 3 particles. As the diameters of the particles were comparable to the wavelength of the incident light, the linear scattering response of these materials does not resemble a Rayleigh scattering profile. Due to a negligible change in the size of Au-LiNbO 3 hybrids relative to the pristine LiNbO 3 , the electromagnetic field inside each type of particle was consistent, resulting in indistinguishable Mie scattered modes between pristine and hybrid particles. 31 These scattering results are comparable to the Mie scattering spectra previously reported for LiNbO 3 based materials having comparable sizes. 32 The presence of distinct resonances in the scattering spectra indicate that both the Au-LiNbO 3 and LiNbO 3 particles act as optical Mie resonators. By considering these optical measurements of individual Au-LiNbO 3 and LiNbO 3 particles, it can be concluded that the Au NPs play no significant role in the Mie resonances in these materials. The features of Mie scattering were not measurable below 400 nm and above 700 nm due to the loss of sensitivity of the measurement setup in these regions. The linear optical spectra of individual Au NPs with diameters of ∼10 nm could not be detected due to their smaller size and limitations of the objective lens to locate individual Au particles at the nanoscale. A further improvement of the dark-field spectroscopy setup in the UV range is challenging because of the lack of objective lenses that have both a high transmission in the UV range and no chromatic aberrations across both the UV and visible range.
We acquired input wavelength-dependent ( Figure 4a) and power-dependent ( Figure 4b) NLO spectra for the prepared mesoporous Au-LiNbO 3 hybrids to ensure the materials are SHG active and to verify the second-order nature of the detected signals. The power dependent SHG signals were collected by varying the power of the laser at the FW of 800 nm while imaging individual Au-LiNbO 3 particles on the substrate. The image emission matrix was averaged to yield a mean value in regions of interest selected to exclude the majority of the dark pixels. The selected regions of interest were maintained unchanged throughout the series of acquired images. The logarithms of the mean values for the SHG signals were plotted against the logarithm of the laser power. The slopes of the linear fits reflect the power dependence of the SHG on the pumping intensity. A slope of ∼2.0 for the Au-LiNbO 3 particles confirm that the signal is generated from a second-order NLO process (i.e., SHG).
The SHG for individual Au-LiNbO 3 particles was also assessed for a series of discrete fundamental wavelengths. The resultant SH responses were each normalized by dividing their intensities by the maximum scattering intensity within each SH band when assessing the frequency doubling behaviour of the SHG process. A tunable SH response was observed at 400, 420, 440, 460 and 480 nm when the product was excited with FWs of 800, 840, 880, 920, and 960 nm, respectively (Figure 4a). These results correlated well to the anticipated frequency doubling of the fundamentals.
The SHG spectra were acquired for a cluster of Au NPs using the same experimental setup to verify whether the strong SHG emission originated from the surface of the gold rather than from the NLO LiNbO 3 core. We used a cluster of solid Au NPs since individual 10-nm diameter gold NPs were difficult to drop-cast and locate using the 63× objective lens. We scanned the FWs over the range from 850 nm to 1 060 nm and collected the SHG response for the clusters of Au NPs ( Figure S7). No SHG response was generated from 400 to 530 nm for the Au NPs. The SHG emitted by the Au-LiNbO 3 hybrids is, therefore, exclusively generated by the LiNbO 3 particles. For NPs made out of centrosymmetric materials such as gold, SHG is mainly emitted from an Angstromthin layer near the particle surface since SHG emission is mainly induced by the surface normal component of the second-order nonlinear susceptibility tensor. 33,34 The Au NPs used in this work are, however, spherical and the SHG generated from the bulk of the gold nanoparticle and through the surface-tangential component of the second-order nonlinear susceptibility is expected to be very low and are ignored in this study ( Figure S7). 34 We measured the SHG spectra of individual pristine LiNbO 3 and hybrid Au-LiNbO 3 particles by sweeping the incident laser over the wavelength range from 850 to 1 070 nm in steps of 5 nm to investigate whether Au NPs play a role in the SHG enhancement of the hybrid structures (Figures 4c and S14a). The corresponding SHG response was recorded from 425 to 530 nm. A comparison was made between the SHG intensity from individual Au-LiNbO 3 particles with individual pristine LiNbO 3 particles examined under the same experimental parameters to quantitatively determine the NLO enhancement. This method to calculate the SHG enhancement has been used in the reported literature, such as where the enhancement of the NLO signals from metasurfaces is calculated by comparing the patterned metasurface to the unpatterned original surface. 35,36 After normalization of the SHG output of the Au-LiNbO 3 hybrids to the SHG output of the bare LiNbO 3 , we were able to calculate the enhancement factors. The obtained enhancement factor was plotted against the FWs of the measurement. The highest enhancement value was obtained at a FW of 1 030 nm, reaching an enhancement of 32 times (Figures 4c and S14a). This evaluation demonstrates a corresponding enhancement in the SHG response at 515 nm by a factor of 32 for the mesoporous LiNbO 3 particles loaded (a) Measured extinction spectra for ensembles of bare LiNbO3 particles (black) and of hybrid Au-LiNbO3 particles (red) compared to the modeled dipole plasmon extinction (gold). These spectra show a plasmonic peak at ∼530 nm for the hybrid Au-LiNbO3 particles, while the bare LiNbO3 particles were optically transparent. (b) Light scattering spectra for individual, bare LiNbO3 particles and individual Au-LiNbO3 hybrid particles measured in a dark field (DF) configuration with white light illumination. The insets in (b) portray grey color images (scale bars = 2 µm) for (black) a single LiNbO3 particle and (red) a single Au-LiNbO3 hybrid particle supported on glass coverslips. The data were normalized to the maximum value of the scattering intensity. (c) Theoretical scattering cross sections (green) from a model microsphere with dielectric constant 6.3 + 0.05i and radius 700 nm in comparison with experimentally measured scattering data from a single, bare LiNbO3 particle (black), normalized to range from 0 to 1 within the selected frequency window.
with Au NPs. We prepared Au-LiNbO 3 hybrids with a lower loading of Au NPs to understand what role the Au NPs play in enhancing the SHG response ( Figures S8,  S9). The average number of Au NPs on each LiNbO 3 particle was ∼6. We acquired the SHG response of these Au-LiNbO 3 hybrids with their lower loading of Au NPs and compared this response to the NLO response of the bare LiNbO 3 particles. Interestingly, no enhancement in the SHG response was observed for these hybrids of Au-LiNbO 3 with only a few Au NPs ( Figure S10). These observations indicated the need to have a higher loading of Au NPs on the surfaces of the LiNbO 3 to enhance the SHG response.

C. Interpretation of scattering signals
With plasmon-induced enhancements greater than 30, the hybrid particles in this investigation show commensurate performance with other NLO nanostructures of a similar size without the requirements of being prepared using precision nanoengineering techniques. For example, an SHG emission enhancement factor of 30 is reported in Ref. 37 and a factor of ∼22 is reported in Ref. 38. However, the plasmonic Au NPs used in this study are not properly tuned to act as efficient nanoantennas, obscuring the phenomena underpinning their ability to magnify SHG.
More precisely, hybrid SHG nanostructures generally employ metal NPs to perform one of two functions relevant to upconversion. The first function, in which plasmons enhance the in-coupling of light to the SHG material, employs the near-fields of highly polarizable plasmons that are tuned to the FW to augment the pump field of an NLO process. 16,37 In the second function, radiative plasmons tuned to the up-or down-converted frequency enhance the out-coupling of light by extracting energy from the NLO material through near-field interactions and then quickly scattering it to the far-field, i.e., via the Purcell effect. 19,38 Carefully designed nanostructures can make use of both strategies, 22 and antenna effects can be significant whether or not the NLO material itself is optically resonant. 21 Our nanoantennas are far too small to act as radiative out-coupling nanoantennas yet have a plasmon resonance that is aligned with the SH output maximum. The most straightforward explanations for the surprising efficacy of these "bad antenna" NPs are that the NPs' enhanced SHG emission is a result of either of the following hypotheses: (i) coupling effects, in which optical energy is routed through the system via more complicated interactions than can be captured by a simple antenna picture; or (ii) superradiant plasmon scattering generated by their collective interaction with the radiation field.
No existing model can fully address either of these questions. Recent theoretical studies have investigated several related phenomena, including radiation from dipole-driven Mie scatterers, 39 NLO scattering from bare dielectric spheres, 40 NLO emission from quantum emitters in idealized Fabry-Pérot cavities, 41,42 and light emis-sion from ensembles of oscillating dipoles. 43,44 However, a novel synthesis of the ideas from each study is required to uncover the operative SHG enhancement mechanism(s) in our nanostructures.
Beginning with hypothesis (i), we first simplify the particle geometries highlighted in Figure 1 in a manner consistent with prior investigations of complicated multicrystalline scatterers. 13 In detail, we model each microsphere as a smooth, isotropic spherical particle of radius a 2 = 500 nm. We also simplify the form of the incoming light, assuming the laser field E 0 (r, ω) to be a monochromatic plane wave of characteristic field strength E 0 and frequency ω 0 that travels in the z-direction and is linearly polarized along x. That is to say, The combination of LiNbO 3 's nearly constant dielectric function 2 = 5.5 + 0.035i (see Ref. 45 and also the Supplementary Information [SI], Section S.2.1) and relatively weak second-order nonlinear susceptibility |χ (2) (ω 0 , ω − ω 0 )| ∼ 10 −6 cm/statV (SI, Section S.2.1 and Figure S11) for observation frequencies ω and laser frequencies in the optical region of interest allows for simplification of the LiNbO 3 's NLO behaviors. Specifically, the nonlinear wave motion of an excited microsphere can be perturbatively expanded as the superposition of Mie resonances at the laser FW, SH, third harmonic, and so on (SI, Section S.2.4).
The FW scattered electric field E sca (r, ω) is simply the component of the total electric field for which the output frequency matches the laser frequency. The FW resonances in general have both appreciable magnetic and electric fields, have nondipolar symmetry, and produce a field E sca (r, ω) = α a − −2
From here, in accordance with previous work, 40 we assume that upconversion in the LiNbO 3 spheres during the experiment is coherent and occurs when a small portion the energy bound in the laser and FW scattered fields is converted into a second material polarization oscillating at the SH (2ω 0 ). This so-called nonlinear polarization field is always proportional to |χ (2) 2 (ω 0 , ω 0 )|E 2 0 (see Eq. [S38]) and generates two SH electric fields of its own.
The first SH field, E 0 (r, ω), describes direct radiation by the nonlinear polarization. It does not interact with the LiNbO 3 sphere directly but does excite the Au NP ensemble. The second, E sca (r, ω) = , is a scattered field similar to the scattered FW electric field E sca (r, ω) but is driven by the nonlinear polarization rather than directly by the laser light. In fact, E sca (r, ω) can be defined identically to its FW counterpart but with the sum over α, i.e. the sum over the set of Mie resonances near the FW, replaced by a sum over the set of Mie resonances β = {T , p , , m } near the SH. The SH scattered fields can also interact with the NP ensemble which, in this case, acts as a set of symmetry-breaking resonant cavities that allow mixing between the otherwise orthogonal SH Mie resonances of the LiNbO 3 sphere.
Whether excited by E 0 or E sca , the excited NPs can transfer energy back to the LiNbO 3 sphere through their near-fields. Thus, the two SH field-NP interactions provide additional pathways for energy to be transferred from the abundant incident power 31 P 0 = ca 3 2 E 2 0 /8 of the laser to the LiNbO 3 microsphere, which then scatters more SH light as a result. Figure 4d provides a sketch of these processes.
Quantification of the resulting SHG enhancements is dramatically simplified by the fact that each Mie resonance β is well-approximated as a discrete oscillator with resonance frequency ω β and damping rate γ β (see SI, Eq. [S16] and Section S.2.2). The plasmons in each NP, here assumed to have radii a 1 = 5 nm, have similarly simple oscillator properties (SI, Section S.2.2), with resonance frequencies ω 1 and damping rates γ 1 . However, even with these conveniences, explicitly describing the numerous energy-transfer pathways between the N NPs of the ensemble and the LiNbO 3 is difficult. Therefore, we first model the near-field energy exchange between the LiNbO 3 and Au by considering a reduced system with only a single Au NP coupled to the LiNbO 3 modes. Letting the NP dipole be oriented along the ν Cartesian axis with a moment d ν (ω), its motion along with the motion of each SH Mie moment ρ β (ω) can be quantified from coupled equations of motion (SI, Section S.2.4), here given in the Fourier domain: In addition to the previously defined constants, each of the moments has a natural frequency Ω 1,β = ω 2 1,β + γ 2 1,β /4, an oscillator strength f ν,β , and a charge e. On the right hand side of Eq. (1), the forces F 1ν (r 0 , ω) and F 2β (ω) describe the driving of the Au NP dipoles and Mie modes, respectively, by the χ (2) 2 -upconverted incident light with the NP centered at r 0 . The analytical form of these forces is complicated and is left for the SI (see Section S.2.4). The phases of the moments' responses to these forces are encoded by their characteristic phase delays η 1 (ω) = 2ω 1 cos ψ 1 + γ 1 sin ψ 1 + 2iω sin ψ 1 and η β (ω) = (Ω * β + ω) exp(−iψ β ), respectively, wherein ψ 1,β are the characteristic phase lag parameters of the NP plasmons and Mie resonances, respectively. Numerical values for these oscillator parameters are given in Tables S1 and S2 and the particles' individual resonance profiles are plotted in Figures S12 and S13.
In addition to the upconversion forces, the moments of the system are coupled with strengths σ βν (r 0 ) = (e 2 /a 3 2 )ê ν · X β (r 0 , k β ), whereinê ν is the ν-oriented unit vector, X β (r, k β ) is the vector spherical harmonic that describes the spatial variations of the fields of the β th mode, and k β = ω β /c. In the case where the couplings strengths are set to zero, the LiNbO 3 and Au NP radiate independently and no energy passes between them. As the magnitude of σ βν (r 0 ) increases, the likelihood that a photon will be exchanged between the LiNbO 3 and Au multiple times before being radiated away also increases.
One can see from the scattering spectra of Figure 3b that no mode splitting is generated by the addition of the Au NPs to the LiNbO 3 sphere in our system. Splitting is readily observed in strongly-coupled NLO polariton systems, 46 such that weak interactions between the Au NPs and microsphere are expected. The theory is in good agreement with this observation, as the weak coupling condition 47,48 The equations of motion of the LiNbO 3 Mie modes can, thus, be solved perturbatively to capture the effects of these exchanges on the power radiated from the microsphere.
where each order in the expansion includes the effects of n exchanges of energy, the power scattered from the β th mode The power scattered by β in the absence of the NP dipoles, P Due to a convenient orthogonality condition (Eq. [S18]), the total power scattered by the microsphere at the SH frequency, P (2ω 0 ), is then simply a sum over the power scattered from each mode, such that P (2ω 0 ) = β P β (2ω 0 ), and the scattering enhancement can be simply defined as P (2ω 0 )/P (1) (2ω 0 ), wherein β (2ω 0 ) is the total SH power scattered from a bare LiNbO 3 sphere. In good agreement with the experiment (see Figure S10), this ratio is very nearly 1 for a model with a single NP due to the weak coupling between the microsphere and the dipole of the plasmons. However, when generalizing to a model with N = 1000 NPs in the ensemble (SI, Section S.1), the enhancement ratio reaches a peak value of ∼30 in the analyzed SH frequency range and reflects the underlying resonant structure of the Mie modes. Figure 4c shows these results in detail, demonstrating the excellent quantitative agreement between the SHG enhancement theory and the experimental data. Importantly, the model provides not only a numerical reproduction of the enhancement peaks but straightforwardly explains their origin. Beginning with the lowest-order nontrivial terms in the perturbation series of the Mie resonance β, one can see that the second-order term ρ (2) β describes the transfer of energy from the upconversioninduced dipoles of the Au NP to a Mie resonance of the LiNbO 3 sphere, a process that is impossible without the presence of the NP ensemble and, as described above, represents a newly discovered pathway through which energy can be transferred from the energy bath of the laser to the LiNbO 3 particle and enhance the overall SHG signal. Accordingly, the third-order term describes the process in which energy is transferred from one Mie mode to another through two separate energy exchanges with the NP ensemble, a process which tends to shuttle energy from broad, more strongly-pumped Mie modes to narrower resonances more weakly excited by the upconversion process.
Higher-order processes can be derived from the analytical solution to ρ (n) β (ω) but contribute increasingly small corrections to the total LiNbO 3 polarization, as is shown in the SI, Section S.1. Further, after comparing the magnitudes of the forces of the system, one can conclude that the Au NP ensemble contributes to the enhancement of SHG primarily by mediating energy transfer from the broader, lower-order Mie resonances to their narrower counterparts. The second-order terms do add a small boost to SHG through the augmentations σ βν (r 0 )F 1ν (r 0 , ω) to the driving forces of each mode (SI, Section S.1), but while their total contribution is important for relatively broad modes β, it is minor for narrow Mie resonances. Figure S14 shows this effect explicitly, with the total enhancement signal well-approximated by a second-order expansion in regions of small enhancement but requiring third-order corrections near the high-order mode resonance positions.
Our second (ii) hypothesis, in contrast to the coupling model developed above, is concerned with the radiation from the Au NPs themselves. Generally, Au NPs less than ∼30 nm in diameter can be assumed to be weak scatterers as their characteristic length scale lies below one tenth of their surface plasmon resonance wavelength, implying that our choice to this point to neglect light emission from the NP ensemble is appropriate. However, small dipole oscillators such as our NPs are known to exhibit superradiance, an effect which is observed when coupling between resonant emitters causes each member of a large ensemble to emit more light than it would on its own. This phenomenon has been observed recently in collections of coupled dipole emitters in many prior investigations 43,44 and thus is of interest here.
To evaluate whether the radiation from each NP in the experimental ensemble is boosted sufficiently by inter-NP interactions to generate overall SHG enhancement from the collections of NPs in our experiments, we combine per-particle optical cross-sections collected from multiple-Mie scattering simulations 49 with analytical models to calculate the time-averaged power radiated at the SH from the "typical" NP in square ensembles of sizes N = 1 through 225. We label this power as P pl (2ω 0 ), such that the power radiated from the ensemble is N P pl (2ω 0 ). With the details left for the SI, Section S.2.3, the ratio of the powers scattered from the ensemble and the bare LiNbO 3 particle is then given by: Here, κ(2ω 0 ) is a spectral profile (Eq. [S34]) that reflects the resonance structure of the LiNbO 3 modes and has a small maximum value of ∼10 −3 . Further, A(N ) is a phenomenological scattering enhancement factor of the Au NPs extracted from the data in Figure S15a that builds in superradiance effects and corresponds to a value ∼300 for N ∼ 1000 that is not greatly influenced by ensemble disorder (see Figure S15b). Combined with the radiative contribution to the damping rate of the Au NPs γ rad 1 ∼ 10 −4 γ 1 , the mean Mie oscillator strength f 2 and damping rate γ 2 , and the rest of the previously defined constants, the superradiant enhancement factor does not overcome the small value of a 3 1 κ(2ω 0 ) to provide an overall scattering boost from the NP ensemble. Explicitly, N P pl (2ω 0 )/P (1) (2ω 0 ) 1 for all relevant wavelengths, varying between ∼2 × 10 −2 and 3 × 10 −2 for 2ω 0 in the visible spectrum. Superradiance, therefore, cannot account for the observed enhancements to SHG.

III. DISCUSSION
The invention of efficient miniaturized NLO devices remains an open and important challenge in the field of nanophotonics. Moreover, with current strategies for the maximization of NLO signal strength at the nanoscale limited by the difficulties and expense of precisely nanoengineering single NLO nano-and microstructures, multiply-resonant nanoantenna systems, and/or individual electric field hotspots, new hybrid nanomaterials that can provide meaningful NLO enhancements without unduly complicated fabrication processes are highly desirable. To this end, we have designed a new class of disordered Au-LiNbO 3 NLO nanostructures that are straightforward to synthesize and demonstrate state-ofthe-art enhancements. These hybrid structures were prepared through a high-throughput solution-phase selfassembly process through which the mesoporous structures of micrometer-scale LiNbO 3 spheres were coated with a disperse layer of ∼10-nm diameter Au NPs with a highly tunable density (∼6-1 000 Au NPs/LiNbO 3 sphere) and an even distribution.
Single nanostructures were then isolated, characterized, and modeled in order to interrogate the mechanisms governing the enhancement of their SHG. Broadband linear and nonlinear optical scattering and extinction measurements of these hybrid structures demonstrate Mielike resonant behaviors within the LiNbO 3 and minimal rearrangement of the LiNbO 3 resonance structure by the Au NPs. The Mie-like behaviors are consistent with prior investigations of multicrystalline NLO microspheres. 13 Further, analytical and numerical models predict that the coupling of each deeply subwavelength, nonradiative NP with its neighboring NPs and the microsphere is weak, such that modifications of the microsphere's spectrum should, in agreement with observations, be minimal. These models also reveal that SHG enhancement occurs within the hyrbid structures through a previously unexplored near-field coupling mechanism, available only at the nanoscale, that enhances the conversion of energy from incoming light into NLO material resonances and then confines it within the nanostructures' narrow resonances.
As a result, the hybrid materials described herein broaden the optical nanomaterial design space by highlighting the NLO enhancement capabilities of subwavelength or otherwise nonradiative nanoantenna systems. For example, the observable SH resonances of the LiNbO 3 are long-lived with γ β /ω β ∼ 10 −3 -10 −2 such that they are candidates for use in sensing and light emission applications where measurable impurity-induced peak shifts and/or tunable monochromatic emission at the nanoscale are desired. We have shown that the SH excitation of such narrow Mie resonances is preferentially and significantly enhanced by weak coupling to a relatively low and disordered volume of a second resonant material (Au) that is resonant at the output SH frequency of the LiNbO 3 sphere. Therefore, new design strategies for sensors and emitters can specifically incorporate imprecise ensembles of very small nanoantennas as needed to optimize the NLO signal intensity.
In addition, the sensitivity of the SH Mie resonances' output frequencies to intrinsic and environmental defects is not degraded by the presence of the Au NPs. As evidence, the strong dependence of the microsphere's Mie resonance spectral positions, widths, and heights on its geometry and dielectric function is preserved in hybrid structures, as is shown in Figure 3b. This preservation is supported by the excellent agreement with the model, which specifies that the weak coupling between the Mie resonances and plasmon dipoles precludes rearrangement or obfuscation of the LiNbO 3 Mie spectrum.
Determination of the exact value of the dielectric function 2 of the structures analyzed in the experiment is difficult, as the mesoporous material is novel, single-particle ellipsometry through e.g. electron beam spectroscopy 50 is in its infancy, and substrate effects can cause nontrivial reorganization of the sphere's Mie spectrum (SI, Section S.2.2). However, as is detailed in the SI, Section S.2.1, the estimate 2 = 5.5 + 0.035i used alongside the chosen value a 2 = 500 nm as well as the HAADF images and linear and nonlinear scattering spectra of Figures 1, 3, and 4 produces uniquely precise agreement between theory, experiment, and prior art 45 in the SH spectral region. Variation of either a 2 or 2 by ∼5% can generate noticeable ∼15-35 nm shifts in the Mie resonance positions of both bare LiNbO 3 and hybrid structures (see Figure 4d and S14c), such that we take ±5% to be a rough estimate of the uncertainty in the microsphere parameters. The small (∼25 nm) phenomenological shifts made to align the theory with the experimental data of Figure 4c fall well within these bounds.
By contrast, the Mie resonances are not particularly sensitive to the parameters of the NP plasmon resonances. The NPs' extinction profiles simply provide a spectral envelope inside which SHG enhancements are possible. Finally, the microsphere's surface coverage by the NP ensemble is limited, allowing for ready interaction between external NPs, cavities, waveguides, molecules, etc. and the LiNbO 3 .
Future advancements in the fabrication and characterization processes of our (or similar) nanostructures that provide greater control over simple parameters like the microsphere's shape and internal nanocrystal density will allow researchers to more precisely tune these systems. Further, the simplicity of the coupling phenomena governing the SHG output intensity enhancements suggests that nanostructures displaying NLO behaviors beyond SHG should be able to take advantage of the same effects. In combination with the flexible perturbative model developed here, our results suggest that nanoscale NLO enhancements can be straightforwardly modeled, designed, and realized in a wide array of nearfield coupled nanoscopic systems displaying all manner of frequency conversion behaviors without the need for expensive nanoengineering.

IV. MATERIALS AND METHODS
Further details regarding the synthesis, characterization, and modeling of the nanostructures used in this work can be found in the SI.
Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Information. After cooling to room temperature, white precipitates were isolated from the solution via a process of centrifugation (Model No. AccuSpin 400, Fisher Scientific) at 8,000 rpm for 15 min and decanting of the solution. These solids were washed by re-suspending them in 10 mL deionized water (18 MΩ·cm, produced using a Barnstead NANOpure DIamond water filtration system). This purification process was repeated for a total of three times. The purified product was dried at 70 • C for 10 h to remove residual water prior to further analyses. The dried precipitates were calcined by heating from room temperature to 600 • C at a rate of 5 • C/min and held at 600 • C for 45 min to induce complete crystallization.

S.1.1.2 Synthesis of Gold-Lithium Niobate Hybrid Nanostructures
An aqueous suspension of porous, monodisperse LiNbO 3 particles (0.5 mg/mL) was added and dispersed in 5 mL of water via sonication. Into this suspension, either 1 mL (for a higher loading of gold nanoparticles) or 0.1 mL (for a lower laoding of gold nanoparticles) of a 5 mM aqueous solution of gold(III) chloride trihydrate (HAuCl 4 · 3H 2 O, 99.9%, Sigma Aldrich) was added, and the mixture stirred at 70 • C for 3 h. This step was followed by the addition of 1 mL or 0.01 mL of 0.5% (w/v) trisodium citrate dihydrate (C 6 H 5 Na 3 O 7 · 2H 2 O, Sigma Aldrich, ≥99%) aqueous solution into the reaction mixture. The reaction mixture was stirred further at 70 • C for 1 h. The precipitates were washed two times with water. The as-synthesized hybrid nanostructures of LiNbO 3 with gold (Au) nanoparticles (NPs) were isolated from the unreacted excess gold salts and unbound Au NPs through centrifugation (Thermo Electron Corporation, IEC microlite microcentrifuge) at 2,000 rpm for 3 min, decanting of the supernatants, and re-dispersion of the isolated solids in DI water with the assistance of a vortexer for 3 min. The morphology and dimensions of the LiNbO 3 particles were characterized using an FEI Osiris X-FEG 8 scanning/transmission electron microscope (TEM/STEM) operated at an accelerating voltage of 200 kV. Samples for TEM/STEM analyses were prepared by dispersing the purified products in ethanol followed by drop-casting 5 µL of each suspension onto separate TEM grids (300 mesh copper grids coated with Formvar/carbon) purchased from Cedarlane Labs. Each TEM grid was dried at ∼230 Torr for at least 20 min prior to analysis. Energy dispersive X-ray spectroscopy (EDS) analyses were performed using the FEI Osiris scanning/TEM, which was equipped with a Super-X EDS system with ChemiSTEM Technology integrating the signal from four spectrometers.
Purity, crystallinity, and phase of the LiNbO 3 particles were characterized using Raman spectroscopy and powder X-ray diffraction (XRD) techniques. Raman spectra were collected using a Renishaw inVia Raman microscope with a 50× SWD objective lens (Leica, 0.5 NA), and a 514 nm laser (argon-ion laser, Model No. Stellar-Pro 514/50) set to 100% laser power with an exposure time of 30 s. The Raman spectrometer was calibrated by collecting the Raman spectrum of a polished silicon (Si) standard with a distinct peak centred at 520 cm −1 . The Raman spectra for the samples were acquired from 100 to 1,000 cm −1 using a grating with 1,800 lines/mm. The XRD patterns of the samples were acquired with a Rigaku R-Axis Rapid diffractometer equipped with a 3 kW sealed tube copper source (Kα radiation, λ = 0.15418 nm) collimated to 0.5 mm. The samples were packed into a cylindrical recess drilled into glass microscope slides (Leica 1 mm Surgipath Snowcoat X-tra Micro Slides) for acquiring XRD patterns for the products.
The optical absorption spectra of the products were measured using an Agilent Technologies ultraviolet-visible spectrophotometer (Agilent 8453, Model No. G1103). For these measurements, the samples were suspended in water and held in 1 cm path length poly(methyl methacrylate) cuvettes (VWRTM, Catalog No. 634-8537). The linear spectra of individual pristine LiNbO 3 and hybrid Au-LiNbO 3 particles were measured using a Zeiss M1m optical microscope operating in a dark-field spectroscopy setup. The light from a halogen lamp was focused onto the sample through a dark field objective (50×, Zeiss Epiplan Neofluar). The scattered light was collected by the same objective and detected by an imaging spectrometer (Princeton Instrument Acton spectrometer equipped with a PIXIS 400 CCD detector cooled to -72°C). The linear scattering cross section spectra were obtained by normalizing the signal from a single particle to the signal coming directly from the lamp, after subtraction of the background in proximity to the particle.
The second harmonic generation (SHG) activity of the individual pristine LiNbO 3 and hybrid Au-LiNbO 3 particles was assessed using a Leica SP5 laser scanning confocal two photon microscope equipped with a Coherent Chameleon Vision II laser and a Zeiss LSM 510 MP confocal microscope. A dilute dispersion of nanostructures was drop-cast onto glass coverslips and brought into the focal point of the microscope. The SHG response was characterized by locating individual particles using a 63× oil immersion objective aperture and scanning the laser excitation from 850 nm to 1,070 nm as the fundamental wavelength.

S.1.2.1 Bare Lithium Niobate Sphere Characterization
Modeling the microsphere as a smooth, isotropic spherical particle, the nearly constant dielectric function of LiNbO 3 in the optical region of interest 1 allows the microsphere's linear response near the SH to be described by the excitation of a set of 21 Mie resonances (SI, Section S.2.1.2). We have established β = {T , p , , m } as a collective index to describe these modes such that the scattered electric field of each mode can be expanded as Using E sca along with the magnetic fields of each mode, B sca (r, ω) = (c/iω)∇×E sca (r, ω), the scattered power averaged over a period of the SH frequency can be calculated from Poynting's theorem as (Section S.2.3.1): Here, α β (ω) = (a +2 2 f β /2ω β ) exp(iψ β )/(ω β −iγ β /2−ω) is the linear polarizability of the mode β, with C β (ω 0 ) ∼ E 0 χ 2 (ω 0 , ω 0 ) a unitless overlap coefficient. In general, C β (ω 0 ) quantifies the efficiency of photon upconversion between the fundamental resonances of the system at ω 0 and a resonance β at the SH, and thus is large only if the fundamental resonances are driven strongly by the incident light. Also, C β (ω 0 ) is zero unless the fundamental and SH resonances satisfy specific symmetry requirements, which are detailed below. Finally, the explicit form of this overlap coefficient is complicated, but described in full in Section S.2.4.2, Eq. (S56).
In general, the resonance frequencies of the Mie modes are independent of p and m, although these indices do determine the strength of the response of a mode to a driving source of a given spatial symmetry. In detail, if we allow α = {T, p, , m} to be the collective index of one fundamental resonance of the LiNbO 3 particle and α = {T , p , , m } to be the index of another, in the case where the sphere has a weak, isotropic, and frequencydependent nonlinear susceptibility χ 2 (ω , ω − ω ) only modes β for which p + p + p is even and m ± m ± m = 0 are driven by SHG. For example, with the microsphere driven by an x-polarized plane wave at the fundamental frequency such that α = {E, 0, , 1} or {M, 1, , 1}, each of the 21 SH modes in the model has an index pair (p , m ) = (0, 0), (0, 2), or (1,2). In contrast, the dependence of a mode's response on T and is determined by its spectral overlap with the source. In the energy window of the observed SHG enhancement, seven sets of modes can be significantly driven by the upconversion process: four sets of electric modes with = 6, 7, 9, and 10, and three sets of magnetic modes with = 7, 10, and 11. See Figure S13 for details.

S.1.2.2 Characterization of the Lithium Niobate Radius and Dielectric Function
Section S.2.2 details the characterization of each of the LiNbO 3 Mie modes, each of which is assigned a resonance frequency ω β , a damping rate γ β , and an oscillator strength f β (equivalently, an effective mass µ β = e 2 /a 3 2 f β ) that determine its spectral position and response magnitude to external stimuli. This characterization is performed using a dielectric function 2 = 5.5 + 0.035i and a radius a 2 = 500 nm estimated from the experiment. The real part of 2 is taken from estimates extracted from a set of nine single-particle scattering experiments conducted on bare LiNbO 3 spheres of radii 350 to 1000 nm. Figure 3c shows results from one representative scattering experiment, with the remainder of the results given Figure S17 and in Table S3. Further, in conjunction with the choice of radius, the choice Re{ 2 } = 5.5 provides the best overlap between the spectral locations of the Mie resonances and the SHG enhancement peaks seen in Figure 4c, as well as the best agreement between the relative peak heights. We note that small (∼ 5%) increases (decreases) to Re{ 2 } and a 2 lead to red (blue) shifts in the Mie resonance positions that can nullify each other, such appropriate choices with the ranges a 2 ± 20 nm and Re{ 2 } ± 0.1 are likely to lead to similar results. The inset of Figure 4c shows the range of Mie resonance energies available with ±4% changes to a 2 , and analogous results when varying Re{ 2 } are shown in Figure S14c.
The imaginary part of 2 , which is expected to be small 1 , is chosen to agree well with the linewidths of the modes observed in the second harmonic data of Figure 4c. The exact value has yet to be characterized via e.g. ellipsometry 3 and is difficult to estimate from the scattering data due to substrate effects (see Figure S13a). The linewidths of the Mie resonances depend sensitively on the rate of internal material losses (see Figure S14d) and the heights of the narrow Mie modes in both scattering and SHG enhancement spectra can vary noticeably with ∼ 5% changes to Im{ 2 }. However, small changes to Im{ 2 } do not modify the peak positions within the SH enhancement spectrum such that we take Im{ 2 } to be a fitting parameter that affects the overall magnitude of the theoretical SHG enhancements, but does not qualitatively affect the results. Further details are given in Section S.2.1.

S.1.2.3 Superradiant Nanoparticle Scattering
Each NP is assumed to be spherical and isotropic with a Drude-Lorentz dielectric function 1 (ω) fit to the Au dielectric data of Ref. 4 (Section S.2.1, Figure S16, and Table S4) and an associated set of surface plasmon resonances within the spectral window of the observed SH emission. Figure 3a demonstrates the excellent agreement between the resulting theoretical absorption cross section of the NP dipole plasmons using this model and extinction measurements performed on single hybrid nanostructures, highlighting the validity of the model and the uniformity of the NP responses in the experiment. We thus restrict our analysis to include only the dipole resonances of the NPs, which we label d i (ω), and allow each particle to lie at a position r i on the surface of the LiNbO 3 microsphere.
The dynamics of the plasmon dipoles are inferred from numerical scattering calculations using both multiple-Mie scattering and boundary element method Maxwell solvers 5,6 . These simulations were performed on arrays of 10-nm diameter Au NPs in square grids with a 20 nm center-to-center spacing, as complete layers on spherical templates containing the N ≈ 1000 Au NPs present in the experiment are too large to compute directly. Moreover, grids reliably approximate the effects of interparticle coupling in the limit that the radii of the Au NPs and their separations are much smaller than the radius of the LiNbO 3 microsphere. Using these scaled down simulations, we extrapolated from smaller ensembles (N ≤ 225) the behaviors of larger collections of Au NPs. For comparative reasons, we used both regular grids and randomized grids of Au NPs to evaluate the effects of disorder on the amount of light scattered by each particle.
Our results ( Figure S15a) show that the maximum scattering cross section of each dipole grows roughly two orders of magnitude from 7.75 × 10 −3 nm 2 as N varies from 1 to 1000, and that the exact magnitude of the growth depends on whether the dipole is oriented perpendicular (maximum of 0.355 nm 2 ) or parallel (1.18 nm 2 ) to the plane of the ensemble. Conversely, the absorption cross section maxima of each dipole converge quickly with ensemble size, as shown in Figure S15a, moving from a single particle value of 13.5 nm 2 to perpendicular and the parallel values of 10.3 nm 2 and 15.7 nm 2 , respectively, for ensembles of N > 50 NPs. In addition, while randomness in the positions of the particles can modulate the per-particle scattering cross section maxima by ±20% in comparison to a regular grid ( Figure S14b), an extreme value distribution fit to the recorded scattering cross section maxima of 100 simulations with randomized grids of N = 16 NPs shows that modifications beyond this range occur in ∼0.4% of NPs.
The regular-grid simulations are, thus, taken to be representative of the ensemble as a whole. Further, Figure S15b shows that while the peak maxima of the scattering observables depend strongly on N , their lineshapes do not. Therefore, with the NPs' scattering strongly modified and their absorption changed relatively little with increasing N , we conclude that the dipole oscillator strengths are best estimated from the absorption data. Explicitly, each dipole is then treated quasistatically with a position-and orientation-dependent polarizabil-ity (see Section S.2.2.2), while the ensemble-induced scattering enhancements are included phenomenologically through modification of the Rayleigh scattering rate of each dipole by a positive factor A(N ). When driven by the scattered fields of a bare LiNbO 3 sphere, the ensemble then radiates a time-averaged scattered power Here, f 1 is the orientation-averaged oscillator strength of the NP dipoles (Eq. [S13]), and K β (Eq. [S29]) is an overlap coefficient that determines the strength of the ensemble's response to the LiNbO 3 fields. See Section S.2.3 for further details.

S.1.2.4 SHG Enhancement via Near-Field Coupling
The first and second-order terms in the solution of the Mie resonance equations of motion in Eq. (1) are which detail the driving of each mode β through direct upconversion and through excitation of the NP by the upconversion process, respectively. Higher-order terms build in the consequences of multiple transfers of energy between the microsphere and the NP and become increasingly tedious to write explicitly, but can be analyzed compactly with from which it is clear that, with |ρ β (ω)|, terms beyond the third order provide only small corrections.
To account for all N NPs, we label each NP with the index i ∈ [1, N ] such that r 0 → r i and d ν (ω) → d iν (ω), but allow the NPs to have identical radii and their plasmons to have identical resonance frequencies, linewidths, phase offsets, and oscillator strengths. Therefore, the coupling term in Eq. (1) is replaced by the sum iν σ βν (r i )d iν (ω) and similar sums in each successive term of the perturbation expansion must be adjusted accordingly. In this work, these sums are carried out by placing three dipoles d iν (ω) at N = 1,000 points arranged on a sphere of radius a 1 + a 2 in a Fibonacci lattice 7 . Calculations with randomized points display only small variations from the results of the regular Fibonacci array in good agreement with the random-ensemble scattering results of Figure 3b and are thus not shown.

S.2.1 Construction of the Material Models from Data
In order to maximize the agreement between our theoretical models and the SH enhancement data, we construct a model dielectric function for Au and a dielectric function and secondorder susceptibility and LiNbO 3 using our own single particle scattering data, as well as from ellipsometry data and atomistic numerical simulations from Refs. 4 and 8, respectively. The Au details are given in Section S.1.2.2 and the LiNbO 3 details are given in Section S.2.1.2.

S.2.1.1 Fitting of the Au Lorentz-Drude Dielectric Model
The oscillator parameters of the dipole plasmons within the Au NPs modeled in this investigation are calculated using a simple Drude-Lorentz model of the dielectric function of gold. This dielectric function is given by and its form and parameters are inferred from ellipsometry data published in Ref. 4. In particular, as shown in Figure S16a, a Drude-model (ω p2 , ω p3 → 0) fit via nonlinear least squares methods to the gold dielectric function provides a good approximation to the dielectric data for photon energies ω < 2.0 eV, but this model cannot capture the effects of interband transitions that become increasingly important at higher energies. Indeed, a Drude model is provided by Ref. 4, but it is not used in this analysis as it produces inaccurate predictions of the dipole plasmon energies of Au NPs, which exist well above the 2.0 eV threshold. To produce a more accurate Au NP model in the region between 2.0 and 2.8 eV, the addition of two Lorentz oscillators to the Drude model dielectric function was found to be sufficient. The NP dipole plasmons have resonance energies near 2.5 eV, well below the upper bound of the single-oscillator Drude-Lorentz model's accuracy near 2.9 eV (see Figure  S16a). The parameters for the Drude-Lorentz dielectric model were inferred by first fitting a simplified Drude model to the real and imaginary parts of single-crystal dielectric data of Ref. 4 at low energies ω Λ 2 . These fits produced estimates ω p1 = 7.20 ± 0.9 eV and Γ 1 = 71.1 ± 13 meV, respectively, which were then used as initial guesses for fits using the more robust Drude-Lorentz model. Focusing on the energy window between 0.8 and 2.8 eV, the Drude-Lorentz fits were made using the regularized function Re{ 1 (ω)} + η(ω)Im{ 1 (ω)} to simultaneously fit the model to both components of the complex dielectric data while avoiding bias toward the data's much larger real part in the energy window between 0.8 and 2.4 eV.
Fits using simplex, differential evolution, simulated annealing, and random search fitting methods 9 produced stable results using weighting functions η(ω) = λ 1 exp(−λ 2 ω) with λ 1 = 40 ± 5 and λ 2 = (6.5 ± 1) × 10 −16 s/rad. The central values for either factor are used in this work. Taking the average of the results of the four fitting routines, the model parameters used in the following discussion are given in Table S4.

S.2.1.2 Fitting of the Lithium Niobate Dielectric Model and Nonlinear Susceptibility
The radii of the samples examined in Figure 4c have radii of a 2 ≈ 500 nm, and we use this value in the following analysis. To estimate 2 , nine scattering spectra were collected using bare LiNbO 3 particles with radii between 350 nm and 1000 nm, with the radii estimated from microscopy analyses with a precision of ±50 nm. Mie theory scattering spectra were compared to the experimental data, but peak splitting and broadening generated by the substrate in the experimental data restricts the utility of the ideal substrate-free model. More precisely, in agreement with Ref. 10, the presence of a substrate in the scattering spectra generates peak broadening that can in most cases be captured to acceptable approximation by a Mie model with phenomenologically added internal damping. However, as is also demonstrated in Ref. 10, the substrate can also induce splitting, shifting, and even peak amplification that are not uniform and are not reproducible by simple Mie theory. These effects can be seen in Figure S17, where single peaks in a Mie theory model are spectrally aligned with multiply-peaked features in the data, and more obviously in Figure S13. The latter compares Mie theory to simulations that include a substrate, showing that a phenomenologically damped dielectric model reproduces some substrate effects (peak suppression, broadening) but not others (enhancements, shifts of narrow peaks).
Thus, estimates of the real part of 2 were performed by aligning the visible broad peaks of the theory and experiment in the spectral window between 720 nm and 380 nm, and the absolute peak widths and relative heights were ignored. Narrow features in the data are also ignored as they are most likely to be shifted by the substrate in a manner not captured by Mie theory, as evidenced by the narrow feature near 2.5 eV in Figure S13a. Results using these comparison guidelines are shown in Table S2.
Within the range of possible values of Re{ 2 } shown therein, the value 5.5 provides the best agreement between the spectral positions and relative magnitudes of the SH peaks in the theory and experiment, assuming a 2 = 500 nm. As is mentioned in the main text and demonstrated in Figures 4c (inset) and S14c, small tweaks to 2 and a 2 will provide similar Mie spectra such that the particular values chosen within a narrow range (± ∼ 5%) are irrelevant to the conclusions of this investigation.
The imaginary part of 2 is taken to be 0.035. This value is likely larger than the true value extracted from earlier ellipsometry experiments 1 , as the linewidths of the Mie resonances of the LiNbO 3 sphere are increased by the presence of a substrate. Because the magnitude of this increase is difficult to infer from first-principles models, we instead choose Im{ 2 } to agree well with the linewidths of the modes observed in the second harmonic data of Figure  4c,b. Within a range of ± ∼ 10% of the chosen value, Im{ 2 } simply scales the magnitude of the SH enhancements, such that, like the choices of Re{ 2 } and a 2 , the precise value of Im{ 2 } does not affect our conclusions.
The second-order susceptibility of the LiNbO 3 , χ 2 (ω , ω − ω ) was fit to the simulated data from Ref. 8 using an anharmonic oscillator model 11 in the approximation that the LiNbO 3 has an isotropic, homogeneous response within the microsphere boundaries. Explic-itly, the susceptibility is given by: wherein Λ and Γ are the natural frequency and linewidth of the Lorentz oscillator, respectively. These features are used to characterize the response of the LiNbO 3 carriers. χ (2) ∞ is a constant offset that approximates the shift of the real part of χ (2) 2 at low energies by higherenergy transitions. Further, s is the characteristic anharmonic oscillator strength with units cm · s 6 /statV and is taken to be very small such that sE 0 Λ 6 with E 0 the characteristic strength of the laser field. Finally, 1 3 is the 3 × 3 identity matrix.
To fit the data, we evaluated the second-order susceptibility in the limit that ω = 2ω , i.e. with the assumption that the driving laser at ω = ω 0 is very narrow and the observation frequency is always at 2ω 0 . With this restriction, the expression for the susceptibility simplifies to χ 2 and can be easily visualized as is shown in Figure S11.
As was done with the dielectric function of Au in Section S.1.2.2, the model was fit to the real part, imaginary part, and absolute value of the data using simplex, differential evolution, simulated annealing, and random search nonlinear least squares methods. An average of the results of each of the four methods was collected for the fit to each function of the data, and the parameter average that produced the best fit to the data in the experimentally relevant energy range was selected. In this case, the fits to the absolute value of the data were superior in the region between ∼1.0 eV -1.5 eV analyzed in Figure 4, returning parameter values of Λ = 4.00 eV, Γ = 748 meV, s = 4.08 × 10 −3 cm · s 6 /statV, and χ The estimation of the dielectric function of the LiNbO 3 is described in the main text. We note here that the corresponding linear susceptibility model of LiNbO 3 to the second-order function described above is a Lorentz-model dielectric χ (1) However, in the case where Λ is sufficiently detuned from the fundamental and second harmonic frequencies, χ (1) (ω) is well-approximated as a dielectric. In our case, Λ− Γ > 2ω 0 such that the constant-dielectric approximation is appropriate and in agreement with Ref.

S.2.2 Estimation of the Mode Oscillator Parameters
In this section, we use the dielectric models we have constructed to infer the parameters for oscillators models of the plasmon and Mie resonances of the coupled Au-LiNbO 3 nanostructures. Bare Au and LiNbO 3 particles are discussed in Sections S.2.2.1 and S.2.2.3, while the effects of weak coupling between members of the NP ensemble are discussed in Section S.2.2.2.

S.2.2.1 Bare Au NP Dipoles
The relation of the oscillator parameters of the Mie and plasmon resonances to the dielectric parameters of their respective particles is done via the particles' response functions. These are shown explicitly in Eqs. (S45)-(S49) (see Section S.2.4.2), and provide a direct link between the Au and LiNbO 3 dielectric functions 1 (ω) and 2 , respectively, the resonance freqeuncies ω, damping rates γ, masses µ, and phase offsets ψ of each mode, and the physical observables of the particles.
For example, with 1 (ω) characterized in Section S.1.2.2, the response function of Eq. (S48) can be used used to construct the absorption cross section, of each of the Au NP dipole modes, in which the plasmon oscillator parameters (subscript 1) and interband resonance parameters (subscript L i ) can be exactly expressed as functions of the dielectric parameters of 1 (ω). These functions are impossible to write explicitly, as the eigenfrequencies Ω 1,L i = ω 1,L i − iγ 1,L i /2 are solutions to a sextic polynomial. Nevertheless, the oscillator parameters are simple to infer numerically by fitting Eq. (S8) to the standard form σ abs 1 (ω) = (4πω/c)a 3 1 Im{[ 1 (ω) − 1]/[ 1 (ω) + 2]} as shown in Figure S12b with the NP radius a 1 set to 5 nm.
It is well-known that plasmon resonance energies are redshifted when the metal nanoparticle supporting them comes in close contact with a dielectric substrate with a large refractive index 12 . As the LiNbO 3 spheres in this investigation are much larger than the Au NPs and have a constant dielectric function larger than 5 (see Section S.2.2.3), some plasmon redshifting is expected. We incorporate the redshifts by lowering the Au NP dipole by 0.23 eV to align the NP absorption maximum with the observed maximum of hybrid-particle extinction measurements in Figure 3a. The oscillator parameters used in this work are given in Table  S1. However, to model the effects of radiation on the plasmon motion, one must modify the results extracted from Eq. (S8). More concretely, in the formal quasistatic approximation under which the absorption cross section is derived, only the rate of Au carrier losses to heat will appear in the cross section expression and any back-action of radiation on the dipole plasmon will be ignored. This back-action generally results in broadening of the plasmon lineshape, such that one can phenomenologically expand γ 1 → γ rad 1 + γ NR 1 to provide the plasmon damping rate with both a radiative (rad) and nonradiative (NR) component.
While approximate, this simple formulation of scattering loss works well in the limit γ rad 1 γ NR 1 . The validity of the approximation can be demonstrated by comparing the Rayleigh scattering cross section in which the interband damping rates have been similarly transformed as γ L i → γ NR L i + γ rad L i to account for radiation broadening and red shifting, to the exact Mie theory scattering cross section of a spherical Au NP with a 5-nm radius. Numerical fits of the modified oscillator model to the Mie scattering lineshape are generally unstable and suggest radiationinduced changes to both γ 1 and γ L i are roughly 1% or smaller, so we instead approximate the radiative damping rates of the dipole plasmons and interband resonances using their Larmor formulae. Assuming γ NR 1 ≈ γ 1 , we have γ rad 1 = 2e 2 ω 2 1 /3c 3 µ 1 = 1.01 × 10 −4 γ NR 1 . The radiation from the Lorentz oscillator resonances is similarly minimal, with γ rad Figure S12a shows the excellent agreement between the Larmor-modified oscillator model and Mie scattering calculations.

S.2.2.2 Ensemble-Modified Au NP Dipoles
To describe in a concise way the scattering observables of the dipole plasmons in the NPs surrounding the LiNbO 3 sphere, it is first necessary to quantify the alterations to the scattering behaviors of each NP caused by its neighbors. In other words, due to inter-NP interactions, we cannot take the masses µ i (or, equivalently, the oscillator strengths f i = e 2 /a 3 1 µ i ) of the N NPs of the ensemble to be their bare values µ 1 (or e 2 /a 3 1 µ 1 ). We can begin this process by letting the polarizability of each NP be a tensor dependent on its position, such that differences between the ensemble perpendicular and parallel dipole oscillations of each NP can be captured. Explicitly, we let are simply the spherical unit vectors evaluated at the angular position of the i th dipole. The tensor elements α ⊥, are the polarizabilities of the dipole components at r i that are oriented perpendicular and parallel to the shell of Au NPs, respectively, and have magnitudes characterized by the oscillator strengths f ⊥, and phases determined by the angles ψ ⊥, .
We can define the terms of α i (ω) with poles at ±Ω 1 as and α ⊥, (ω) = a 3 1 f ⊥, e iψ ⊥, /2ω 1 (Ω 1 ∓ ω), respectively. Therefore, upon excitation by an impinging field E(r i , ω), the dipole set up in the i th NP is Here, E(r i , ω) is the (real) magnitude of the electric field at r i ,ˆ i is the complex polarization unit vector that describes the phases and orientations of the field components, and Θ(ω) is the Heaviside function. Letting d (±) (ω,ˆ i ) be the dipole terms valid at positive and negative frequencies, one finds the magnitude of d (+) is We can, therefore, define as the orientationdependent dipole plasmon oscillator strength.
Withˆ i oriented strictly perpendicular to the ensemble, one can see that Figure S15c shows the simulated absorption cross section maxima of NPs in ensembles excited with these two choices ofˆ i , where it can be seen that for N 50 the absorption cross section of each particle is independent of N . Analytically, the absorption cross sections of ensemble perpendicular and parallel dipoles are simply σ abs ⊥, (ω) = (4πω/c)R abs Im{α ⊥, (ω)}, where R abs = 1.16 phenomenologically accounts for the absorption contribution of the Lorentz oscillators in the region ω < 2.8 eV. With cos ψ ⊥, ≈ cos ψ 1 ≈ 1, the cross section maxima are then σ abs ⊥, (ω 1 ) ≈ 4πR abs f ⊥, a 3 1 /cγ 1 such that the oscillator strengths implied by Figure S14c are f = 2.37 × 10 30 s −2 and f ⊥ /f = 0.668.

S.2.2.3 Bare Lithium Niobate Mie Resonances
Inference of the oscillator parameters of the Mie resonances β of the LiNbO 3 microsphere is more straightforward, as their damping rates γ β contain only losses to radiation. Thus, only a single observable is needed to extract the full set of parameters from each mode. We use the response functions shown in Eq. (S49) to fit γ β as well as the resonance frequencies ω β , effective masses µ β (equivalently, the oscillator strenghts f β = e 2 /a 3 2 µ β ), and phase offsets ψ β . Explicitly, we let (S14) Note that, in contrast with Eq. (S49), the above expressions only consider a single resonance for a given β. See Section S.2.4 for details. The results of the fits are shown graphically in Figure S15b,c and are tabulated in Table S2.

S.2.3 Derivation of the Scattering Observables
We provide here a more thorough derivation of the superradiant scattering enhancement ratio. First, in Section S.2.3.1, we describe the scattered fields and power of the modes of a dielectric sphere with a diameter on the order of a wavelength of the scattered light. Second, in sections Section S.2.3.2 we develop the scattering observables of the ensemble of NPs surrounding a dielectric sphere. In Section S.2.3.3, we complete the calculation of the ratio of scattered powers of a NP-dressed and bare LiNbO 3 dielectric sphere.

S.2.3.1 Mie Modes of the Lithium Niobate Sphere Driven by a Plane Wave
The time averaged radiated power from each of the LiNbO 3 SH modes can be calculated from Poynting's theorem. Beginning with wherein (2) signifies that the fields arise from a second-order scattering process, ρ β (ω) are the magnitudes of the moments of the modes β, and X β (r, ω) are the regularized vector spherical harmonics for T = M and E, respectively, the calculation of the Poynting vector is straightforward.
Here, h (1) (x) are the spherical Hankel functions of the first kind, P m (x) are the associated Legendre polynomials, S p (x) = cos(x)δ p even + sin(x)δ p odd , and k = ω/c. Further, with from Faraday's law, and the time-averaged scattered power from the dielectric sphere can be rapidly simplified from wherein the surface integral is taken to be across a sphere of radius r → ∞, · τ is the time-average operator over a period τ , and the Fourier transforms of the fields F sca (r, t) = F sca (r, ω) exp(−iωt) dω/2π have been used. Letting the moments of the sphere be driven by the second-harmonic upconversion of an incoming plane wave of frequency ω and electric field strength E 0 , we can say ρ β the linear polarizability of the β th mode and C β (ω) an overlap coefficient defined in Eq. (S56). The final result is Importantly, the radiated power from the Mie modes contains no cross-terms due to the orthogonality condition imposed by Eq. (S18). Further, due to the good agreement between an oscillator model of each mode and its electromagnetic response (see Section S.2.2), the scattered power can also be modeled mechanically. To do so, it is important to first define a generalized coordinate q β (ω) to represent the displacement magnitude of the moments ρ β (ω) that obeys the reality condition q β (−ω) = q * β (ω). We will choose the definition is a real quantity. With q β (t) defined, we model losses to radiation as a weak, frequency-independent linear damping process. In this case, where A β is a unitless proportionality constant that relates the mechanical quantities to the electromagnetics. This leads directly tō via Eq. (S19) and the identity q β (t) = (−iω) exp(−iωt)q β (ω) dω/2π such that A β = (1 − δ p1 δ m0 )(e 2 c 3 /2πω 3 0 a 6 2 µ β ω β γ β ) = (1 − δ p1 δ m0 )(f β c 3 /2πω 3 0 a 3 2 ω β γ β ). Finally, to arrive at Eq. (S1) and the notation of the main text we substitute the explicit form of A β and make the simplificationP 2 (2ω 0 ) → P (1) (2ω 0 ) to the notation.

S.2.3.2 Plasmon Dipoles Driven by the Scattered Lithium Niobate Electric Field
The radiated power by the ensemble of Au NPs can be straightforwardly derived using a mechanical model parameterized using the values outlined in Sections S.2.2.1 and S.2.2.2. More explicitly, the time-averaged scattered power by the i th dipole plasmon of the Au NP ensemble when the collection is driven by an external field can be calculated in a straightforward manner using a well-known mechanical model of dipole radiation wherein x i (t) = d i (t)/e is the coordinate characterizing the magnitude of the dipole plasmon located at r i and oriented alongx i and A i (N,ˆ i ) is a phenomenological enhancement factor that builds in the ensemble-induced scattering enhancements seen in Figure 4c. Similarly to the orientation-dependent oscillator strengths, we take A i (N, where A ⊥, (N ) are the enhancements experience by ensemble-perpendicular and parallel dipoles, respectively, in an ensemble of N NPs.
With the radiation back-force F rad i (t) treated as a damping force, one can let F rad the mass of the i th plasmon, such that an average over a period τ = 2π/2ω 0 of the second harmonic frequency gives: Inserting the identityẋ i (t) = (−iω)x i (ω) exp(iωt) dω/2π twice and letting x i (ω) = α i (ω) · E (2) sca (r, ω)/e, one finds Application of the identity of Eq. (S19) and neglect of terms proportional to ρ β (−2ω 0 ) ρ β (2ω 0 ) provides Finally, we average over the power scattered from each dipole i to define the power scattered from the typical dipole. Numerically, one can show that terms proportional to are at least two orders of magnitude smaller when β = β than when β = β . Thus, we can safely neglect the cross terms of the sum and defineP 1 (ω) = P i (ω) i such that (S28) Rectification of Eq. (S28) with Eq. (S2) and the notation of the main text can be achieved by substituting appropriately the dimensionless constant and simplifying the notation such thatP 1 (2ω 0 ) → P pl (2ω 0 ). Here, we have taken similar to the definition of the angle-averaged oscillator strength f 1 . Thus, all that is left to characterize inP 1 (2ω 0 ) are the phenomenological constants A ⊥, (N ). This can be achieved by replacing E sca (r i , ω) with a field E 0 (r i , ω) = E 0 π[δ(ω − 2ω 0 ) + δ(ω + 2ω 0 )]ê ⊥, withê ⊥, =r i andθ i , respectively, in Eq. (S26). The results can then be compared to the simulations detailed in Figure 4c, which provide the simulated scattering cross sections of a grid of spheres upon which a plane wave polarized parallel or perpendicular to the ensemble plane is impinged, such that only ensemble-parallel or perpendicular dipoles are excited. The use of E 0 (r i , ω) in the theory faithfully approximates the simulated arrangement in the limit where the spacing between NPs is much smaller than the impinging light wavelength and the LiNbO 3 sphere radius a 2 .
Explicitly, the power scattered by the i th dipole is and, after an angular average, where the superscript PW signifies that the driving source in a plane wave andP PW 1 (2ω 0 ) = P PW i (2ω 0 ) i . The associated scattering cross section toP PW 1 (2ω 0 ) that can be directly compared to the simulation is given by wherein R sca = 6.00 is a scaling constant that builds in the contribution of the Lorentz resonances of Au at low energies. With all other constants already well-characterized, we find the best fits to the simulated scattering data are A ⊥ (N ) = 2.71N 0.489 and A (N ) = 2.43N 0.730 such that A(1000) = 301. Finally, we find that due to the similar average magnitudes over r i of the harmonics X β (r i , k β ), K β varies little between modes, with max β {K β } = K M,0,10,2 = 10.2 × 10 −3 and min β {K β } = K E,0,10,0 = 5.04 × 10 −3 .

S.2.3.3 Ratio of Time-Averaged Scattered Powers
One can see from Eqs. (S1) and (S2) that the superradiant scattering enhancement ratio NP 1 (2ω 0 )/P 2 (2ω 0 ) is greatly simplified by defining the spectral profile the associated spectral profile wherein f 2 = f β β and γ 2 = γ β β are the average Mie oscillator strength and damping rate, respectively. As K β is small for all of the relevant modes of the system, so too is κ(ω) for ω in the optical region, with a value between 1.67 × 10 −3 and 1.78 × 10 −3 .

S.2.4 Derivation of the Near-Field Enhanced Scattering
In this section, the derivation of the coupled-oscillator model of the hybrid Au-LiNbO 3 nanostructures is developed from first principles. Section S.2.4.1 details the calculations of the relevant field quantities and Section S.2.4.2 translates these field quantities into an osicllator model.

S.2.4.1 Nonlinear Optics First Principles
The oscillator model of the modes of the LiNbO 3 microsphere is generated from the solutions to the coupled nonlinear wave equation for the total electric field wherein J 0 (r, ω) is the current density of the laser. The polarization fields P (1) (r, ω) = χ 2 (r)E(r, ω) and build in the first-and second-order material response of the LiNbO 3 sphere, respectively. The sphere's electric susceptibilities are allowed to vary in space such that χ 2 (ω , ω − ω )Θ(r ≤ a 2 ) with 1 3 the 3 × 3 identity matrix. The sphere's dielectric function is then given by 2 (r) = 1 + 4πχ (1) 2 (r). From here, the electric field can be perturbatively expanded with the assumption that |χ (2) 2 (r; ω , ω−ω )| is a small quantity. Explicitly, we let E(r, ω) = E 0 (r, ω) sca (r, ω), wherein the first term is the laser's electric field and the terms in the sums contribute n th -order corrections to the vacuum-like and scattered fields E vac (r, ω) = ∞ n=2 E (n) 0 (r, ω) and E sca (r, ω) = ∞ n=1 E (n) sca (r, ω), respectively, that are set up by the polarized sphere. Subtracting the laser field equation ∇ × ∇ × E 0 (r, ω) − (ω 2 /c 2 )E 0 (r, ω) = (4πiω/c 2 )J 0 (r, ω) from Eq. (S35), one finds The currents in right-hand-sides of Eq. (S37) are bound currents that are only nonzero where the first-and second-order susceptibilities of the LiNbO 3 sphere are nonzero, respectively, such that build a model that includes only these two sources and temporarily neglects interactions between the NPs and microsphere. Ignoring the effects of radiation, the scattered potential of the NP is (S42) for χ 1 (r, ω) = χ 1 (ω)Θ(r ∈ V 1 ) and r not on the boundary of the NP's volume V 1 . Assuming the laser wavelength λ 0 = 2πc/ω 0 is much longer than the extent of the NP, the incident potential can be written as Φ 0 (r, ω) ≈ r · lim k→0 E 0 (r, ω) Assuming that the solutions to Eq. (S42) are near zero at ±ω 0 , i.e. the NP does not respond strongly near the fundamental frequency of the laser, the term proportional to Φ 0 (r, ω) on the RHS can be dropped. Thus, with the dielectric function 1 (r, ω) = 1 (ω)Θ(r ≤ a 1 ) + 1Θ(r ≥ a 1 ) defined such that the NP lies at the origin, one finds: wherein the (scalar) Green's function G NP (r, r , ω) is the standard solution to the Poisson equation in spherical coordinates 3 . Much like the dyadic Green's function for the LiNbO 3 fields, G NP (r, r , ω) = G sca (r, r , ω) + G 0 (r, r , ω) is separable into a scattering part and a vacuum-like or free-space part. The former encodes the resonances of the NP, while the latter serves simply to satisfy the principle of superposition and is otherwise unimportant.

S.2.4.2 Construction of the Oscillator Model
The above definitions are clarified by the forms of both G LNO and G NP , which are cumbersome but manageable. The former, in the case where both the source charges and the observer are outside the sphere's surface, is given by for a sphere centered at the origin. Each of the mode functions (a +1 1 /r +1 )P m (cos θ)S p (mφ) describes the spatial variation of the observables of an electric multipole mode with a characteristic response function 14 that describes its oscillations in time. Here, and m give the order and degree of each mode's corresponding spherical harmonic and p the reflection symmetry of each mode across the x-axis.
The Green's function of the LiNbO 3 sphere can be similarly expanded using a set of socalled quasinormal geometric resonances of both magnetic and electric multipole symmetry 2 with response functions that depend on the same spherical harmonic symmetry parameters p, , m. For an observer inside the sphere, this Green's function is: where again the sphere is centered at the origin. The mode functions M p m and N p m are identical to the regularized vector spherical harmonics given in (S16) but with the spherical Hankel functions replaced with spherical Bessel functions, h (1) (x) → j (x). The Green's function for an observer outside the sphere is, similarly: with the forms of the unwieldy response functions A p m (ω), B p m (ω), C p m (ω), and D p m (ω) described in Ref. 13. Further, G 0 (r, r ; ω) is the Green's function of free space.
To compare the responses of the NP and microsphere, we can see that B < p m (ω) obeys the simple relation lim k→0 ( √ 2 ) B p m (ω) − 1 = (1 − 2 )/( 2 + + 1). This gives us a foundation from which to make quantitative comparisons between the oscillator parameters assigned to each response in either particle. In detail, using the Drude-Lorentz dielectric model of Au from Section S.1.2.2, we can see that the dipole response function g 1 (ω) describes the oscillations of three coupled modes, one which has primarily free-electron character (the dipole plasmon) two others which are mostly comprised of the material's interband resonances. Explicitly, wherein the mass µ, resonance frequency ω, complex eigenvalue Ω, and phase offset ψ of the plasmon are labeled with subscript 1, and the oscillator parameters of the fictitious oscillator are labeled with subscripts L i . Explicit values of these parameters are given in Section S.2.2.
As is also shown in Section S.2.2, there are many modes of the LiNbO 3 sphere that have significant response magnitudes in the range of laser energies in which SHG is observed (∼2.3-2.4 eV, see Figure 4c). To simplify their description, we use an expansion of the LiNbO 3 mode response functions: is an overlap integral over the radial components j (kr) of the fundamental mode functions. Further, the functions I α 1 α 2 α 3 are overlap triple integrals wherein X T p m (r, k) = M p m (r, k)δ T,M + N p m (r, k)δ T,E .
Using the definition of the dipole polarizabilities α ν (ω) given in Section S.2.2.2, we can define the motion of each component of a single Au dipole and each Mie resonance in the absence of their mutual coupling as Finally, to reproduce the equations of motion of Eq. (1), we must define the external forces acting on the particle moments and introduce the coupling forces to Eq. (S55). The former can be quickly written as F 1ν (ω) = eê ν · E (2) 0 (r 0 , ω) and F 2β (ω) = eE 0 πδ(ω − 2ω 0 )C β (ω). The latter arise from the interaction energy one is immediately delivered the equations of motion in the main text with σ βν = (e 2 /a 3 2 )ê ν · X β (r 0 , k β ). A detailed description of the solutions to the equations of motion are given in Section S.1, and an analysis of the contribution of the various modes of the LiNbO 3 sphere to enhancement signal is shown in Figure S14.
In closing, we note that, as the solutions to the equations of motion discussed in the main text and highlighted in Section S.1 involve a second perturbation expansion of the equations of motion of each Mie resonance, two sets of superscripts arise that correspond to two separate expansions. In order to clarify the notation, we can see that no terms in the solution to the wave equation are kept beyond second order, such that we can replace superscripted variable names with script names. Further, we can drop the superscripts of the first-order terms altogether, such that, elsewhere in the SI and main text, we let E Figure S1: Scanning transmission electron microscopy (STEM) images. LiNbO 3 particle images in (a) bright-field and (b) dark-field modes.    nanostructures. This analysis included measurements obtained from 125 independent Au NPs. The variance of 1.5 nm is reported as one standard deviation from the calculated mean of 10 nm. Figure S6: Au NP absorbance. Ultraviolet (UV)-visible absorbance spectrum of ∼10-nm diameter Au NPs suspended in an aqueous solution. This spectrum indicates the plasmonic band for the nanoparticles was centered at ∼520 nm, which is close to the 530 nm plasmonic band for hybrid Au-LiNbO 3 hybrid particles. Figure S7: Minimal SHG from Au NPs. The SHG analyses for pristine ∼10-nm diameter gold nanoparticles recorded by sweeping the excitation laser wavelength from 850 to 1,080 nm with a step size of 5 nm. These results indicate a lack of SHG response in the centrosymmetric Au NPs. Figure S8: Images of low NP-loading hybrid nanostructures. Hybrid Au-LiNbO 3 particles prepared using 0.1 mL of a 5 mM HAuCl 4 solution to perform the in situ synthesis of the Au NPs. The resulting products were characterized by scanning transmission electron microscopy (STEM) operating in: (a), (c) a high-angle annular dark-field (HAADF) mode; and (b), (d) a bright-field mode. Figure S9: Elemental analysis of low NP-loading hybrid nanostructures. Energy dispersive X-ray spectroscopy (EDS) analysis of the hybrid Au-LiNbO 3 particles prepared using 0.1 mL of 5 mM HAuCl 4 . These images show (a) a HAADF image obtained by STEM techniques, (b) an overlay of the Au NPs (as detected by EDS) on the HAADF image of the assemblies, (c) an EDS map of the Au NPs, and (d) an EDS map of Nb within these assemblies. Figure S10: SHG ehancement spectra comparison. Evaluation of the enhancement in the SHG signal as a function of the fundamental wavelength (i.e., ranging from 850 to 1,070 nm) of the incident laser for the individual hybrid Au-LiNbO 3 particles prepared with a lower loading of Au NPs. After normalization of the SHG output of the individual Au-LiNbO 3 hybrids to the SHG output of the bare LiNbO 3 , enhancement factors were plotted against the FWs of the measurement. Figure S11: LiNbO 3 second-order susceptibility fit. Comparison of the real and imaginary parts of the magnitude χ (2) (ω, ω) ≡n·χ (2) 2 (ω, ω)·n of the second-order susceptibility in the limit where the output frequency is twice the input frequency (n is any real unit vector).    The grids were arranged as in (a) but with each NP randomly shifted between -a 1 and a 1 in the in-plane directions and up to ±2a 1 along the out-of-plane axis such that no two particles overlap. The observed probability density (blue) is fit to an extreme value distribution (black) to estimate the proportion of the particles with cross sections lower and higher than the average from the ordered grid of 16 particles (red). (c) Demonstration that the absorption cross sections of the Au NPs approach stable values as the number of particles in the ensemble increases, both for parallel (blue) and perpendicular (red) incident light. Points show ensemble averages of the per-particle cross section maxima for ensembles of size N = 1 to 225, while the solid lines show the extracted averages of 15.7 nm 2 for parallel light and 10.5 nm 2 for perpendicular. The averages are drawn over the range of data points from which they are calculated. (d) Probability density of the spectral position of the scattering maxima of the NPs from 100 simulations of N = 16 NP ensembles. The energy of the scattering maximum of a single particle (see Figure S12a) and the average of the energies of the scattering maxima from a simulation of an N = 225 NP ensemble are shown for reference. Inset: A zoomed-in picture of the probability density, with identical quantities and units plotted on the axes as the main figure.  Figure S18: Input and SHG lineshapes. Normalized lineshapes of the input laser (red) and SHG (blue) power spectra.