THz Generation via Optical Rectification in Nanomaterials: Universal Modeling Approach and Effective χ¯¯(2)$\bar{\bar{\chi }}^{(2)}$ Description

Optical rectification (OR) at the nanoscale has attracted an increasing interest in the prospect of providing efficient ultracompact terahertz (THz) sources. Here, a universal modeling approach capable of addressing both isotropic and anisotropic all‐dielectric nonlinear nanomaterials on an ultra‐broad spectral range, covering the highly dispersive phonon‐polariton window, and different orientations of the crystallographic axes with respect to the geometry of the structure is reported. This analysis is exemplified by considering two study cases, that is, nanopillars of AlGaAs and of LiNbO3. A close comparison between the two cases is established in terms of THz generation efficiency from 4 to 14 THz. Phonon‐polariton contributions to the OR process are disentangled from the electronic one, and a model order reduction based on the reciprocity theorem is applied and validated on both the considered configurations. These results, combined with the inspection of the THz near‐field features, pave the way to the design and optimization of nonlinear metasurfaces for THz generation and detection at the nanoscale.


Introduction
The invention of the laser in 1960 led to a revolution in the development of optical systems that were able to overcome the limitations of high-frequency electronic devices (e.g., in terms of power spectral density or data transmission rate) and that opened a new era in many fundamental and applied research fields such as atomic physics or nonlinear optics. [1]ince then, different spectral domains within the electromagnetic spectrum have been successfully addressed, but the terahertz (THz) region lying between the microwaves and the mid-infrared has undoubtedly been the least explored one, owing to important limitations imposed by material properties or device dimensions at such frequencies (roughly between 0.1 and 30 THz).In the last two decades, though, significant progress has been experienced by the THz region due to the increasing interest in investigating a vast variety of physical phenomena that overlap with such a spectral range, [2] including: i) collective lattice vibrations (phonons) in crystals; ii) molecular fingerprints of organic materials; or iii) charge carrier dynamics in semiconductors.[5][6][7][8] For instance, the large contrast of THz absorptance of different materials enables the use of such radiation in imaging systems at airport security checks or for medical purposes (non-invasive monitoring of damaged human tissues or characterization of pharmaceutical tablets thanks to the low THz photon energy in the order of meV).Another exciting example involves nonlinear THz spectroscopy, [9] where conventional time-domain measurements are performed as a function of the THz field strength or where THz pumping schemes are employed not only for exciting infrared-active phonons or molecular motions in polar media but also to analyze the coupling of such driven modes with electronic degrees of freedom by means of an optical probe. [10]The aforementioned usages of far-infrared (FIR) radiation demand the availability of high THz pulse energies [9,10] (in the order of millijoules), electric field peaks of tenths of MV cm −1 or output powers higher than 10 mW in continuous-wave (CW) operation; [6] all this, if possible, not being at the expense of reduced bandwidth, increased size, and non-effective costs.As a result, the THz community has witnessed a very rapid increase, particularly in the research of THz sources, which can be classified according to different criteria like the resonant/non-resonant character of the excitation [11] the pulsed/CW operation regime [3,4] or in terms of the underlying physical mechanism. [10]Following this last approach, THz emitters based on photoinduced transient currents, [12] spin-charge transfer, [10] and ionized airplasma [13] are extensively explored in the literature.In addition, a relevant framework for THz generation is the one that deals with the optical rectification (OR) of a femtosecond near-infrared (NIR) pulse; OR is a second-order nonlinear process that consists in the difference-frequency generation (DFG) between the intrapulse spectral components separated by a THz frequency.The first OR measurement in an electro-optical crystal (namely, potassium dihydrogen phosphate or KDP) was reported by M. Bass and colleagues as early as 1962, [14] succeeded by more sophisticated schemes based on Cherenkov-type geometry, tilted pulse front and quasi-phase matching in bulk and waveguiding structures. [11,15,16]Recently, nonlinear optics has benefited from the advance in fabrication techniques [17] to manufacture low-dimensional artificial materials (metamaterials) that allow to circumvent the efficiency-limiting constraint of optical-THz momentum matching in larger microstructures.In this spirit, THz generation by OR in plasmonic (i.e., metallic) metasurfaces has been proven to provide enhanced performance [18] due to the resonances held in their constituents or meta-atoms, though such nanostructures still suffer from small interaction volumes, low damage thresholds, and intrinsic ohmic losses.These challenges can be tackled by means of semiconductor materials, [19] and in particular with those possessing a strong  (2) nonlinearity; among these, aluminum gallium arsenide (AlGaAs) is a promising candidate for second-or third-harmonic generation [20] as well as for FIR emission. [21,22]Actually, the transverse-optical (TO) phonon resonances of AlGaAs fall within the THz band, giving rise to a resonant linear and nonlinear response that boosts the conversion efficiency up to seven orders of magnitude [22] with respect to their metallic counterparts.Moreover, ferroelectrics [23] such as lithium niobate [24][25][26] (LiNbO 3 or LN) are also attractive in view of their large nonlinearity, broad transparency region (from 400 nm to 5 μm) and extensive commercial availability.
Such an evident breakthrough established by nonlinear metaoptics for THz generation at the nanoscale has stimulated the realization of a large number of OR studies, comprising mainly experiments [27][28][29] and some fewer theoretical descriptions. [30,31]he main goal of this work is to provide a consistent and universal model of THz radiation from all-dielectric resonant nanoscatterers that could be used to design and optimize the functionalities of THz metasurfaces.The structure of the paper is as follows: In Section 2, the recipe for evaluating the pump-induced nonlinear THz currents is reported, based on a generalized Miller's rule for the  (2) susceptibility tensor of non-centrosymmetric media.This recipe is followed in the subsequent Section 3 using finite element method (FEM) simulations, accompanied by an analysis of: i) THz emission tunability with crystal orientation and ii) THz near-field response.Section 4 aims at reproducing the results of previous sections from an effective model of the nano-objects, exploiting the Lorentz reciprocity theorem.We conclude with a brief discussion and some general remarks in Section 5.

Universal Modeling of Optical Rectification in χ( 2) Nanomaterials
A key step toward the design of THz sources based on OR in nanomaterials is to model the interaction of an isolated nanostructure with an intense optical pump beam.The main requirements for an efficient interaction can be identified by considering the wave-equation for the electric field including nonlinear terms up to the second order, which in the frequency domain reads as follows ( The wave-equation above describes, in the electric-dipole approximation and for non-magnetic media, the spatial variation of the THz-field E(r, Ω) generated by the OR of a NIR pulse.The induced nonlinear polarization density P NL acts as the local source of E(r, Ω), which depends on the tensor product (:) between the χ(2) nonlinear susceptibility and the 3D fields in the linear (, Ω − ) regime; note that, as |Ω| << ||, the frequencies of the two pump components interacting via DFG are opposite in sign and so P NL follows the intensity profile of the NIR pulse.With this, on the one hand, a resonant response to the pump beam is desirable, and on the other, a strong χ(2) will further enhance the nonlinear source within the interaction volume.Interestingly, the material response at THz frequencies is governed by both electrons and ions in the crystal lattice.Pursuing a classical approach introduced with Garret's, [32] Lax and Nelson's, [33] and Boyd and Pollack's [34] pioneering works, we phenomenologically describe the  (2)    tensor-elements of non-centrosymmetric media as a function of the linear ( χ(1) ) electronic and ionic susceptibilities evaluated at the frequencies of the three interacting fields.The idea is to model the medium as a collection of nonlinearly coupled harmonic oscillators, giving rise to an anharmonic contribution to the crystal potential energy that yields the following second-order response where (1)a  ( 3 ) is the partial (a = e electronic; a = i ionic) linear susceptibility at  3 along the principal axis  and  abc  is a modelparameter also known as the generalized Miller's delta, related to the expansion coefficients of the anharmonic potential (and so material dependent).The above expression is valid for any kind of second-order interaction; particularly, for the DFG between two spectral components belonging to the transparency region of the medium in question, we have | (1)e  ( 3 ) (1)e  ( 1 ) (1)e  ( 2 ) (1)i n  ( 3 ) (1)e  ( 1 ) (1)   The first term in Equation ( 3) corresponds to the purely electronic contribution to the DFG nonlinearity (thus real and dispersionless), while the second term accounts for the contribution of the N IR-active phonon modes of the crystal at THz frequencies ( 3 ).In the following, we apply this so-called generalized Miller's rule to the cases of AlGaAs ( 43m point group) and LN (3m point group), which are of great interest for electrooptics and harmonic generation.In order to find the values of the  abc  parameters (all lattice contributions are modeled with the same  i n ee  ≡  iee  ), we have fitted the THz dispersion curve of (2) to the experimentally accessible low-frequency (microwave) and high-frequency (NIR) limit values reported in literature. [34,35]As shown in Figure 1, the ionic contribution is responsible for enhancing χ(2) by one order of magnitude with respect to the background electronic nonlinearity at the zone-center TO-phonon resonance frequencies.As a matter of fact, note that for AlGaAs  iee   eee  < 0, while for LN (except for the weakest  (2)   yyy ) the electronic and ionic contributions are in phase; this actually explains the reversed dispersion of the χ(2) of AlGaAs in the THz with respect to that of a classical Lorentzian χ(1) as well as the fact that its optical-mixing nonlinearity is stronger than the electro-optic one (as opposed to LN).
Knowing the spatial distribution of the electric field at the pump frequencies  1 =  and  2 =  − Ω, together with the corresponding (2) , allows to solve the waveequation (Equation ( 1)) for the generated THz field.Notably, the focus of the present work will be placed on the far-field (r = r ff >>  3 ) characteristics of the scattering media, which, as detailed in the subsequent Section 3, will consist of nanopillars with dimensions from ≈  3 ∕200 to ≈  3 ∕40 in the 4 − 14 THz frequency range.

Optical-to-THz Conversion Process
Let us consider the particular case where the incident electric field consists of a superposition of two monochromatic plane waves with NIR angular frequencies  1 and  2 ; in such a scenario, the nonlinear polarization density in Equation ( 1) oscillating at the difference-frequency  3 =  1 −  2 takes the form of where u  is the unit vector along the crystallographic direction  and E , (r,  1,2 ) are the amplitudes of the ] time-harmonic field components of the pump.The spatial dependence in Equation ( 4) arises from the fact that the source P NL (r,  3 ) is induced inside the nonlinear medium, which, in our case, involves a nanostructure that inhomogeneously distributes the fields at  1,2 within its volume.
In particular, we will concentrate on the analysis of the optical response of nanocylinders having an aspect ratio of H∕R = 2 and made of i) AlGaAs and ii) LiNbO 3 .
The zincblende AlGaAs is characterized by a cubic conventional unit cell where III-and V-group atoms are located at two FCC lattices displaced by 1/4 of the diagonal (see Figure 2a).It is thus optically isotropic, and its direct band-gap wavelength is  Γ ≃ 0.75 m (for 18% Al content [35] ).Conversely, LN (below its Curie temperature T c = 1483 K) belongs to the trigonal crystal system, and its hexagonal primitive unit cell is shown in Figure 2c.Along the [0001] or c-axis, Nb and Li cations sit inside oxygen octahedra but are not placed neither in the center nor at an oxygen plane; [36,37] this displacement between oppositely charged sublattices gives rise to a spontaneous polarization along the caxis and explains the ferroelectric nature of LN.It is important to stress that the (uniaxial) anisotropic properties of LN (e.g., the χ(2) tensor elements of Section 2) are commonly expressed in an orthogonal XYZ basis that is defined as follows: X is chosen along one of the three equivalent directions of the hexagonal unit cell base ([2 11 0] in Figure 2c), Z is chosen along [0001], and Y is then fixed to form a right-handed frame.Hence, with this convention, Z is the extraordinary (e) axis and XY is the ordinary (o) plane.
In order to construct the nonlinear source of Equation ( 4), we first study the linear response of our nanocylinders to the plane  [38] software.The linear scattering efficiency in the NIR region (black curve) is shown in (b) for AlGaAs and in (d) for LiNbO 3 , together with the corresponding multipolar decomposition. [39]The insets in (b) and (d) show, in the xz-plane, the distribution of the near-field enhancement due to the magnetic dipolar mode excited at  1 = 1.53 m.
wave excitation at frequencies  1,2 by means of finite element method (FEM) simulations performed in COMSOL Multiphysics 6.1.The air-surrounded nanopillars are designed in such a way that the z (lab) axis is orthogonal to the crystal cut planes ((001) for AlGaAs and XY for LN) and that the x (lab) axis coincides with the [110] and X crystal directions, respectively (see Figure 2a,c).The NIR dispersion of the scattering efficiency (i.e., scattering cross-section normalized to the illuminated geometrical area R 2 ) due to a normally incident x-polarized plane wave is shown in Figure 2b,d for AlGaAs and LN configurations, respectively.Both dielectric nanopillars exhibit multiresonant behavior, resulting in strong field localization and enhancement in the bulk and near the surface.Decomposing the induced currents on a multipolar basis enables us to link such resonant scattering to the excitation of the so-called Mie modes [17] (large displacement currents appearing from phase-retardation) within the nanostructures.It can be observed that while AlGaAs yields a stronger field concentration than LN owing to a higher NIR refractive index (n AlGaAs ≃ 3.3 > (n LN ) o,e ≃ 2.2), the latter might provide an advantageous performance for spectrally broader pump signals.
We now aim at investigating the capability of such dielectric nanoantennas to operate as potential THz sources.To do so, we set the wavelength of one of the incident wave's spectral component  1 = 2c∕ 1 at the strongest magnetic dipolar (MD) resonance around 1.55 m, which has already proven to be efficient for enhancing nonlinear phenomena in similar structures.Then, we accordingly shift from the MD peak the angular frequency  2 <  1 of a second pump spectral component in such a way that, via DFG,  1 −  2 =  3 results in a THz frequency; by sweeping the value of this shift, we are able to address light generation in a finite THz bandwidth.In this way, combining in Equation ( 4) the linear field solutions of Figure 2 with the nonlinear susceptibility models of Figure 1, we perform full-wave FEM simulations following the method described in ref. [22].The results are reported in Figure 3, in terms of the conversion efficiency  DFG , defined as the radiated THz power normalized to the square of the input NIR power, to be used as a figure-of-merit (FOM) of the DFG process.
A highly dispersive THz generation is observed both for the Al-GaAs and LN nanoantennas.Two local maxima of the conversion efficiency ( DFG ≃ 10 −8 W −1 ) are obtained at around f 3 =  3 ∕2 = 8.4 THz and 11.1 THz for the case of AlGaAs, while the use of LN yields an efficiency peak of  DFG ≃ 10 −9 W −1 at f 3 = 11.64THz.Even if, at first sight, one could attempt to ascribe such a THz generation increase to the TO-phonon-enhanced χ(2) susceptibility, the aforementioned THz frequencies are shifted from those of the bulk media's f TO =  TO ∕2 (recall Figure 1).Actually, as already demonstrated in ref. [22] for AlGaAs, the latticeinduced negative real permittivity bands (Reststrahlen bands) at the THz region allow the excitation of localized surface phononpolaritons [40] (SPhPs) in the nanoscatterers, causing a resonant emission of THz fields at the SPhP-Fröhlich frequencies (see Figure 3c).This view is reinforced if one ignores the phononic contribution to the linear THz permittivity of AlGaAs and LN by setting it equal to its real high-frequency value  ∞ (Rayleigh THz scattering): the yellowish solid lines in Figure 3a,b    DFG would be largest at the TO-phonon frequencies, while the SPhP-related peaks would disappear.The latter is true for both dielectrics under study, although an interesting difference emerges from the stronger absorption lines of LN: at the lowest A 1 -phonon frequency around 7.5 THz  DFG is reduced by a factor of ≈ 600 between the yellow and blue curves in Figure 3b, due to the large Im[ e ] ≃ 190; [41] such a drop is of a factor of ≈ 30 in the AlGaAs nanopillar at its lowest 8 THz GaAs-like phonon mode, where Im[] ≃ 65. [42] The influence of TO-phonons in efficiency is further evidenced if we analyze the hypothetical  DFG ( 3 ) coming from a purely electronic χ(2) (red curves in Figure 3a,b): the ionic contribution to the nonlinearity accounted for in Equation (3) enhances the THz emission by up to 2-3 orders of magnitude at f TO , both for AlGaAs and LN nanopillars.
At this point, it is worth noting that there is a substantial change in the  DFG dispersion of Figure 3a with respect to that reported in ref. [22], on account of the different modeling approaches of the linear and nonlinear THz properties of AlGaAs.On the one hand, in the present work, the linear THz permittivity is described by an additive-Lorentzian model fitted to the empirical results reported in ref. [42].On the other hand, the χ( 2) is modeled as a function of the Lorentzian χ(1) via the general-ized Miller's rule, which invokes Faust-Henry coefficient values that differ from those employed in ref. [22] (see Sections S1 and S2, Supporting Information).In spite of the quantitative change due to a more penalizing scenario for the χ(2) of AlGaAs, we find that our current results emphasize the key role of SPhPs in the THz generation process, whose response is still able to peak  DFG .Last but not least, we have observed that the configuration explored for the LN nanopillar is fairly elucidated if, instead of all of the 11 non-vanishing nonlinear coefficients, only the ZZZ tensor component is taken into account in the simulations (blue dashed line in Figure 3b).A comparison with the full model calculations (blue solid line in Figure 3b) indicates that the  (2) ZZZ component is the dominant one.This will facilitate the interpretation of the far-field emission properties discussed in the following section.

Effects of Crystallographic Axes Rotation
The time-varying nonlinear polarization in Section 3.1 yields a bulk nonlinear polarization current equal to J NL (r,  3 ) = −i 3 P NL (r,  3 ), oscillating at THz frequencies.Ergo the modeled scatterers behave as sub-wavelength antennas radiating THz electromagnetic waves; in the far-field region, these are transverse fields with an associated emission pattern that depends on the (, ) direction.
The nonlinear currents induced by the optical pumps within the AlGaAs and LN nanopillars modeled in the previous section are shown on the left-hand side of Figure 4a,d.In both cases the x (and y) components of J NL (r,  3 ) in different regions of the pillars are the same in magnitude but (unlike the z component) oppositely directed, yielding a destructive interference of their far-field contributions.By analyzing the 3D emission diagrams calculated in COMSOL via the Stratton-Chu transformation (right panels in Figure 4a,d), it is confirmed that the farfield THz radiation is equivalent to that of a z-directed electric dipole current element.Such an in-plane emission could be useful for exciting THz surface waves in larger antenna arrays, [22] but it is also desirable to investigate the potential of the designed nanoantennas to scatter THz light toward other directions such as the forward-backward ones, which can be appealing for microscopy, sensing or telecommunications applications. [43,44]This issue can be tackled by growing the nanocylinders along different crystallographic directions of the nonlinear media.For Al-GaAs, we have studied the case where the <100>-direction is parallel to the nanopillar's axis (i.e., the so-called <100>-cut); for the zincblende structure, there are, however, two other possible non-equivalent crystal directions along which the nanoantennas can be fabricated: <110> and <111> (see Section S3, Supporting Information for more details).The FEM simulation results for the nonlinear THz emission of such configurations are shown in Figure 4b,c.The spatial distribution of the local nonlinear currents is now altered, and the THz far-field emission changes correspondingly: the THz radiation null is now tuned toward the x-axis (<110>-cut) and towards the (, ) = {(41 • , 90 • );(−139 • , −90 • )} directions (<111>-cut).For LN, in addition to the already studied Z-cut configuration, the X-cut, and Y-cut LN are also technologically relevant cleavage planes.The former is employed as a substrate in LN thin films [37] and it is often preferred in nonlinear metasurfaces [25,29] to exploit the strongest  (2)   ZZZ coefficient, while the latter is interesting due to its piezoelectricity and suitability for surface acoustic wave devices. [37,45]Interestingly, both the X-and Y-cuts also enable us to remove the THz emission null from the forward-backward directions, as demonstrated in Figure 4e,f.All of the above considerations correspond to the case where the optical pumps are polarized along the x-axis of the lab frame and are useful to reveal the possibility of engineering the nonlinear THz emission from nanoscale resonators.Their performance could in principle be further advanced by optimizing the polarization and/or angle of incidence of the pumps (apart from the geometrical parameters of the scatterers).For instance, the <111>-cut AlGaAs nanopillar not only allows to enhance the optical-to-THz conversion efficiency with respect to the <100>cut design (see Figure S2, Supporting Information) but also provides the possibility to differently orient the THz far-field pattern for different pump polarizations at constant efficiency, similarly to the case of SHG. [46]

Nonlinear Response in the Near-Field
So far we have discussed the far-field of the THz radiation generated by optical rectification in AlGaAs and LN nanoantennas.However, also the near-field THz response of χ(2) nanomaterials would be of interest in view of possible applications in the context of photonic chips, both for the generation and the detection of THz radiation.Integrated hybrid optical-terahertz photonics can pave the way to an unprecedented miniaturization of THz functionalities and to their integration in a platform being much more compatible with the fiber-optics technology, which is an undisputed driver of the photonic market (for an overview, see ref. [47], and references therein).Even though an exhaustive analysis of the nonlinear near-field response of χ(2) nanomaterials is beyond the scope of the present work, here we discuss basic features of the near-field formation for two different configurations of the nonlinearly interacting infrared and THz beams.
First of all, notice that despite the fact that the far-field radiation patterns of Figure 4 are almost identical to each other (apart from the rotation due to the different orientation of the crystallographic axes) and correspond to the radiation pattern of a point-like electric dipole, the spatial distribution of the local nonlinear currents J NL (r,  3 ) look differently.This is an indication that in the near-fields non-radiative contributions are present depending on the specific configuration of the nanopillar.
To highlight this issue, Figure 5a-d shows the 3D patterns of the near-field (colormap: electric field in normalized units; arrows: field orientation) nonlinearly generated at selected THz frequencies for (Figure 5a,b) AlGaAs and (Figure 5c,d) LN nanopillars with <100>-cut and Z-cut, respectively.Note that to avoid numerical artifacts in the near-field analysis all the edges have been rounded with a fillet radius of R∕4 (R being the nanopillar radius), which is also consistent with the typical fabrication deviations.For AlGaAs, we picked the frequency of the GaAslike TO-phonon at 8.02 THz (panel (a)) and that of the SPhP resonance at 11.1 THz (panel (b)).Note that, despite the different nature of the selected resonances (material resonance versus quasi-static nanostructure resonance), the near-field patterns retrieved at these THz frequencies look almost identical and the field is concentrated around the centers of the top and bottom bases of the pillar, with an orientation that is typical of an electric dipole aligned to the vertical z-axis.For LN, we considered the near-field patterns at two similar frequencies, with the lower sitting at the A 1 (TO)-phononic peak, that is, 7.43 THz (panel (c)), and the other peak of the SPhP resonance, that is, 11.64 THz (panel (d)).Note that the two patterns look now different from each other, which is due to the higher dispersion of the different  non-vanishing LN χ(2) tensor components (as detailed in Figure 1).Most interestingly, this contrast in the χ(2) structure of both materials is responsible for the dramatic difference that can be appreciated in the near-field patterns of the AlGaAs and of the LN pillar at the SPhP resonance (compare panels (b) and (d) in Figure 5): for the LN (AlGaAs) the internal field is high (low) at the center of the nanopillar and the field outside is concentrated on the rounded edges (on the top and bottom bases).
Near-fields in χ( 2) nanomaterials are of interest not only for the generation of THz radiation but also for its detection at the nanoscale.Our modeling approach is ready to explore this scenario.Actually, the optical rectification considered in previous sections is the conjugate process of the electro-optic (EO) effect, [48] which in turn is well known to be at the heart of the electro-optic sampling technique for THz detection (a simplified sketch of the latter is shown in Figure 5e).In a photonic integrated platform, this technique demands two key ingredients to be implemented: [47] i) an efficient electro-optic medium to be illuminated with an intense optical or infrared beam; ii) a resonant THz antenna capable of concentrating the THz field to be detected in the electro-optic medium, usually placed in the antenna gap.A χ(2) nanoantenna of the kind considered in the present study is capable of providing both such functionalities when operated at the SPhP resonance.To investigate this possibility we simulated the difference-frequency generation process, which, in a frequency domain description, will be one of the contributions to the electro-optic effect; [48] we thus considered an infrared beam at  1 =  MD as the optical probe and a THz wave at  3 ∕2 = 11.1 THz as the signal to be detected impinging on an AlGaAs nanopillar (both z-polarized and traveling along x), and we calculated the infrared light generated at the frequency The near-field distribution of the optical wave resulting from the electro-optic effect in the AlGaAs nanopillar (which in Figure 5 is named the EO signal to distinguish it from the incident optical probe) is shown in Figure 5f.As anticipated by the THz generation process, there is a field localization due to the magnetic dipole Mie resonance at  2 that enhances the near-field in the forward (x) direction, which could be beneficial for coupling the EO signal to the next stage of the detection system (e.g., optical waveguides).Moreover, Figure 5g shows that the nanoantenna receiver simultaneously acts as a THz antenna when the frequency of the incoming THz field matches the localized SPhP resonances of the material, further enhancing the EO signal at the forward point Q.These observations unveil the potential of the proposed all-dielectric nanopillars not only as efficient THz sources but also as ultracompact nanoscale THz detectors.

Model Order Reduction Based on Reciprocity Theorem
The nonlinear scattering problem of OR in a dielectric nanoantenna is sketched in the left panel of Figure 6a, and can be summarized as follows: thanks to the second-order nonlinear susceptibility χ(2) of the nanoantenna medium, the two incident optical fields (pumps) with angular frequencies  1 and  2 will give rise to a local current density J  3 a , oscillating at the difference angular frequency  3 =  1 −  2 (belonging to the THz spectral region) and thus emitting a DFG electric field E 𝜔 3 a at this frequency.Under suitable hypotheses, the modeling of this process can be reduced to a set of three linear scattering problems, one per each of the angular frequencies under consideration, and an effective nonlinear susceptibility tensor [49] χ(2) can be subsequently introduced, with benefits in terms of computational speed and physical insight.The first two steps, namely the calculations of the scattered fields at  1,2 , E(r ′ ,  1,2 ), upon illumination with incident fields E  1,2 inc , are the same as in the numerical studies reported in the previous sections.For the third step, we need to connect the actual nonlinearly generated local currents (J b , that can be determined from a linear scattering calculation, as sketched in the right panel of Figure 6a.This can be done by resorting to the Lorentz reciprocity theorem (LRT), stating that where the integration is performed on the whole 3D domain Ω where all the currents and fields are defined, and this domain is supposed to be filled with a reciprocal medium, having  ij =  ji (which, in general, can be inhomogeneous and anisotropic).For the current J (locally a plane wave).We then followed the same procedure detailed in ref. [50], with the main difference here being represented by the relationship between the scattered physical field E 𝜔 3 a and the local effective source.For the latter, in view of the radiation pattern retrieved by FEM calculations of Section 3.2 from our sub-wavelength scatterers, we assumed an effective point-like electric dipole with current dipole moment μ 3 (see Figure 6b), whose components are given by with j, m, and n the cartesian coordinates of the xyz lab frame, V c the nanoantenna volume and χ(2) the effective nonlinear susceptibility tensor, to be evaluated as The optical-to-THz conversion efficiency is finally determined as ), with the THz power emitted by such effective point-like source given by the Larmor's formula: [51] The results of these calculations are shown with cyancolored circular markers in Figure 3a,b for AlGaAs and LiNbO 3 nanocylinders, respectively.For the latter, we considered the  (2) ZZZ only component which turned out to be the dominant χ(2) tensor element for the DFG process under consideration, as discussed in previous Section 3.1.Note that the effective linear modeling retrieves a very good quantitative match with the full nonlinear model calculations (cf.cyan markers with blue curves in Figure 3a,b).The reduced model (valid for any kind of crystallographic orientation) is thus capable of an accurate prediction of the OR process by resorting to linear simulations that can be implemented in any electromagnetic computational tool.It is worth mentioning that, contrary to the case when the wavelength is comparable to the resonator's characteristic dimension, evaluating a single illumination condition at the linear THz step is enough to determine the dipolar emission in our case.The exploitation of such a model can thus significantly speed up the design of novel configurations of THz nanomaterials and boost developments in the field.For instance, the strong second-order response of AlGaAs/LiNbO 3 nanoantennas combined with sophisticated nanofabrication techniques bring all-dielectric metasurfaces to the spotlight of potential feeders for THz m-sized antennas; [25] the presented reduced model could be of great aid for numerical calculations of such extended platforms as it would allow replacing the nanoantennas by point dipole sources with current moment μ 3 in the computational domain.

Conclusion
We have theoretically investigated the generation of THz radiation by the difference-frequency generation between two spectral components of an NIR pump pulse in AlGaAs and LiNbO 3 Mie-resonant nanopillars.It is found that optical phonons play a substantial role in the frequency conversion process and the reason is twofold: on the one hand, the bulk nonlinear susceptibility χ(2) ( 3 ,  1 , − 2 ) becomes resonant at the TO-phonon frequencies and, on the other hand, the lattice-induced negative permittivity at THz frequencies gives rise to a resonant emission of THz light, due to the excitation of localized SPhPs.The dipolar THz emission from the nanopillars has allowed us to describe them as effective point-like THz sources placed at the center of the nanoantennas; we have demonstrated that such an effective picture, which solely demands linear simulations, is able to exactly reproduce the dispersion of the optical-to-THz conversion efficiency derived from nonlinear FEM calculations.The presented reduced model thereby not only provides a clear physical understanding of the phenomenon but also simplifies the numerical implementation of the down-conversion process in more complex mesoscopic THz photonic structures consisting of ensembles of such dipolar nanocylinders.Finally, in comparison to well-developed bulk schemes for OR in χ( 2) -crystals, the present approach evades the phase-matching requirements and limits the effect of phonon absorption at THz frequencies via nanostructuring, which in turn makes spectral and spatial mode overlaps decisive for efficient THz generation.These can be optimized by exploiting the degrees of freedom of the pump beam, as well as by controlling the crystallographic direction along which the nanoantennas are grown.Functionalities such as THz sensing or even detection could profit from the local field enhancement achieved due to the multiresonant character of the proposed nanoantennas.For instance, electro-optic crystals used in conventional THz set-ups (with dimensions in the order of 100 m) could be substituted by our much smaller meta-atoms, which can be the core of an optimized design with improved detection characteristics (bandwidth, sensitivity, dynamic range etc.).With this, we envision such resonators to be potential candidates for their integration on thin-film-based AlGaAs/LN photonic circuits [52,53] that have recently been proven to be essential for addressing society-driven applications.

Figure 2 .
Figure 2. Resonant linear optical response of all-dielectric nanopillars.(a) and (c) show the designed geometry and pumping configuration for the AlGaAs (R = 200 nm) and LiNbO 3 (R = 300 nm) nanopillars, respectively; the crystallographic unit cells are produced by using VESTA[38] software.The linear scattering efficiency in the NIR region (black curve) is shown in (b) for AlGaAs and in (d) for LiNbO 3 , together with the corresponding multipolar decomposition.[39]The insets in (b) and (d) show, in the xz-plane, the distribution of the near-field enhancement due to the magnetic dipolar mode excited at  1 = 1.53 m.

Figure 3 . 2 )
Figure 3. Generation efficiency comparison in the region between 4 and 14 THz.The simulated conversion efficiency curves for the a) AlGaAs and b) LN nanocylinders are derived for the following scenarios: full nonlinear problem (blue solid lines), NIR nonlinearity (red solid lines) and Rayleigh scattering at the THz (yellow solid lines) with ( 3 ) =  ∞ .An extra dispersion curve is included in panel (b) (blue dashed line), belonging to the case where only the  (2) ZZZ component is assumed to contribute to THz generation.The black dotted lines in (a) and (b) represent the real part of the dispersive ( 3 ) for both materials, while the cyan-colored filled circles belong to the  DFG values derived from the effective susceptibility method described in Section 4. Panel (c) shows the resonant THz extinction due to the excitation of SPhPs by an incident z-polarized THz plane wave.The inset within shows, for the AlGaAs nanopillar, the near-field enhancement of E z together with the direction of the scattered field (white arrows) in the xz-plane at around f 3 = 8.4 THz.The internal field is almost uniform, while the scattered one is localized close to the top and bottom surfaces (more pronounced on the cylinder edges).

Figure 4 .
Figure 4. THz nonlinear emission characteristics for different orientations of the crystal axes, with normally incident NIR pumps at  1,2 ≃  MD and polarized along x. a-c) Local nonlinear current density (normalized to its maximum value in each case) and corresponding 3D far-field radiation pattern of the AlGaAs nanopillar for the <100>-cut (a), <110>-cut (b), and <111>-cut (c), all at f 3 = 11.1 THz.d-f) Same for the LiNbO 3 nanopillar at f 3 = 11.64THz, in the case of Z-cut (d), X-cut (e), and Y-cut (f).

Figure 5 .
Figure 5. Nonlinear response in the near-field for both generation and detection schemes.a-d) Near-field pattern (|E| normalized to its maximum value) of the nonlinearly generated THz radiation at selected frequencies for <100>-cut AlGaAs (a,b) and Z-cut LiNbO 3 nanopillars (c,d).e) Sketch of the conventional electro-optic detection in a thick χ(2) crystal.f) Near-field pattern of the emitted infrared signal (EO signal) generated via DFG between the incident infrared beam and the THz radiation to be detected, when the latter is in resonance with the 11.1 THz SPhP of the <100>-cut AlGaAs nanopillar.g) Field spectrum (|E| normalized to the peak at 187.54 THz) of the DFG infrared near field evaluated in point Q of the panel (f).Arrows represent the orientation of the electric near fields outside the nanopillar.

Figure 6 .
Figure 6.Sketch of the reduced modeling approach.a) Conceptual representation of the two scenarios linked by the Lorentz reciprocity theorem.b) Definition of the effective problem, which replaces our nanoantenna by an effective point-like electric dipole at the THz frequency  3 , with current dipole moment μ 3 , located at the center of the nano-object (origin of the lab frame).

𝜔 3 a
) and scattered fields (E  3 a ) to a configuration of far-field currents, J  3 b , and subsequently scattered local near fields, E  3

𝜔 3 b 3 b = p 𝜔 3 b
, we took the one corresponding to a linearlypolarized uniform plane wave with electric field E  3 inc at the position of the nanoantenna; this can be done by selecting E  3 inc to be the far-field emitted along the direction of the radiation Laser Photonics Rev. 2024, 18, 2300669 maximum of an ideal point-like dipole source J  (r − r ff ), that is, E