Metallic nanostructures as electronic billiards for nonlinear terahertz photonics

Optical properties of metallic nanoparticles are most often considered in terms of plasmons, the coupled states of light and quasi-free electrons. Here we predict that confinement of electrons inside the nanostructure leads to another, very different, type of resonance, which determines the optical properties in the frequency range significantly below the plasmonic resonance. We demonstrate that closely placed confinement-induces resonances typically join into a single composite"super-resonance"which produces giant nonlinearity at low frequencies. Our simulations show how such nonlinearities can be used for efficient down-conversion of optical pump to terahertz and mid-infrared frequencies in sub-micrometer devices based on nanoparticle composites. We discuss the interaction of these quantum-confinement-induced resonances with the conventional plasmonic ones, as well as the unusual quantum level statistics, adapting here the paradigms of the electronic billiard theory.

Introduction.Light propagating in the vicinity of metallic surfaces or metallic nanostructures is strongly coupled to the "electronic fluid" formed by quasi-free electrons in the metal, resulting in joint electron-photon excitations which are called surface or particle plasmons, depending on the geometry.Plasmons, and the corresponding plasmonic resonance (PR) are at the very heart of optics of nanostructures.
PRs appear by matching of the incoming light to the intra-particle fields, leading to strong surface charges and resonance peaks of linear and nonlinear response of metallic nanoparticles at certain frequencies [1][2][3].That is, PRs are defined via the matching condition of the fields, rather than electrons themselves, and have no direct relation to the electron confinement inside the nanostructure.Plasmonic resonances are located at quite high frequencies, commonly in the visible range.
For very small nanoclusters and nanoparticles (few hundreds of atoms and below), the response is molecularlike [9,14], typically including a cacophony of resonances replaced by a single PR [14,28,29].At lower frequencies, a prominent molecular-like resonance is the HOMO-LUMO transition [14,[30][31][32].It describes an excitation of a single localized electron, and depends heavily on the molecular structure [30].
As we shift to larger nanoparticles, the HOMO-LUMO transition effectively fades out [14] because its decay rate grows exponentially with the size of the nanostructure [31,32].However, the PR is not the only one which remains.Signatures of a single resonance well below the plasmonic one but above the HOMO-LUMO transition were observed experimentally [33,34] and later confirmed theoretically [8,35] involving ab-initio simulations [35].This resonance is however mostly overlooked since then, and there is no consensus in explanations of its nature.Whereas in [8,34] the confinement-based argument were put forward, [33,35] tries to explain it without leaving the plasmonic framework by introducing "restoring force on the electrons".
Here we show that confinement-based resonances under normal conditions (at least for nanostructures with regular enough geometry) merge together and form a sin-FIG.1.
A spherical gold nanostructure (left) of 2.75 nm diameter and the band level structure (right) appearing due to the electron confinement (quantum billiard resonances).Fermi energy EF and work function EW as well as exemplary virtual transitions contributing to linear (orange arrows) and nonlinear (red arrows) properties are indicated.One of the wavefunctions ψ(x, y, z) is visualized inside the nanostructure (a sector of the sphere is cut for visualization purposes).Even for the simplest nanoparticles, energy level structure is quite complex.
gle broadband "super-resonance" located commonly in THz and MIR frequency range, with both the width and position widely controllable with the nanostructure geometry.Such super-resonance must be considered as one of the universal signatures of metallic nanostructures, yet fundamentally different from the PR.
We study the statistical features of confinement-based resonances by adapting and modifying the paradigm of neighboring-level statistics from the field of quantum billiards [36].Electrons confined in a nanostructure certainly represents a type of quantum billiard.Yet, up to now, with only very few exceptions [37,38], nonmetallic billiards such as semiconductor quantum dots [39][40][41][42][43][44][45] were considered.As we show, the metallic nature of our billiard provides an unique opportunity to observe certain features the level statistics directly in the optical properties.
Furthermore, the super-resonance provides a broadband nonlinear response, leading to giant nonlinearities.We demonstrate how these nonlinearities can be used for an efficient optical rectification and differencefrequency mixing in the nanocomposites, enabling broadband conversion from optical to MIR or THz ranges by sub-micrometer devices.
The model.We used a simple analytical approach of a single-particle-in-a-box [6,8,16,34] (see also more details in Supplementary) with electron fully confined inside a nanoparticle.For calculation of optical properties, we consider electron energy structure characterized by Fermi energy E F and work function of E W , as illustrated in Fig. 1 right.This approximation works well for metals with a simple Fermi surface such as alkali metals (Li, Na, Ca, Rb), it is an acceptable simplification for metals with somewhat more complicated Fermi surfaces such as Cs, Cu, Ag or Au, and is barely applicable at all for other metals.The linear (χ (1) ) and nonlinear (χ (3) ) susceptibilities were calculated using a version of the standard perturbative iterative approach [46] which takes into account selection rules following from the Fermi-Dirac statistics (see Supplementary).The linear and nonlinear optical properties can be described via a sum of contributions of virtual transitions from inside the Fermi sea to outside and back [see orange (linear) and red (nonlinear) arrows Fig. 1].We also assumed fast population decay time T 1 = 50 fs and dephasing time T 2 = 5 fs [16].
To be more specific, in our numerical simulations we consider gold since it is a very widespread material, and is suitable for composites due to low imaginary part of susceptibility.Yet we note that conclusions we draw below are basically metal-independent (taking into account precautions mentioned above).In gold, a simple ideal-metal picture discussed above neglects several linear and nonlinear effects, such as interband transitions, influence of finite temperature, and hot-electron nonlinearities.However, these effects are negligible for lowfrequency response in THz or MIR range driven by femtosecond pulses, and our model remains adequate in this regime (see Supplementary for justifications and detailed estimates).
Quantum billiard (confinement-based) resonances.The typical level structure obtained by the above model for an exemplary spherical gold nanoparticle of the diameter d = 2.75 nm is shown in Fig. 1.These levels originate from electron confinement in the nanostructure.Even in the presented case of a very simple particle, the levels look quite irregular.This is a familiar picture in the framework of of quantum billiard theory [36], were the statistical properties of the level distribution play one of the central roles.For instance, one can consider the neighboring level statistics (NLS), which allows to distinguish between integrable (regular) and non-integrable (chaotic) billiards.For the regular billiards, such as spheres, the probability density P (ω) of neighboring-level distance ω obeys Poissonian statistics: ln P (ω) ∝ −ω.This is also true in our case: NLS corresponding to Fig. 1 is plotted in Fig. 5 Judging from such statistics, one might expect a conglomerate of resonances near zero frequency, but this is not the case.An optical response χ (1) resulting from the electron confinement for few exemplary nanostructures is shown in Fig. 5.Note that, in addition to the confinement-based impact shown in Fig. 5, the full linear response includes also the so-called Drude part, representing the action of quasi-free-electrons (see Supplementary for more details).
The clearly observed feature of the confinement-based linear response is the presence of a single resonance-like peak in the IR/THz range at a non-zero resonance frequency ω conf which quickly decreases with increasing particle size.Such peak was observed experimentally [33,34] and theoretically [8,35].It is easy to see that this resonance coincides well with the minimally possible allowed transitions close to the Fermi energy which can be, for a spherical nanoparticle, analytically estimated as (see Supplementary for details), cf. also [34]): where m e is the electron mass.This makes it fundamentally different from the PR, being located at a significantly lower frequency.Yet, why does only a single resonance arise?Can we influence it width and position?These and related questions will be addressed in the following paragraphs.Closer consideration allows to establish that, starting already from quite small nanosphere diameters d ≈ 2.5 nm, the resonance near ω conf is composed of many transitions with nearby frequencies [see Fig. 5(a)].These transitions merge into one single "super-resonance".This is especially well observable if we consider much larger T 2 (which would correspond to low temperatures [47]).In this case, many separated resonances are indeed visible in the optical response, as shown in Fig. 2(a,b) by black curves.The particular structure of the transitions depends significantly on the geometry [cf.Fig. 5(b) for a cylinder].For instance, the position and width of the super-resonance for a cylinder are shifted in comparison to a sphere with the same volume and diameter by the noticeable amount of 35% and 23%, correspondingly.Nevertheless, Eq. 1 remains a valid, yet rough estimation of the position of the resonance.
Both the position of the super-resonance ω conf and its structure can be analyzed using the level-distance statistics similar to NLS.A naïve approach would be to calculate such statistics using all dipole-allowed transitions between the confinement-based levels, shown in downward green rectangles in Fig. 2(a), and covering extremely broad range around 10 eV (see green bars in Fig. 2(a) and also Supplementary).However, it must be modified to include only transitions from below to above of E F , obviously corresponding to the Pauli principle and absence of population above the Fermi level.In addition, in the statistics we weight the transitions by the square of the corresponding dipole momentum, thus taking into account the known tendency of the transition dipole momenta to rapidly decrease, on average, with the energy difference.The resulting modified statistics is shown by the red and blue bars in Fig. 5 and agrees nicely (for large T 2 ) both with the position of the super-resonance and with its width for small T 2 .The super-resonance has a certain "natural" width, for the case of nanospheres it can be estimated as ω conf /4 (see Supplementary).Analysis of the position of the super-resonance for nanospheres (see Supplementary) indicates that the energy of the participating states is located mostly in radial (rather than angular) motion.
Nonlinearities.The above described low-frequency resonance is expected to lead also to strong nonlinearities; in our case, χ (3) = 0 as calculated using the approach described above.We note that such approach to calculate Kerr nonlinearity was already utilized in [16], however, instead of discrete spectrum, approximation of continuous density of states was used.Nevertheless, we checked that our calculations are in quantitative agreement with [16]; they are also in agreement with experimental measurements for short pulses (see [11] and references therein).An example of χ (3) (ω; ω, ω, −ω) for the four-wave-mixing (FWM) process ω + ω − ω = ω (corre- sponding to the Kerr nonlinearity) is shown in Fig. 3(a) for several diameters.The low-frequency resonance we observed in χ (1) is also well-visible here.Whereas in the linear response the Drude part dominates (see Supplementary), in the nonlinear response it is fully absent.
We now try to exploit this low-frequency resonance.Motivated by detection and spectroscopic applications of THz and MIR radiation, we focus on the FWM providing a signal in THz and MIR range, generated from a sub-100 fs pump pulse.Nonlinear susceptibilities χ (3) (δ; ω, ω, −2ω + δ) and χ (3) (δ; ω, ω + δ, −2ω), leading to generation of low-frequency signal at frequency δ as a result of a FWM process in a two-color pump at frequencies ω (fundamental) and 2ω (its second harmonics) are presented in Fig. 3(b-c).Both Kerr rectification nonlinearities presented in Fig. 3 are several orders of magnitude higher than the Kerr nonlinearity of the fused silica ∼ 2 × 10 −22 m 2 /V 2 .In Fig. 3(a), where pump-frequency dependencies are shown, the billiard super-resonance described above is very clearly visible.This is not a unique property of metallic nanostructures.Giant nonlinearities in semiconductor nanostructures due to billiard resonances were recently predicted in [48].
Efficient frequency difference generation.As an interesting application we consider the process of optical rectification and difference frequency generation, governed by three nonlinearities χ (3) (δ; ω 0 , ω 0 , −2ω 0 + δ), χ (3) (δ; ω 0 +δ, ω 0 , −2ω 0 ) and χ (3) (δ; ω 0 , ω 0 +δ, −2ω 0 ), with a two-color optical pump at around ω 0 and 2ω 0 and signal δ ω 0 in THz and MIR.We solve the propagation equations, assuming slowly varying envelope approximation and taking into account dispersion relations, but neglecting nonlinear effects for the pump waves because of very very small propagation distance (see Supplementary for details).Both χ eff and the linear susceptibility χ (1) eff are calculated from given linear and nonlinear properties of the nanoparticles (χ NP ) and host (χ h ) using the effective medium approach [49] (see Supplementary).By calculation of the linear properties the full linear susceptibility containing both confinement-based and Drude parts are included.As a host material, we take fused silica which possesses strong losses in the range between 30 and 40 THz (see Fig. 4), but otherwise is quite transparent [50].We consider the filling factor of f = 0.01 and neglect the nonlinearity of the host.Resulting effective linear quantities are shown in Fig. 4 and demonstrate the usual PR resonance at around 2.4 eV with the width of around 30 THz.The shortest pulses still supported by this resonance are around 30 fs in duration.
Assuming an exemplary pulse durations of around 30 fs, we must consider two regions for the pump where conversion works significantly different.For the signal in the THz range (δ/2π ≤ 30 THz), the frequencies jω 0 and jω 0 + δ (j = 1, 2) are both located within the spectrum of the pump.In contrast, for the signal in MIR range δ/2π > 30 THz, the components ω 0 + δ, 2ω 0 + δ are not within the pump spectrum anymore.This leads to different treatment of these two frequency ranges for the selected pulse duration (see Supplementary).
The resulting field amplitude at 0th harmonic A 0 is given in Fig. 4(b) for different parameters and for the pump amplitudes A 1 = A 2 = 10 10 V/m.This pump for 30-fs pulses corresponds to a fluence around 0.3 J/cm 2 , which is yet below the damage threshold of gold (around 0.5 J/cm 2 [51]) and of fused silica (around 1 J/cm 2 [52]).One can see that in THz range the signal amplitude reaches 5×10 8 V/m corresponding to efficiency of around 10 −3 .In MIR range, the amplitude can exceed 10 9 V/m, delivering efficiencies above the percent level.Moreover, the maximal efficiency is achieved at 100 nm propagation distance for THz signal and 1 µm for MIR signal.From Fig. 4(b) one can also see that the most efficient conversion is achieved for the pump frequency ω 0 centered at the PR (solid lines in Fig. 4).In this case, the coupling of the pump to the signal is most efficient.

Discussion and conclusions.
We showed that confinement-based energy levels in metallic nanostructures, representing an integrable (or close to integrable) quantum billiards, typically join together into a single "super-resonance", which position and width can be controlled by the geometry of the nanostructure.Whereas we focused here on (almost) integrable quantum metallic billiards, we anticipate richer resonance structure and control possibilities if truly chaotic billiards are considered.We analyzed the super-resonance, using the level statistics extended in comparison to typically used in quantum chaos theory.In the linear regime the ballistic super-resonances is "hidden" behind the much stronger Drude response, yet it manifests itself strongly in a giant nonlinearity.This nonlinearity can be in addition enhanced by interaction with plasmons and effectively used to down-convert light to THz and MIR range with high efficiency already after 100-nm distances, despite of huge linear and nonlinear losses.Our confinement-based framework might also be helpful in a deeper understanding of the recent experimental work on efficient THz generation in nanostructures [53,54] and paves a way to extend newly proposed electronic meta-devices [55] into the nonlinear regime.IB, AD and UM acknowledge support from the Deutsche Forschungsgemeinschaft (DFG) under Germany's Excellence Strategy within the Cluster of Excellence PhoenixD (Photonics, Optics, and Engineering -Innovation Across Disciplines) (EXC 2122, pro-jectID 390833453).AH acknowledges support from European Union project H2020-MSCA-RISE-2018-823897 "Atlantic".

SUPPLEMENTARY Simplified Hamiltonian and wavefunctions for the case of spherical particles
To approach the problem analytically, we consider a spherical metallic particle of the radius a (diameter d = 2a).Since we are interested in low frequencies, we neglect the interband transitions.We neglect the temperature effects assuming that N electrons in conduction band occupy all levels below Fermi energy E F (see Fig. 1 of the main article).The energy structure is calculated assuming one-electron approximation and corresponds to that of a single free electron in a infinite-strength spherical potential of the radius a.To take into account finiteness of the potential, the levels above the work function E W are disregarded.The validity of these approximations is justified below.
The corresponding single-particle eigenproblem can be then formulated as Ĥ0 ψ = Eψ, with H 0 = − 2 /2m e ∆ + V (r), V (r) = 0 for r ≤ a and V (r) = ∞ for r > a (m e is the electron mass).With the approximations above, we neglect various effects of electron-electron and electronion interaction such as interband transitions, electron heating, the change of the eigenstates due to the finite height of the potential, and other effects, which play only minor role at low frequencies and ultrashort sub-100-fs pulse durations.The validity of this approximation is discussed in Sec.below.As it will be shown there, our simple model is rather adequate for the parameters we consider, despite of its simplicity.The advantage of this approach is the possibility to determine the energy structure analytically.The corresponding eigenfunctions are combinations of spherical harmonics.The energies are defined as where n, l are quantum numbers, and α nl is the nth zero of the Bessel function of order l.In contrast to the Coulomb potential, there is no degeneracy in l.
The eigenfunctions of the problem described in the main article are: where n, m, l are quantum numbers, R nl (r) = √ 2j l (α nl r/a) √ a 3 j l+1 (α nl ) , Y m l are spherical functions (−l ≤ m ≤ l), j l is the spherical Bessel function of order l, and α nl is its lth zero.
The radial part of the matrix element is, for the allowed transitions l − l = ±1, and zero in other cases (here e is the electron charge).

Cylindrical geometry
In this section we determine the eigenvalues for the cylindrical geometry.We assume here that the main axis of the cylinder oriented along z-direction, and the light is assumed to be linearly polarized also in z-direction.Note the difference in the denotations with the part where propagation is considered: there, z-direction is the direction of the light propagation.
The eigenfunctions in (r, θ, z)-coordinates are [56]: where C nlm is the normalization factor, h is the height of the cylinder, αlm is lth zero of the Bessel function J m (x) of order m.These eigenfunctions are described by three integer quantum numbers: n describes the localization in z direction, m and l -along the orthogonal directions.
The energies of these eigenstates are given by the expression Because the light is assumed to be linearly polarized in z-direction, only z-components of the dipole moments play a role, and they can be calculated as: where we omitted for simplicity a pre-factor which comes from the normalization of the wavefunctions.

Derivation of Eq. 1 of the main article
As we see in the main article (see also below), the main role in the super-resonance is played by the transitions close to the Fermi energy E F .For not very small nanostructures, this corresponds to relatively large l and n.For large n and l an analytical estimation is possible.Based on this, the energy difference between the allowed transitions ∆l = ±1 near the energy E is Near the Fermi energy E ≈ E F , substituting Eq. ( 3), we obtain Eq. 1 of the main article, with ω conf identical to ∆E.
We note that Eq. ( 10) is valid for n l.In contrast, for n ≈ 1 we have α nl ≈ (2n + l), that is, the positions of the resonances would be shifted in this latter case by the factor ∼ π/2 to the lower frequencies.The resonances shown in Fig. 2 of the main article as well as in Fig. 5 coincide well with Eq. ( 10), indicating, that large n are involved (see also Sec. ).
In addition to the derivation of the position of the super-resonance, we are able to approximately deduce its "natural" shape and width.For this, see Sec. below.

Linear and nonlinear susceptibilities of a single nanostructure
The corresponding expressions for χ (1) and χ (3) are obtained by the method of iterations [46]: The evolution of the density matrix ρ in the presence of damping can be, under suitable approximations, written as where Ĥ is the full Hamiltonian and ρ (eq) is the equilibrium value for ρ, Γ describes the decay.The Hamiltonian consists of the action of the potential well H 0 (see the main article) as well as the action of the field V = erE (in dipole approximation).In terms of the eigenfrequencies of Ĥ0 , ω mn = (E m − E n )/ , Eq. ( 12) can be rewritten in terms of perturbations as: ρmn = −iω mn ρ mn − i [V, ρ] − γ mn (ρ mn − ρ (eq) mn ), (13) where ρ mn = m| ρ |n , ρ (eq) mn = m ρ (eq) n , γ mn = m| Γ |n , |m , |n are eigenstates of Ĥ0 corresponding to eigenvalues E m , E n (note that here, in contrast to previous sections, we denote by n and m the "multiindices", fully describing the eigenstate).One can obtain ρ iteratively, in the form of ρ = ρ (0) + θρ (1) + θ 2 ρ (2) + . .., where θ is a formal small parameter, assuming thereby that V is a small perturbation "of order of θ".As an initial approximation we obtain and ρ (n) is related to ρ (n−1) as: The nonlinear polarization of n th order is defined via i , E i are the components of the vectors P (n) , E, respectively.P (n) is given in terms of ρ (n) as P (n) = −eN tr ρ (n) r , where N is concentration of the particles.The expression for χ (n)  is obtained by comparing the two expressions for P (n) above.For the first order susceptibility χ (1) ij (ω p ) we thus obtain: We note once more that a and n are "multiindices", that is, every of them describes a particular set of quantum numbers n, l, m fully characterizing the system; ω mn = E mn / , γ mn = δ mn γ (δ mn is a Kronecker symbol), γ = 1/T 2 ; T 2 and T 1 and are given in the main article.For the isotropic case we consider here, the indices i, j in Eq. ( 16) are disregarded.According to Eq. ( 14), ρ (0) ll describes the unperturbed populations (see more details below).The first term in Eq. ( 16) describes the Drude dispersion.It must be introduced into Eq.( 16) as an additional phenomenological term since its proper firstprinciple treatment is possible only if electron-phonon interactions [57] are taken into account, which is not the case in our approach.
For the third-order susceptibility we have: where P I denotes permutations of the frequencies ω p , ω q and ω r with the Cartesian indices h, i, k permuted simultaneously.
As it was mentioned in the main article, every term in the expressions for χ (1) and χ (3) can be seen as a sum over all transitions through the intermediate virtual states [46], with the initial and final state being the same.Since we do not consider effects of finite temperatures here, the initial populations are taken in the form: Moreover, in Eq. ( 18), due to the Pauli principle, we keep only the transitions over the intermediate virtual levels which are outside of the "Fermi sea", that is, with , where E 0 is the ground state.In Eq. ( 16), in contrast to Eq. ( 18), this pre-selection happens automatically.To take into account the finite depth of our potential, we also do not consider levels with E > E W +E F +E 0 , where E W is the work function.For this article, we have taken E W = 5.1 eV, E F = 5.53 eV.
For practical computations, in Eq. ( 16) and Eq. ( 18) we use the radial parts of the dipole moments given by Eq. ( 5), and average over angular parts [16].This is possible if assuming that only the transitions with l 1 are relevant, which is indeed the case even for smallest diameters we consider, as one can see in Fig. 1 of the main article.In this situation, averaging over the angular dependencies for −l ≤ m ≤ l gives [16] the constant factor A = 1/3 for the linear susceptibility χ (1)  and A = 2/15 for third order susceptibility χ (3) .Fur-thermore, when calculating χ (3) , we take into account a population-induced correction factor T 1 /T 2 (see [16]).As a result, the susceptibilities obtained in Eq. ( 16), Eq. ( 18) are corrected as: Linear susceptibility for different diameters Whereas in Fig. 1(a) of the main article only two particular examples of the linear susceptibility for spherical nanoparticles are shown, in Fig. 5(a,b) more examples are given, to illustrate further the dependence of the superresonance position on the particle size.

Influence of other mechanisms
Although the approaches presented in our article is rather universal and largely material-independent, in the main paper we, to be specific, considered the parameters of gold as our basic case since gold has the most practical importance.However, gold is quite a "complicated" metal in the sense that many other mechanisms contribute to nonlinearity.Besides, in all metals, in addition to quantum-confinement (billiard-part) the Drude part of χ (1) is contributing.In this section we will clarify in more detail the question how other mechanisms influence the overall nonlinear and linear response.Drude part of χ (1)   In the Figure 2 of the main article we show the linear response without the influence of the Drude part of χ (1)  (the first part in the left side of Eq. ( 16) above, also described by Eq. ( 17) above).The influence of the Drude part is much larger that the confinement-based part, especially at small frequencies in MIR and THz range.In particular, Fig. 5(c,d) shows the Drude part alone and together with the quantum confinement part of the linear susceptibility for few particular diameters.One can clearly see that the Drude part absolutely dominates in the linear response at low frequencies.However, this does not mean that that the confinement-based superresonance considered in this paper "disappears" as we take into account Drude.It can be still deconvoluted and separated from the Drude part [34].Besides, at it was shown in the main text, it is directly visible when considering the nonlinear optical properties such as Kerr effect.

Effect of interband transitions on the nonlinear response
In the main article, we considered rather simplified bandgap consisting only of one band.Whereas such approximation is very well suitable for some materials such as alkali metals, for other materials such as gold it could be claimed to be a rather bad approximation.However, at least in the particular case of gold, as soon as we consider low-frequency response, it can be shown that the nonlinearity due to interband transitions is much smaller than the one due to the billiard resonances.
In the case of gold [59], there is a strong resonance in the optical response around 2.4 eV, responsible to the transitions from the 5d valence band to the 6sp conduction band, as well as a number of less pronounced resonances in the range between around 2 to 10 eV.That is, the effects of the interband transitions might be pronounced for the photon energies above 1 eV.Even at the frequency resonant with the interband transition the impact of confinement-based resonances to the rectificationlike FWM processes we consider is one or two orders of magnitude larger that the effect of the interband transitions, as we will see in the next paragraphs.
Let us first consider the Kerr nonlinearity.Experimentally measured Kerr susceptibility for the photon energy around 2.3 eV, that is, close to resonance of the above mentioned transition is (for short, 100fs-scale pulses) χ (3) (ω; ω, ω, −ω) ∼ 10 −18 m 2 /V 2 (see for instance [11]; note that in many other references long, picosecond pulses are considered, demonstrating higher nonlinearity as discussed in the subsection below).These measurements, of course, include all effects simultaneously, in particular intraband transitions and confinement-induced effects.We see that our calcula-tions give the same order of magnitude of susceptibility for this frequency [see Fig. 3(a) of the main article] with only the confinement-based nonlinearity included.This means that the interband transitions do not dominate the intraband even at the interband resonance frequency; the impact of confinement-based resonances is of the same order of magnitude or higher.
On the other hand, as we decrease the frequency from the interband resonance towards THz range, the influence of the interband resonance quickly decreases whereas the influence of the confinement-based resonance increases [see Fig. 3(a) of the main article].Therefore, we come to the conclusion that in the THz range the Kerr nonlinearity χ (3) (ω; ω, ω, −ω) is indeed dominated by the confinement-based (intraband) transitions.
The same is also true for the FWM processes responsible for rectification-like effects considered in the main article: the influence of the interband transitions must be significantly smaller than of the intraband (confinementbased) ones.In order to make our estimation more quantitative at this point, we use the estimation technique described in [60].We start from the Kerr process χ (3) (ω; ω, ω, −ω) assuming it to be fully resonant to the interband transition ("worst-case scenario," where the interband action is maximal).We consider then the processes χ (3) (δ; ω + δ, ω, −2ω) and χ (3) (δ; ω, ω, −2ω + δ).Because of the missing resonant terms in Eq. ( 18) (which gives a factor ∼ ω 2 in comparison to the interband resonant case) and also because of smaller population of the virtual levels (factor ∼ T 2 2 ), the nonlinear susceptibility is reduced by a factor (ωT 2 ) 2 ∼ 10 2 comparing to the Kerr susceptibility at the interband resonance.Since, as it was established before, even at the interband resonance the confinement-based Kerr nonlinearity is at least of the same order of magnitude that the interband-induced Kerr nonlinearity, we therefore conclude that for the FWM processes the intraband (confinement-based) resonances are at least by the factor of 100 larger than the interband ones.Based on our calculations of the intraband nonlinearities [see Fig. 3(b) of the main article], we can estimate the interband effect to the nonlinear susceptibility for the considered FWM processes as 10 −20 -10 −21 m 2 /V 2 for ω at the interband resonance (and even lower away from that resonance).This is much less than the confinement-based impact as shown in Fig. 3 of the main article.
As an alternative and fully independent method to estimate the impact of the interband transitions we introduce a gap directly into our numerical model.That is, we modify our single-band structure as the following: where Λ is the lattice constant (for gold Λ ≈ 4 Å), E g ≈ 2.4 eV.This modification mimics the bandgap which opens near the edges of the Brillouin zone.The comparison of two calculations, that is, using the single band Eq. ( 2) model and two band model Eq. ( 21), are shown in Fig. 6 for different exemplary diameters.One can see that, whereas the interband transition does provide some limited modification to the confinementbased dynamics for very small nanoparticles d 3 nm, this influence quickly decreases and becomes negligible for larger diameters.
Thermal, hot-electron and related effects Taking into account temperature introduces several effects, which were neglected in the main article.First of all, it leads to an additional hot-electron contribution into nonlinearly [11], which can overcome, by several orders of magnitude, the nonlinearities considered in this article up to now.The hot electron mechanism involves laserinduced intraband excitation in the conduction band, followed by the energy dissipation of the excited electrons.This process leads to a modification of Fermi-Dirac distribution which depends on the frequency and intensity of the pump, leading thereby to frequency-dependent nonlinearity [61][62][63].The key role in the quick thermalization is played by the electron-electron and electron-phonon interactions.
This nonlinearity has relatively slow, sub-picosecondscale, turn-on time [64], and therefore its influence quickly decreases with decreasing of the pulse duration [11].For the pulses considered here (10-30 femtoseconds) this nonlinearity plays a negligible role.Indeed, the experimentally measured Kerr nonlinearity for gold nanospheres [11,65] for the pulses of 100 fs duration corresponds, by the order of magnitude (∼ 10 −18 m 2 /V 2 , see also discussion in the previous subsection) to our calculations at around the same frequency (cf.Fig. 3(a) of the main article).Since in our calculations we do not take into account the hot electron nonlinearities, we come to the conclusion that for short pulses such nonlinearities are pretty much negligible in comparison to the confinement-based resonances, or at least do not play a dominant role.This is even more true if we consider lower frequencies towards THz range, since the confinement-based nonlinearity has a resonance at low frequencies, whereas the thermal nonlinearity is not expected to demonstrate a resonant behaviour.
A part of the above-described thermalization process is the electron-electron interaction.Fast thermalization, indeed, is the primarily consequence of the electronelectron interactions [63,66].Electron-electron interactions are trackable in the linear properties of the nanostructure (see for instance [61]), however, the corresponding modification is rather minor and even this small modification starts to be visible at the time scale of few tens of femtoseconds.Another thermal effect is the overall non-rectangular shape of the electron distribution near the Fermi zone edge as soon as the temperature is nonzero (in our calculations in the main article we assumed zero temperature).The effect of nonzero temperature in the vicinity of Fermi-level is shown in Fig. 7 for an exemplary diameters d = 2.75 nm and temperature T = 300 K.It is obtained by modifying ρ (0) ll from Eq. ( 19) to the Fermi-Dirac distribution for the finite temperature.One can see that this modifies only slightly the linear response.The nonlinear response is also modified quite moderately.

General definitions
In this section we describe different variants of the level statistics, extending the discussion related to Fig. 2 of the main article.To understand the mechanisms, governing the formation of one single super-resonance it is very constructive to consider the level statistics, varying the selection rules included into that statistics.Commonly in the quantum billiard and quantum chaos theory [36] one considers the so called neighboring level statistics.That is, we consider the difference between the neighboring levels δω i = ω i+1 − ω i , where the eigenfrequencies ω i are obtained by ordering of the eigenfrequencies in the increasing order, that is, we order them in such a way that ω j ≤ ω i for j < i.In the case of nanospheres the corresponding eigenvalues are ω nl = E nl , cf.Eq. ( 2), reaaranged accordingly.Note that i, j here are single indices (and not multiindices as in Sec. ).Then, the probability P (ω)dω that δω i is located in the range between ω and ω + dω is calculated.The easiest way to visualize such statistics is to use the histogram technique.For nanospheres this traditional neighboring level statistics is presented in Fig. 2a of the main article (yellow bars).
The neighboring level statistics do allow to determine the universality classes of different billiards.It, however, does not take into account the properties of eigenfunctions and of the transition rules, so it seems to be of little use for the optical properties.In order to improve usability for optics, we extend the statistics to take into account all transitions, not only neighboring, and impose additional selection rules.That is, we consider the statistics of energy differences ω ij = ω i − ω j for i > j (assuming the ordering of ω j as discussed above).In addition, when constructing the probability density P (ω) of ω ij being in the interval [ω, ω + dω], we take into account the "strength" of the transition by weighting P (ω ij ) with a weight where µ ij is the dipole momentum of the corresponding transition i → j, and ∆ ij is defined as that is, takes into account the Fermi-Dirac distribution Eq. ( 19) (we call it below "the Fermi sea condition").The resulting statistics is shown in Fig. 2 (red and blue bars) and repeated for convenience as the inset in Fig. 8(a) -for the larger frequency range.We note that the traditional statistics discussed in the previous paragraph is a partial case of this more general approach; namely, we obtain the traditional neighboring level statistics assuming w ij = δ i,j=i+1 .

Natural shape and width of the super-resonance for nanospheres
To find out the "natural width" of the super-resonance for the case of nanostructures, we use a more more precise version of Eq. ( 10): which includes corrections s nl to the values of the roots.All transitions contributing to the super-resonance are characterized by the same value of 2n + l before and after the transition, that is, 2n + l = 2n + l + 1.However, values of n and l can be different.Therefore the energy difference between the eigenfunctions |ψ nl and |ψ n l will be modified as follows: For the typical range of n and l actual for our particular situation, we approximate s nl as being distributed in the range [-0.25,0].For the probability distribution of the difference s nl − s n l , given by P (s n l − s nl ) ∼ P (s)P (s + s n l − s nl )ds, we obtain a triangular sym-metric shape with the maximum at zero and full width of 1/2, that is we have a constant FWHM of 1/4 of the distribution of the difference s nl − s n l .Using Eq. ( 11) we obtain that this translates to the width ≈ ω conf /4 in frequency, where ω conf is the position of the super-resonance as defined in Eq. 1 of the main article.If we consider the dipole-momentum-weighted statistics as described below, this simplified conclusion will be modified by the fact that the dipole momentum depends on the energy level difference.Such dependence will result in sharper lower-frequency shoulder of the statistics and smoother higher-frequency shoulder.Both such shape and the width of the peak in the statistics are in a surprisingly good agreement with the findings shown in Fig. 2 of the main article (for the case of large T 2 ).Of course, reducing the T 2 to room-temperature values additionally increases the width of the super-resonance beyond the natural width of ≈ ω conf /4.

Influence of different mechanisms on the statistics
Varying the selection rules incorporated into the weighting Eq. ( 22), we can study how these rules influence the statistics.Different variants of statistics are shown in Fig. 8.In particular, in Fig. 8(a) we consider P (ω ij ) weighted by Eq. ( 22) with ∆ ij = 1, that is not taking into account the Fermi sea condition -i.e.all transition, not only the transitions from the below to above the Fermi sea level, are allowed.This situation describes dielectrics rather than metals, and should be contrasted to the "metallic case", that is, the situation where also the "Fermi sea condition" is satisfied [inset to Fig. 8(a) as well as Fig. 2a of the main article].One can see that in the former case the resonance is much more broad and, in addition, noticeably shifted to the lower frequencies.
On the other hand, if we remove all restrictions at all, that is, consider w ij = 1 for all i, j, we will see very broad distribution over the scale of many eV Fig. 8(b).The statistics changes not too significantly if we take into account only allowed transitions (but not yet distinguishing between the strengths of the transitions, that is, assuming w ji = 0 if µ ij = 0 and w ij = 1 otherwise, and also not taking into account the Fermi sea condition).Yet, in this case the peak, corresponding to the low-frequency super-resonance (at around 0.5 eV) does already appear [see Fig. 8(c)].This peak becomes even sharper if we in addition take into account the Fermi sea condition as shown in Fig. 8(d).But, in addition to this low-frequency peak, in Fig. 8(d) we see also many other peaks at higher frequency.These many peaks disappear as we take into account the "strengths" of the transitions (µ 2 -weighting as given by Eq. ( 22)), see inset to Fig. 8(a) and Fig. 2a of the main article.
Therefore, we conclude that there is, in general, a lot of possible confinement-based transitions at low and high frequencies.The key role in the formation of a single super-resonance is played by the dipole moments, that is, by the symmetry and composition of the wavefunctions, and only to a lesser extent by the Fermi sea condition.As we take them into account, we are left with a single bunch of closely spaced resonances, which, taking into account that every of these resonances are broadband, merge into a super-resonance.Whereas the consideration here was focused on the case of nanospheres, we remark that, most probably, this particular situation is quite geometry-independent, at least for regular billiards.Indeed, in this case we expect that distantly spaced eigenfunctions have very different number of oscillations in every spatial direction, making µ ij small.This, however, is not necessarily the case for the irregular, chaotic billiards which have also rather chaotic eigenfunctions.Finally, as Fig. 2(b) of the main article shows, we can tune the positions of the resonances in such a way that the resulting super-resonance is broadened.

Effective properties of nanocomposite
The effective linear properties of the nanocomposite for arbitrary frequency ω are calculated using the effective medium approach [49] using the linear properties of the host ε h (ω) (in our case SiO 2 ), and nanostructures ε NP (ω) [which is given by ε NP (ω) = 1 + χ (1) , where χ (1)  is calculated according to Eq. ( 16)], as: where f is the filling factor, x(ω) is defined as We note that at Mie resonance |x| is especially large.Finally, the effective nonlinear susceptibility χ eff for every process is calculated for given linear and nonlinear properties of the nanostructures and host as follows: NP x(ω 0 )x ( ω 1 ) 2 x(ω 2 ), (28) where χ NP is the nonlinear susceptibility given by Eq. 5 or Eq. 6 of the main article, with χ (3) in those expressions calculated using Eq. ( 18).

Propagation equations
Assuming slowly varying envelope approximation and neglecting nonlinear effects for the pump waves, the governing equations are: where c is the speed of light in vacuum, ε 0 is vacuum permittivity, A i , α i k i , i = 0, 1, 2 are correspondingly the (a) by yellow bars and coincides well with the Poissonian distribution (black dashed line).

FIG. 2 .
FIG. 2. Billiard resonances and weighted level statistics for a golden sphere of few exemplary diameters d (a) and a cylinder (b) of the diameter d and height h.The blue and red curves in (a) and (b) show Im (χ(1) ) as the function of frequency for different sizes (see legend).Bars of the same colors show the corresponding weighted level statistics P (ω) normalized to some P0 (P0 is not the same for different curves).Green bars show the statistics for all allowed transitions attached, for visibility, to the upper x-axis (that is, technically, 1 − P (ω)/P0 is shown; see also larger frequency scales in Supplementary).Black lines show Im (χ(1) ) for T2 increased 100 times.Yellow bars show NLS and dashed black line indicates the Poissonian statistics.

FIG. 3 .
FIG. 3. Nonlinear susceptibility for the Kerr nonlinearity (a) and for different FWM processes leading to the optical rectification (b,c) in dependence on the signal δ (c) and pump ω (b) frequency, shown for the nanostructures of different diameters.

FIG. 4 .
FIG. 4. Efficient generation of THz and MIR light in a composite of gold spheres with 2-color pump.(a) effective refractive index n eff and effective losses α eff in dependence on frequency ω for a gold nanostructures with r = 10 nm immersed into fused silica (f = 0.01).Vertical lines show two variants of the 2-color pump at ω0 = 1.55 eV (λ0 = 530 nm, dashed line) and ω0 = 2.4 eV (λ0 = 530 nm, solid line).The horizontal lines connect the spectral components of two-color pump.(b) The generated field amplitude for different propagation distances L and nanostructure diameter d (see legend) and the pump as described in text.Solid vertical line in (b) separates THz from MIR band.

FIG. 5 .
FIG. 5. (a,b) Linear susceptibility due to confinement-based resonances only, in dependence on frequency for spherical nanoparticle of different diameters.(c,d) The Drude contribution as well as the full linear susceptibility (Drude + confinement) for few selected diameters.

FIG. 7 .
FIG. 7. Real (a,c) and imaginary (b,d) part of the linear (a,b) and nonlinear (Kerr) (b,c) susceptibility for the structure d = 2.75 assuming sharp transition at E = EF (blue lines) and taking into account Fermi-Dirac distribution for T = 300 K (orange lines).