Extremely large third-order nonlinear optical effects caused by electron transport in quantum plasmonic metasurfaces with subnanometer gaps

In this study, a third-order nonlinear optical responses in quantum plasmonic metasurfaces composed of metallic nano-objects with subnanometer gaps were investigated using time-dependent density functional theory, a fully quantum mechanical approach. At gap distances of ≥ 0.6 nm, the third-order nonlinearities monotonically increased as the gap distance decreased, owing to enhancement of the induced charge densities at the gaps between nano-objects. Particularly, when the third harmonic generation overlapped with the plasmon resonance, a large third-order nonlinearity was achieved. At smaller gap distances down to 0.1 nm, we observed the appearance of extremely large third-order nonlinearity without the assistance of the plasmon resonance. At a gap distance of 0.1 nm, the observed third-order nonlinearity was approximately three orders of magnitude larger than that seen at longer gap distances. The extremely large third-order nonlinearities were found to originate from electron transport by quantum tunneling and/or overbarrier currents through the subnanometer gaps.

In addition to the enhanced electromagnetic fields in plasmonic systems with subnanometer gaps, they also produce accompanying quantum mechanical effects. In the linear response of an isolated nanodimer system composed of two metallic nanoparticles with a small gap, quantum effects have been shown to affect the optical characteristics as previously reported in theoretical [26][27][28][29][30] and experimental [31][32][33] studies. The effects depend strongly on the relationship between the Fermi energy and the potential barrier at the gap. When the nanodimer is separated by a large gap distance, the Fermi energy is sufficiently lower than the potential barrier. At smaller gap distances comparable to the radius of the constituent nanoparticles, the potential barrier starts to decrease; however, the barrier height is still quite high compared to the Fermi energy. When the gap becomes smaller than 0.4 nm, the Fermi energy is approximately equal to or slightly greater than the potential barrier. In this case, electrons may cross the barrier via quantum mechanical tunneling and/or overbarrier currents through the gap. These currents produce charge transfer through the gap that suppresses the plasmonic enhancement.
In nonlinear responses, the quantum current flowing in the nanodimer has been reported to affect the electric field enhancement 34 and the harmonic generation efficiency 35 . However, to the best of our knowledge, there have been no prior theoretical or experimental reports that discuss how such currents flowing across the gaps contribute to the nonlinearity of the metasurface with subnanometer gaps. Although there have been a few recent measurements, the gap distance was 0.6 nm in the smallest case, where the nonlinearities of the metasurface should still be solely determined by the strong optical enhancement at the gaps 25 . Therefore, it is highly intriguing to explore whether higher nonlinearity can be achieved by the currents flowing through the gaps of the metasurface.
In this study, quantum plasmonic metasurfaces composed of metallic nanospheres with subnanometer gaps were theoretically investigated through a fully quantum mechanical calculation using time-dependent density functional theory (TDDFT) 36,37 combined with a two-dimensional (2D) coarse-graining approach 38 to the electromagnetism. We investigated their third-order nonlinearities for off-resonant incident pulses and clarified their dependence on the fundamental frequency and the gap distance. Our results show that, for gap distances of ≥ 0.6 nm, third-order nonlinearities monotonically increase as the gap distance decreases owing to the optical enhancement at the gap. In particular, when the frequency of the third harmonic generation overlaps with the plasmon resonance frequency, strong resonant third-order nonlinearities were observed. At smaller gap distances down to 0.1 nm, we found the appearance of extremely large third-order nonlinearity that is not assisted by the plasmon resonance. At the gap distance of 0.1 nm, the observed third-order nonlinearity is approximately three orders of magnitude larger than those at longer gap distances. It was found that the extremely large third-order nonlinearities originate from electron transport by quantum tunneling and/or overbarrier currents through the subnanometer gaps. Our findings demonstrate a new way to increase the nonlinearities of metasurfaces that should be enormously useful for downsizing all-optical switches.

Results
Studied system and theoretical approach. Figure 1a displays the studied system where metallic nanospheres with a diameter a are periodically arrayed in the xy plane with a gap distance d and a period length l. The incident light of a planar pulse propagates toward the negative z direction with the x-electric and y-magnetic components, E i and H i , respectively. The time profile of E i is described in the Supplementary Information. To treat the quantum mechanical effects with a moderate computational cost, we employed the jellium model (JM) in which ionic structures are replaced by a positive background charge of a spherical shape with a shape boundary. Although this JM includes a considerable simplification, the plasmonic motion of electrons in the nanoparticles are well described 26,28,39 . Previous studies reported quantitative agreement between the JM and measurements of plasmonic systems with subnanometer gaps where the quantum electron tunneling played a key role 31,32 . In the JM, the medium is specified by the Wigner-Seitz radius, r s , that specifies the average charge density, n + = ((4π)r s 3 /3) −1 . We employed a value of r s = 4.01 Bohr corresponding to the Na metal. We set a = 3.1 nm, where each nanosphere included 398 electrons that constitute the closed shell structure. This size is sufficiently large to ensure that the nanoparticle exhibits a well-developed plasmonic resonance 29 . www.nature.com/scientificreports/ To calculate the optical responses of the metasurface, we employed the TDDFT that has been extensively used to investigate the optical properties of molecules 40 and solids 41 at a first-principles level. We combined the TDDFT with a 2D coarse-graining approach in which the light-matter interaction in two-dimensional materials is aptly described by coupling to the Maxwell's equations 38 . Adiabatic and local density approximations were used for the exchange-correlation potential 42 . All the calculations were carried out using SALMON, an open-source code (https ://salmo n-tddft .jp/) developed in our group 43 . The Supplementary Information contains a detailed description of the adopted numerical approach.
Linear optical response. Before discussing nonlinear optical responses, we briefly look back on the linear response of the metasurface. Figure 1b shows the linear optical absorption rate of the metasurfaces as the gap distance d was varied. When d is sufficiently large, the optical absorption shown by a bold red and yellow band appears at approximately 3 eV, which is not significantly different to the plasmon resonance of a single nanosphere. As the gap distance d decreases to 0.2 nm, the frequency of the absorption starts to be red-shifted, and the magnitude increases. These features originate from increasing interactions between the nanospheres. The plasmonic charge densities induced on the spheres are strongly attracted to each other, forming the bonding dipolar plasmon mode. However, at the locus of d ≤ 0.2 nm, the plasmon resonance rapidly decays and hybridizes into multiple plasmon modes including the bonding octopolar-and the void-plasmon modes that are caused by quantum tunneling and/or overbarrier currents flowing through the subnanometer gaps. At the locus of d ≤ 0 nm where the constituent spheres start to overlap geometrically, there is no potential barriers that prevent a conduction current from flowing throughout the metasurface. These trends had already been established in our previous study 44 . Nonlinear optical response. Next, we explored the nonlinearities in the optical response of the plasmonic metasurfaces. First, we focused on the metasurface with d = 2 nm where the nanospheres are sufficiently separated from each other, and the response reflects the characteristic features of a single nanosphere, as seen in Fig. 1b. In the 2D coarse-graining approach adopted here, we employed a macroscopic description in which the metasurface shown in Fig. 1a was treated as a uniform thin-film of zero-thickness 38 . In this approach, the evolution of the electric field was described by where E i and E t are the incident and the macroscopic transmitted electric fields, and J[E t ] is the 2D macroscopic electric current density that includes nonlinear signals and is produced by the field E t . The current J[E t ] was calculated from the microscopic electron dynamics for which we utilize the TDDFT. The boundary condition on the film, E t = E i + E r , determines the macroscopic reflected electric field, E r . The details of the theory are explained in the Supplementary Information. These macroscopic reflected and transmitted electric fields, E r and E t , should be observed in actual measurements. Figure 2a displays the time profiles of E i and E r . In all subsequent results, electric fields divided by the maximum amplitude of E i will be shown and denoted as E . The full duration of the incident pulse was set to 55 fs with the envelope shaped by a cosine-squared function. The mathematical expression for E i is described in the Supplementary Information. The top panel displaying the black solid line is E i , where the fundamental frequency ω i is set to 0.96 eV, far from the plasmon resonance at d = 2 nm, which is approximately 3 eV as seen in Fig. 1b. The middle panel shows E r calculated for the three different incident intensities I = 10 11 , 10 10 , and 10 9 W/cm 2 that are plotted by the blue, red, and green lines, respectively. Since we applied the off-resonant incident pulse, E r was much smaller than E i , indicating a high transparency. After t ≈ 20 fs, small nonlinear signals that depend on I were observed, as seen in the magnified box. To distinguish between the nonlinear signals more clearly, we calculated, E r , the difference between E r and the lowest intensity reflection of 10 9 W/cm 2 . The bottom panel shows �E r (t) for the two cases of I = 10 11 and 10 10 W/cm 2 , where the third harmonic generations ware clearly detected as the intensity increased. In all subsequent results, we used a single intensity of I = 10 10 W/cm 2 that corresponds to the intensity of the pulse used in the previous experimental study on the plasmonic metasurface with subnanometer gaps 25 . To clarify the dependence on the fundamental frequency ω i , we examined three different ω i values for the same metasurface, characterized by d = 2 nm. Figure 2b shows the resultant power spectrum of the normalized reflected electric field | E r (ω)| 2 , where the plasmon resonant frequency ω r is indicated by the vertically drawn gray line.
The blue, red, and green lines correspond to ω i = 0.82 (= ω r /3.5), 0.96 (= ω r /3), and 1.15 (= ω r /2.5) eV, respectively. At the first harmonic generation, | E r (ω)| 2 slightly increased with ω i because the frequency comes slightly closer to the plasmon resonance ω r . In contrast, for the third harmonic generation appearing around 3 eV, the highest nonlinearity was achieved at ω i = 0.96 eV. This is because the third-order signal appears closely to ω r and thus is plasmonically enhanced. There is a time-delay in the third-order harmonic generation with respect to the incident and the reflected fields, as seen by comparing the top, middle, and bottom panels of Fig. 2a, due to the inherent time requirement for resonant enhancement. The enhanced nonlinearity assisted by the plasmon resonance was also reported in the previous study that dealt with an isolated nanodimer theoretically 31 . Finally, we note that the fifth harmonic generation is visible only at ω i = 0.96 eV.
Dependence on gap-size. We now move on to the main subject of the present study, clarifying the thirdorder nonlinearities of metasurfaces with various subnanometer gaps. To quantify the third-order nonlinear efficiency, we introduce a quantity, R , the detailed definition of which is described in the Supplementary Information. In simple terms, R NL indicates the nonlinear reflectivity caused by third harmonic generation for the case with a gap distance d and the incident pulse with a fundamental frequency ω i . Figure 3a summarizes the resultant R NL for gap distances from − 0.2 to 2 nm and the fundamental frequencies ω r /3.5 < ω i < ω r /2.5, where ω r is the plasmon resonant frequency. When the fundamental fre-

Scientific Reports
| (2020) 10:21270 | https://doi.org/10.1038/s41598-020-77909-y www.nature.com/scientificreports/ quency satisfies ω i = ω r /3, it is marked in Fig. 3 by a black pentacle. At d = 0 nm, ω i is widely sampled from 0.47 to 1.1 eV, and the pentacle is not added to this case because the plasmon resonance is hardly distinguished, as seen in Fig. 1b. Figure 3a allows us to find three distinctive trends that are indicated by the different colored lines. The first trend is indicated by the blue lines representing d = 2, 1.2, and 0.6 nm. Here, the peaks appear close to the pentacles. The peak values gradually increase and the peak frequencies are red-shifted as the gap size decreases. These trends are in accordance with the linear responses shown in Fig. 1b, indicating that the large third-order nonlinearities are assisted by the plasmonic resonance of the bonding dipolar mode. The second    NL ] shows a prominent maximum at d = 0.1 nm. To quantify the observed extremely large third-order nonlinearity of the studied metasurfaces, they were compared to the nonlinearity of a SiO 2 thin film which is conventionally used in all-optical switches 25 . As described in the Supplementary Information, we estimated the R NL of a SiO 2 thin film with the same thickness, a, as the present metasurface. The value was calculated to be 1.43 × 10 −13 , eight orders of magnitude lower than the studied metasurface.
The origin of high nonlinearity. To clarify the physical origin of the different features observed for gap distances above and below d = 0.6 nm in Fig. 3b, we examined the nonlinear electric current density flowing in the metasurface. In the 2D coarse-graining approach employed here, the reflected field E r (t) is equivalent to the 2D macroscopic electric current density J(t) up to a constant factor, E r = −(2π/c) J[E t ] , and J(t) is given from the microscopic current density j x, y, z, t by J(t) = dxdy/l 2 dzj x, y, z, t . To investigate the physical mechanism that produces the nonlinear behavior, we introduce the nonlinear current density defined as To elucidate the behavior causing the extremely large third-order nonlinearities, we focused on the temporal and spatial distribution of the nonlinear electric current density. The upper panels of Fig. 4b-e illustrate the spatio-frequency distribution of the nonlinear microscopic electric current density j NL (x, ω) 2 . The fundamental frequency, ω i , is set to the value that gives max[ J NL ] for the metasurfaces with d = 2, 0.4, 0.1, and − 0.2 nm. The vertical axis is the x-axis of the metasurface that is parallel to the polarization direction of the incident pulse. The horizontally drawn pink dashed lines mark off the periodic length l with the left schematics showing the constituent nanospheres. The lower panels indicate the nonlinear 2D macroscopic electric current density J NL (t) . In Fig. 4b, for d = 2 nm, J NL (t) solely consists of the plasmonically assisted third harmonic generation whose spatial distribution is fully confined to the sphere, as seen from the distribution of j NL (x, ω) 2 . This is explained by the position of the Fermi energy which is much lower than the potential barriers at the gaps. Although small second harmonic components are visible in j NL around the edge of spheres, they vanish after the spatial integration  Figure 4c displays the case of d = 0.4 nm, where the potential barrier is slightly higher than the Fermi energy. The nonlinear current is still mainly composed of the third harmonic component. The amplitude is larger than that of the d = 2 nm case owing to the stronger electric field enhancement at the gaps. It was noted that a slight nonlinear component of the fundamental frequency is visible in J NL at 10 ≤ t ≤ 30 fs. We confirmed that this weak signal is also seen in j NL (x, ω) 2 flowing across the gaps by quantum mechanical tunneling, although it cannot be seen on the color scale in the upper panel. As seen in Fig. 4d, J NL (t) at d = 0.1 nm shows drastic changes compared to the above-mentioned trends at larger gap distances. At this distance, the Fermi energy is slightly higher than the potential barriers, changing the mechanism of the current from tunneling to overbarrier. The nonlinear current includes both n = 1 and n = 3 components and propagates through the subnanometer gaps. As seen in Figs. 3b and 4a, this gap distance produced the largest third-order nonlinearity. At d = − 0.2 nm shown in Fig. 4e, the third harmonic generation current rapidly decreases. Here, the tunneling and/or the overbarrier currents are unified into a conduction one  Fig. 3b. In the upper panels, the vertical axes denote the x-axis of the metasurface which is parallel to the polarization direction of the incident pulse. The horizontally drawn pink dashed lines mark off the periodic length l along the x direction.
Scientific Reports | (2020) 10:21270 | https://doi.org/10.1038/s41598-020-77909-y www.nature.com/scientificreports/ flowing throughout the metasurface because nanospheres are directly connected. We noted that the magnitude of the nonlinear current decreases compared to the case of d = 0.1 nm. From these observations, we conclude that the extremely large nonlinearity is caused by electron transport through the gaps via the tunneling and/or the overbarrier mechanisms. While performing this study, we used various approximations. They are outlined below in conjunction with the limitations of this study. Since we assume the Na nanospheres described by the JM, it does not include any d-electron effects that appear in typical plasmonic materials such as the noble metals. Therefore, in our simulation, the nonlinearity caused by the d-electrons is ignored. Moreover, the JM ignores the ionic structure of metallic nanoparticles that may cause strong electric-field enhancements at the apexes of clusters, known as the lightning rod effect 29 .
Although we have focused on very small nanospheres with a diameter of a = 3.1 nm, actual plasmonic nanoparticles used widely in measurements are larger in size, reaching 10-100 nm 1, 2 . Furthermore, despite the great advantage that the optical characteristics of metasurfaces can be finely tuned over a wide range by the shape of the nanoparticles, this study remained limited to the elementary geometry of spheres. These geometrical differences would affect electron tunneling and/or overbarrier currents, modifying the nonlinearities. To explore such effects, one prospective candidate method is quantum hydrodynamic theory (QHT) that can describe the nonlinear light-matter interaction of plasmonic systems [45][46][47] . In particular, a recently proposed QHT study has revealed that it can be directly derived from the TDDFT with the JM, demonstrating good agreement in the linear response regime 47 . Nevertheless, the computational cost of the QHT is significantly lower than the TDDFT because the QHT is an orbital free approach. Therefore, the QHT is expected to allow investigations of much larger systems with various geometries. We consider that the application of the QHT to metallic metasurfaces and the comparison with the TDDFT results in the nonlinear regime is an important topic for a future study.

Discussion
In conclusion, we have presented a theoretical investigation of third-order nonlinearities in quantum plasmonic metasurfaces with subnanometer gaps, mainly focused on third harmonic generation. Nonlinear optical responses of a metasurface composed of metallic nanospheres have been examined using the TDDFT with the JM, where an off-resonant incident pulse ensuring negligibly small cumulative thermo-optical effects was employed. We have calculated the third-order nonlinearities as functions of the fundamental frequency of the pulse and the gap distance of the metasurface. It has been shown that, for gap distances of ≥ 0.6 nm, the thirdorder nonlinearities monotonically increase as the gap distance decreases. This is caused by enhancement of the screening charge densities that are induced at the interfaces of the nanospheres. When the frequency of the third harmonic generation overlaps with the linear plasmon resonance, large third-order nonlinearities are observed. At further smaller gap distances down to 0.1 nm, we find the appearance of the extremely large thirdorder nonlinearity that is not assisted by the plasmon resonance. In particular, at the gap distance of 0.1 nm, we have achieved a third-order nonlinearity three orders of magnitude larger than that of the longest gap distance case. The extremely large third-order nonlinearities originate from electron transport by quantum tunneling and/or overbarrier currents through the subnanometer gaps. At gap distances of d ≤ 0 nm where the spheres geometrically overlap, the tunneling and/or the overbarrier currents are unified into a usual conduction current that flows throughout the metasurface. In that case, the third-order nonlinearities were observed to decrease. Our findings suggest a new way to increase the nonlinearity of metasurfaces, which is expected to be enormously useful for downsizing all-optical switches.