Second harmonic generation from metallic arrays of rectangular holes

The generation process of second harmonic (SH) radiation from holes periodically arranged on a metal surface is investigated. Three main modulating factors affecting the optical response are identified: the near-field distribution at the wavelength of the fundamental harmonic, how SH light couples to the diffraction orders of the lattice, and its propagation properties inside the holes. It is shown that light generated at the second harmonic can excite electromagnetic modes otherwise inaccessible in the linear regime under normal incidence illumination. It is demonstrated that the emission of SH radiation is only allowed along off-normal paths precisely due to that symmetry. Two different regimes are studied in the context of extraordinary optical transmission, where enhanced linear transmission either occurs through localized electromagnetic modes or is aided by surface plasmon polaritons (SPPs). While localized resonances in metallic hole arrays have been previously investigated, the role played by SPPs in SH generation has not been addressed so far. In general, good agreement is found between our calculations (based on the finite difference time domain method) and the experimental results on localized resonances, even though no free fitting parameters were used in describing the materials. It is found that SH emission is strongly modulated by enhanced fields at the fundamental wavelength (either localized or surface plasmon modes) on the glass metal interface. This is so in the transmission side but also in reflection, where emission can only be explained by an efficient tunneling of SH photons through the holes from the output to the input side. Finally, the existence of a dark SPP at the fundamental field is identified through a noninvasive method for the first time, by analyzing the efficiency and far-field pattern distribution in transmission at the second harmonic.

Second harmonic generation is investigated in the context of extraordinary optical transmission. Two different regimes are studied in which enhanced linear transmission occurs either aided by surface plasmon polaritons or through localized electromagnetic modes. It is shown that the radiation pattern and the transmission efficiency of second harmonic light strongly depends on the electromagnetic modes excited at both the fundamental and second harmonic frequencies. Second harmonic generation is constrained by the spatial distribution of these fields and their ability to couple with radiation. In particular, the symmetry properties of SH fields are crucial in the near-to-far field coupling process. For instance, it is shown that light generated at SH can excite electromagnetic modes otherwise unaccessible in the linear regime under normal incidence illumination. By reciprocity, these modes can only couple to off-normal diffraction orders in the lattice.

I. INTRODUCTION
Second harmonic generation (SHG) is a non-linear process that creates a single photon at λ/2 through the interaction of two photons of wavelength λ. 1 Since its discovery, SHG has become a current toolkit in optics with potential applications in several areas: from research in biological imaging 2,3 to recent advances in quantum information through parametric down conversion [inverse second harmonic (SH) process]. 4 The SH fields originate from the bulk and the surface of many substances. In non-centrosymmetric materials, like quartz crystals, SHG mostly originates from bulk. In centrosymmetric materials, like metals, bulk generation is forbidden by symmetry, 5 but SH fields may generate at the surface where the inversion symmetry is broken. [6][7][8][9][10] Initially, most of the experiments investigated SHG from metal surfaces, but the advance in nanofabrication techniques and optical characterization at the nanoscale 11-16 has turned the attention to SH effects in metallic nanostructures. For instance, some theoretical predictions on SHG from spherical nano-particles, developed over more than two decades, [17][18][19][20][21][22] have been experimentally tested at single particle level only recently. 23 Nano optics know-how opens the possibility to future SHG-based applications for optical characterization. [24][25][26][27][28][29] The SHG process is very weak in a flat metal surface, so only high-intensity lasers provide enough output. Nanostructured metal films locally enhance the intensity of the incident field, which is useful to obtain SHG at less demanding laser powers. Taking advantage of the strong electromagnetic (EM) fields found in holey films, several attempts have been conducted to exploit their nonlinear response for harmonic generation. [30][31][32][33][34][35][36][37][38][39][40] Holey metal films have been widely investigated in the linear regime especially since the discovery of extraordinary optical transmission (EOT). 41 However, despite all the work done many questions are still unsolved about the nonlinear response of these systems. For instance, Nieuwstadt et al. 32 found SH enhancement due to localized modes occurring close to the cutoff wavelength of the fundamental harmonic (FH) field, λ c . That enhancement was explained in terms of slow EM modes localized in the holes at FH. But the conclusion has been challenged by the same authors, 42 and recent experiments have shown that the time delay at FH is similar for different aspect ratios. 40 In this work, we investigate SHG from two dimensional rectangular hole arrays (RHAs) in a square lattice carved on gold films. The source of FH fields is an external beam illuminating the system at normal incidence from the air side (see Fig. 1).
The paper is organized as follows: in Sec. II we describe the theoretical approach, which is based on the finite difference time domain (FDTD) method. 43 In Sec. III we discuss the main mechanisms governing SHG in RHAs and how these are dictated by symmetry considerations and the range of geometrical parameters. In Sec. IV we investigate the response at SH due to localized resonances at FH (we compare our numerical results with the experimental measurements in Ref. 32) and extend our analysis to surface plasmon polaritons (SPPs). We end up with the conclusions in Sec. V.

II. THEORETICAL APPROACH
We use the FDTD method for calculations. The linear response of gold is described by the Drude-Lorentz model presented in Ref. 44, which provides suitable functional fit to experimental data valid in the visible range. 45 To calculate the optical response at SH frequency we follow a perturbative approach, for which the fundamental field is not affected by the SH field (non-depletion approximation). This is an excellent approximation given that the radiated intensity at FH is roughly ten orders of magnitude larger than at SH. We solved simultaneously the system of coupled equations governing the propagation of the FH and SH fields. The equations for the FH field obviously coincide with the linear Maxwell equations when no SH field is present, while the SH field obeys the Maxwell equations with sources determined by the FH field through the second order polarization vector, P (SH) . The induced polarization at the SH frequency mainly originates from surface contributions in centrosymmetric materials, and it can be expressed as: 7,8,46,47 where n and t stand for normal and tangential to the surface respectively, and χ ijk are the non-vanishing components of the second-order susceptibility tensor. The contribution from bulk to SH is negligible compared with the surface one. [46][47][48] The FH electric field is taken at the metal surface, and from it P (SH) is calculated at the same location. 47 In FDTD, every spatial location in the system can be visualized as a cube with the electric field components pointing along the edges and the magnetic field components being normal to the faces. 43 The interface between two different materials is composed by adjacent faces, the electric field components lying on it. When a face rests on a metal surface the electric field is computed with the piece linear recursive convolution method 49 . Outside the metal, the electric field is updated as corresponds to a lossless dielectric. This method applies both to SH and FH fields, which avoids spurious asymmetries in the SH field. Every metal face belongs to a given FDTD cell and has associated a single value of P (SH) . The tangential electric field in Eq. 1 is obtained from the average of the electric fields located at the face edges. The normal component of the electric field in Eq. 1 is obtained from the average of the electric fields located at the edges perpendicular to that face, inside the metal.
We intentionally work at fixed wavelength and change the geometrical parameters defining the array, using realistic values similar to that of experiments. 32 In this way, we separate the spectral contributions to SHG arising from the geometry to those arising from the spectral dependence of the nonlinear polarizability (which is not well known).
For the sake of simplicity, we assume that the nonlinear response of gold is weakly dispersive around the frequency range investigated, and for definiteness we take χ In many circumstances, Eq. 1 can be simplified by neglecting the contributions to SH from χ (2) tnt and χ (2) ntt , provided SH radiation keeps unaltered under this approximation. This is the case here and up to 80% of the total radiated photons at SH originate from χ (2) nnn , as we will show later on. Under this approximation |E (SH) | ∝ χ (2) nnn , so our results can be extrapolated to other experimental measurements. The dominant tensor element preserves the full symmetry of the SH field and provides qualitative and quantitative results.
Finally, let us mention that we checked our code against analytical results obtained for flat metallic surfaces 50 obtaining accurate results for a numerical mesh size of 5 nm (not shown), which has been used in the following RHA calculations.

III. FAR-FIELD EMISSION AT SH: GENERAL PROPERTIES
Holey metal films display a complex optical linear response, which is mainly controlled by the geometrical parameters of the array, the optical properties of the metal and the refractive index of the surrounding medium. 51 Here, we describe the main differences between the linear and the nonlinear response, the most important: FH and SH fields have opposite parity symmetry. At the same wavelength, light inside the holes can couple to waveguide modes of different symmetry whether it originates from an external source (laser beam) or it is nonlinearly generated at SH. These EM modes are unaccessible in the linear regime at normal incidence, while they are allowed for oblique incidence. In the last case however, such a mode can not be isolated from other modes, which can have opposite symmetry. Therefore symmetry determines the ability of the local SH fields to couple with the propagating diffraction orders of the lattice. For instance, emission of SH light is only allowed for off-normal propagation at normal incidence. In addition, the radiation pattern and intensity at FH is essential to understand SHG, but also the propagation properties of the waveguide modes excited at SH inside the holes.
The intensity of SH radiation depends on both the material properties (linear dielectric constant and χ (2) ) and the geometry of the system in a complicated manner. In what follows, we describe three main factors that modulate the nonlinear response: i) coupling of SH light with the lattice diffraction orders (controlled by the symmetry of SH fields); ii) local field distribution at FH; and iii) the propagation properties of SH light inside the holes.
In the calculations, the FH beam is a truncated planewave at λ FH = 830 nm. The source illuminates the system at normal incidence from the air side with the electric field pointing along the x-axis. The intensity used, 0.1 MW/cm 2 , has been estimated from the linear power measurement reported in Fig. 3 of Ref. 32. The whole system is on a glass substrate (n glass = 1.5).
A note about geometry: roughly, the optical response of a RHA at λ FH has localized character for configurations with n glass p << λ c ≈ λ FH . On the other hand, the scattering properties are dominated by the coupling to surfaces modes at larger periods. In general, some degree of hybridization between localized modes and surface waves always exists. 52 In this section, the period is varied and the hole shape is fixed (a x = a y = 280 nm). These geometrical parameters are chosen so that we can explore both optical regimes keeping the incident wavelength unchanged.
A. Symmetry of SH fields: far-field coupling Figure 2 shows the computed electric field on a plane placed at z = 125 nm inside a single hole of an array with p = 500 nm: for (a) λ FH = 830 nm and (b) λ SH = 415 nm. The amplitude is superimposed to the vector field map.
The FH field, which is represented in Fig. 2(a), corresponds to the transversal electric mode of the waveguide TE 0,1 , for these range of parameters and wavelengths typical in EOT. 51 This EM mode is the fundamental one at FH, being its cutoff wavelength λ TE0,1 c = 712 nm. In Fig. 2(c), the electric field distribution of the TE 0,1 mode is obtained analytically for a perfect electric conducting (PEC) infinite waveguide, with the considered cross section. The electric field pattern for PEC looks quite alike to the one for gold. Usually, the identification between EM modes in real and PEC metals is justified. At first approximation, the finite conductivity of real metals translates to an effective enlargement of the hole size due to the EM field penetration into the metal, 53 which explains the deviations between (a) and (c) in Fig. 2.
The symmetry properties of the waveguide modes excited can be characterized by the action of the operator Π j , which provides the parity symmetry of the electric field along the j-direction over the orthogonal plane that crosses the center of the hole.
The parity along the y-direction, Π y , is identical for the modes excited at FH and SH. The y-component of the incident electric field is zero, so it does not impose any constraint through that direction.
The parity along the x-direction, Π x , is different for the FH and SH fields. The x-component (y-component) of the TE 0,1 mode has even (odd) parity, that is, (−x, y), for both gold and PEC. The incident field, a linearly polarized plane wave, only can excite EM modes with that symmetry at normal incidence. The electric field at SH has opposite symmetry through the action of Π x [see Fig. 2 Therefore the SH fields can not couple to the TE 0,1 mode at λ SH , but they might couple to the TE 1,1 and TM 1,1 modes. These are the least decaying modes that can be excited at SH, λ TM1,1 c = 508 nm and λ TE1,1 c = 585 nm. None of them individually leads to the vector field pattern of Fig. 2(b) (numerically calculated at SH). Instead a superposition of both is excited at SH as shown in Fig. 2(d), corresponding to the analytical calculations for a PEC waveguide. The relative amplitude and phase of each mode has been adjusted to reproduce the numerical result. A different aspect ratio or hole shape would end up in a different balance between the corresponding waveguide modes. 38,54 The cutoff wavelengths of the modes excited at SH are larger than λ SH , so light propa-gates through the holes with an amplitude only affected by the absorption in the metal.
The change of symmetry can be readily explained. The mirror symmetry of SH fields results from the properties of the second order polarization vector. Accordingly to the approximated expression P (SH) ≈ χ (2) nnn |E (FH) n | 2 n for Eq. 1, the direction of P (SH) at a given point is approximately determined by the unitary vector normal to the surface at that location. By definition, n is positive at the left side wall of the hole (x = −a x /2), and negative at the right side wall (x = +a x /2). The symmetry of the charge density at SH switches from odd to even because |E (FH) n | 2 is equal at both sides [ Fig. 2(a) and Fig. 2(c)]. Symmetry is crucial in the SHG process and determines the ability of the SH field to couple with the propagating diffraction orders in the substrate and the incident half-space. In a periodic structure only propagation through those directions given by the Bragg's condition are allowed at normal incidence, set by the parallel to the surface reciprocal lattice vector, G i,j = 2π p (i, j) for a square lattice (i and j are integers). Using the dispersion relation of light in the substrate we calculate the angle θ i,j for the (i, j)-order in transmission with respect to the vertical direction, expressed as function of the fundamental wavelength: Propagating modes are then associated to absolute values of sin (θ ij ) equal or less than unity, yielding evanescent modes otherwise. In the reflection side we can calculate the corresponding angles by taking n glass = 1. The situation is different in the case of SHG. Note that the 0th diffraction order field is constant on a given x − y plane, so it has the same parity symmetry that of the incident field at normal incidence. As a result, SH and 0th diffraction order fields have opposite parity symmetry, so the overlap between both EM modes is zero. Therefore, there is no SH radiation in a half-space with all the diffraction orders being evanescent except the 0th diffraction order. In our calculations, SH emission in transmission at the 0th diffraction order, J T SH (0, 0), is sixteen orders of magnitude (at the level of noise due to numerical round off) less intense than through other directions, for all periods. The same occurs in the reflection region, confirming that our FDTD implementation fully respects the symmetry of SH fields. Similar arguments based on symmetry considerations explain why specular radiation is forbidden for regular arrangements of nanoparticles. 55 Moreover SH photons with G (SH) i,j = 0 are forbidden under the normal incidence condition, so only the evanescent waves on the lattice at FH (SPPs, non-propagating near-field scattered by the holes,...) can create SH fields obeying momentum conservation. 56 Interestingly, the different diffraction orders that contribute to the same θ i,j do not need to have the same intensity. As example, we define two SH powers in trans- mission per unit cell for the 1st diffraction order, which are represented in Fig. 3(a) as a function of the period. The parallel component, J T SH (±1, 0) = J T SH (+1, 0)+ J T SH (−1, 0), accounts for the SH fields radiated on two planes defined by vectors G ±1,0 , which are parallel to the electric field of the FH source. The perpendicular component J T SH (0, ±1) = J T SH (0, +1) + J T SH (0, −1) represents the same but for the planes given by the corresponding diffraction orders, which are perpendicular to the electric field of the FH source in this case. Parallel and perpendicular components at SH are different, a result that will be explained later on. The 2nd diffraction order J T SH (±1, ±1), the sum of the four possible combinations of the (±1, ±1) diffraction orders, is represented in Fig. 3(b). In this case, because the reciprocal lattice vectors are parallel to the diagonals of the unit cell each order carries the same SH intensity. Finally, the results for the 3th diffraction order are shown in panel Fig. 3(c). The 3th diffraction order becomes evanescent for periods shorter than 553 nm, as predicted by Eq. 2. Note that the angle of emission increases as the period size decreases, for a given diffraction order (θ is shown for a few periods in Fig. 3).
The emission of SHG into the transmission region splits in different angular contributions, which can be controlled by the period size. On the other hand, the profile of total SH emission in transmission as a function of the period is characterized by a dip followed by a peak, both features related to the excitation at FH of SPPs (on the glass-metal interface). They are determined by the field intensity and distribution at FH. This very influent factor on SHG is discussed next.
B. Field distribution at FH: local source of SH fields In Fig. 4(a) the FH power transmitted per unit cell, J T FH , is calculated for several periods ranging from 610 nm to 400 nm. The dip at p ≈ 540 nm corresponds to an EOT minimum (indicated with a vertical line, as a guide to the eye). The EOT maximum at FH is reached for p ≈ 500 nm. The corresponding SH powers radiated in transmission and reflection per unit cell, J T SH and J R SH , are shown in Fig. 4(b) and Fig. 4(c) with solid circular symbols. Figure 4(b) also shows with empty squares the SH power emitted in transmission calculated by neglecting χ (2) ntt and χ (2) tnt contributions (but keeping χ (2) showing that is an excellent approximation as advanced in the introduction.
Second harmonic generation fundamentally depends on the specific details (phase and intensity) of E (FH) on the metal surface, unlike for example two-photon luminescence that essentially depends on the local intensity. 57 Let us analyze the EM modes supported by the investigated structures, which are ultimately related to the features of the linear transmission spectrum of Fig. 4(a).
The linear transmission spectra for periods p = 560 nm, p = 540 nm and p = 500 nm are shown in panels (a), (d) and (g) of Fig. 5. The vertical lines depict the wavelength of the external source λ FH . Three EOT peaks can be distinguished in each figure, resulting from the excitation of surface EM modes of the corrugated structure. 51 Each surface mode has associated a full EOT feature, characterized by the typical profile of a Fano resonance. 58 At the EOT minimum the SPP of the holey film is hardly affected by the presence of holes, 59 so we can identify every peak in Fig. 5 with the help of the flat surface dispersion relation for SPPs, k SPP . 60 The frequency at which a SPP can be excited at normal incidence, λ SPP , is given by the condition of momentum conservation at the surface and can be approximately calculated by folding k SPP into the first Brillouin's zone, i.e., k SPP = |G i,j | (where G i,j is a reciprocal lattice vector). For example, the EOT minima in Fig. 5(d) (p = 540 nm) are located at wavelengths: 830 nm, 640 nm and 574 nm. These energies correspond to three different bounded EM modes. The EOT peak at the near infrared is due to the (±1, 0)-SPP of the glass-metal interface, being λ SPP = 844 nm. The prominent feature at visible is mainly due to the glassmetal (±1, ±1)-SPP, with λ SPP = 636 nm. Slightly blueshifted, we can also see a less intense and more narrow peak which corresponds to the air-metal (±1, 0)-SPP, with λ SPP = 578 nm.
It is reasonable to think that the field distribution of SPPs supported by a RHA determines which surface in the system is the main source of SH radiation. Our formalism allow us to switch off the generation of SH fields from a given surface (by forcing χ (2) = 0 at that surface). This allows us to obtain insight into the predominant SH process. We have calculated J T SH and J R SH as all SHG would be exclusively generated from three different regions separately: transmission side (output), reflection side (input) and hole walls. The output and input contributions include SHG from the edges and corners of the holes. The radiation in transmission at SH originating from the output interface, J T,out.

SH
, is shown in Fig. 4(b) with empty circular symbols, while the remaining SH radiation (input+walls) is shown with triangles. The corresponding results for the reflection region are shown in Fig. 4(c). Note that the superposition principle   In transmission, the glass-metal (±1, ±1)-SPP is responsible for radiation at SH in transmission from p = 610 nm to p = 550 nm, given that J T SH ≈ J T,out.

SH
. However, the contribution to SH from walls and input surface is not negligible [see triangular symbols in Fig. 4(b)], which is coherent with near-field at SH [Fig. 5(c)]. For periods ranging from p = 540 nm to p = 400 nm J T SH ≈ J T,out.

SH
, so generation at the output region accounts for most of the SH emission in transmission. This result is also coherent with the SH near field maps shown in Fig. 5. On the output surface, the glass-metal (±1, 0)-SPP develops within that period range. This SPP has a clear fingerprint in the FH nearfield at the glass-metal interface, seen in Fig. 5(e) and Fig. 5(h) (see relative scale). At the EOT minimum the near-field is intense enough to generate strong local SH fields [ Fig. 5(f)], which are comparable to those generated at maximum transmission (only six times more intense) [ Fig. 5(i)]. The optical response at SH in transmission has then a straightforward explanation: overall the generation of SH radiation in transmission is controlled by enhanced fields at FH on the glass-metal interface because of the excitation of surface waves bound to that surface, as expected for an asymmetric dielectric configuration 61 .
In the reflection region, the radiation process at SH presents three different regimes [see Fig. 4(c)]. From 610 nm to 550 nm, J R,out.

SH
≤ J R,in.

SH
. The interpretation is that both walls and input side contributions to SH are important in reflection. For that period range, the system is having access to the (±1, ±1)-SPP at FH, which is bound to the glass-metal surface. In the input side, the FH field is enhanced at the holes [ Fig. 5(b)]. This field is characterized by a combination of localized and surface modes. In fact, the near-field within that period range is affected by the optical response of a single hole as reported in Ref. 52, given the close proximity between the FH wavelength and the hole cutoff (≈ 712 nm). Therefore, the FH field is distributed at both sides of the metal layer and inside the holes, yielding a more complex optical response at SH, specially in reflection. From period at which the EOT minimum occurs for the chosen incident wavelength, up to p = 435 nm we find that J R SH ≈ J R,out.

SH
. The FH field is asymmetrically distributed and concentrates on the output surface as it corresponds to the glass-metal (±1, 0)-SPP. To explain this behavior we need to understand how SH light propagates inside the holes, which is discussed in the following paragraphs. Finally, for p < λ SH all off-normal diffraction orders are evanescent in air, so SH radiation in reflection is zero within the round off precision in our numerical calculations.

C. Propagation inside the holes at SH: light absorption
In Sec. III A we advanced that propagation inside the holes is only limited by light absorption in the metal, given that λ SH < λ TM1,1 c < λ TE1,1 c , for the investigated hole dimensions. Within the period range where SHG is caused by the glass-metal (±1, 0)-SPP, the SH fields created at the glass-metal interface can go through the holes and be emitted into the reflection region, which explains that most of the SH emission in reflection originates at the output surface. Moreover, J T SH and J R SH intensities have the same order of magnitude [compare Fig. 4(b) with Fig. 4(c)], except for the shortest periods where J R SH = 0. The amount of SH light in reflection is not negligible, so light absorption is not limiting the SH emission process.
Un punto más en la transición? To further investigate the consequences of light absorption inside the holes we compare gold with a series of hypothetical metals, for which the real part of the dielectric constant is essentially that of silver but the imaginary part, ε i , is modified. The value of χ (2) in these simulations is the one of gold. The SH powers in transmission and in reflection over total SHG are shown in Fig. 6, as a function of the imaginary part of the dielectric constant at SH. We have chosen p = 500 nm, which illustrates the optical response of the system at the glass-metal (±1, 0)-SPP. The skin depth, the cutoff wavelengths of the waveguide modes at SH, and the near field pattern at FH are practically the same in all the series. We have also found that J T SH ≈ J T,out.
SH and J R SH ≈ J R,out.

SH
for all ε i considered (not shown). In the lossless case J R SH > J T SH . The difference between both of them reduces with increasing absorption inside the holes. Eventually the inequality in reversed, and J T SH > J R SH , providing a clear evidence of the critical role played by metal absorption.

A. SPP related effects
In this section, the linear and the nonlinear response around the glass-metal (±1, 0)-SPP is analyzed in more detail. We define two optical properties to characterize emission: SH efficiency in transmission and the change of polarization of the first diffraction order.
In Fig. 7(a) we show the results for SH efficiency in transmission, defined as: This coefficient is independent of the illumination intensity at FH. The vertical line, indicating the period for which the wavelength of the considered incident beam λ FH corresponds to an EOT minimum at FH, coincides with the maximum value of α. Efficiency at the EOT minimum is five orders of magnitude larger than the minimum value obtained among the periods investigated. The FH field distribution is again the key point. The optical properties of a SPP in a holey metal film are different at different wavelengths. For convenience, we distinguish between the response at the EOT minimum and the response at the transmission maximum using different labels in each case. At the EOT minimum the SPP of the holey surface remains "unperturbed", showing the same optical response than that of a SPP on a flat metal surface, as explained in Sec. III B. Therefore, the influence of holes can be ignored assuming that the EM field at their centers is negligible. A mode like this is weakly coupled to the far-field, so we call it dark-SPP. The SPP of the corrugated structure is highly "perturbed" at the EOT maximum, so its dispersion relation deviates from the one of a flat metal surface. The holes scatter light and we call it bright-SPP. Efficiency reaches the highest value at the EOT minimum because the dark-SPP is weakly coupled to the far-field in transmission at FH [ Fig. 4(a)], while its near-field at the output surface creates enough SH photons in the transmission region [ Fig. 4(b)]. Figure 3(a) shows that the far-field pattern of SHG is not symmetrically distributed, being different the so-called parallel and perpendicular components of SH radiation for the 1st diffraction order [J T SH (±1, 0) = J T SH (0, ±1)]. We define the change of polarization for the 1st diffraction order as: shown in Fig. 7(b). The parallel component dominates the far-field emission for linearly generated light at the same wavelength of SH [J T FH (0, ±1)/J T FH (±1, 0) ≈ 0]. We could expect φ SH 1 ≈ 0, however φ SH 1 reaches values even higher than the unity. We observe that J T SH (±1, 0) dominates for short periods, while for large periods the balance reverses. In between, a sudden change in the SHG spatial distribution exists, at the period for which λ FH approximately matches the EOT minimum.
As explained, the coupling with surface modes at FH is very weak for short periods. The radiative and even the non-propagating field polarizaton follows that of the incident field, which is the typical response found in isolated holes. 62 As the period increases, the bright-SPP starts to participate in the generation process. Its wavevector is directed along the x-direction, however this plasmon mode is efficiently scattered by the holes and may provoke, for every interaction with them, a strong depolarization of the evanescent FH field thus increasing φ SH 1 . The same argument can be sustained for the largest periods investigated, for which the glass-metal (±1, ±1)-SPP is excited at FH. A singular behavior occurs at the EOT minimum, at the dark-SPP. This mode weakly scatters light at FH, so SH radiation can be only generated along the direction of its wavevector, i.e., the x-axis.
Similar results are expected in experimental setups by tuning λ FH , for a fixed period. In that situation, a measure of the change of polarization might be useful for surface assays, where the balance between the radiated intensity through different directions would determine the quality of a given sample. Surface defects, shape imperfections, debris on the surface, and whatever other experimental realization far from a perfect hole array would have an impairing impact on φ SH 1 . In addition, the EOT minimum (which in simulations appears as a sharp and narrow dip) is quite sensitive to sample imperfections and the size of the system. 63 Therefore, SH efficiency wold provide additional information on high quality structure characterization. 16 Deviations from the expected α profile [ Fig. 7(a)], like broadening, would indicate a mis-alignment of the FH beam, structure imperfections generated during the fabrication process or the presence of chemical products on the surface.

B. Localized resonances and related effects
The influence of localized resonances in SHG from RHAs was first studied in Ref. 32. In this work enhanced SH emission occurring close to the cutoff wavelength of the FH field was reported. That enhancement was explained in terms of slow EM modes localized in the holes at FH. But the conclusion has been challenged by the same authors, 42 and recent experiments have shown that the time delay at FH is similar for different aspect ratios. 40 Here, we compare our numerical results with the experimental measurements and discuss these results at the light of our findings.
In that experimental work several hole arrays were investigated, each consisting of 20x20 rectangular holes milled in a square lattice. The system was deposited on a glass substrate. The film thickness, period and hole area were kept fixed at h = 160 nm, p = 410 nm and S = 3.4x10 4 nm 2 respectively. The period was precisely chosen to avoid hybridization between localized modes and surface waves at FH. 52 The different arrays were designed with holes of different dimensions, characterized by the aspect ratio, AR. The system was illuminated at normal incidence from the air side, being λ FH = 830 nm and with FH peak powers ranging from approximately 1.0 mW to 50.0 mW (intensities from 1.5 KW/cm 2 to 0.075 MW/cm 2 ). We use the same parameters than in the experimental work, but our x-y axes are rotated 90 degrees with respect to Ref. 32. With our choice, AR = a y /a x , a x = S/AR and a y = √ AR x S. Like in previous calculations, the FH beam is a truncated plane-wave at λ FH = 830 nm. The FH source illuminates the system at normal incidence from the air side with the electric field pointing along the x-axis and delivers 0.1 MW/cm 2 . We present in Fig. 8(a) the computed power emitted in transmission at FH (square symbols) and SH (full circular symbols), as a function of the aspect ratio. From the power per unit cell calculated with FDTD, we can directly compare our results with the experimental ones, taking into account that the samples covered 20x20 unit cells. The simulations reproduce the general trend of experiments [ Fig. 3 in Ref. 32]. For instance, both linear and nonlinear transmissions peak at AR ≈ 2 with similar SHG power values. In contrast, the measured linear transmission is flatter as a function of AR, and our result on SHG is characterized by a broader peak compared with the experimental case. In any case, the main features of SHG in RHAs are captured by our numerical implementation.
This system can be analyzed using the ingredients previously presented in this paper. First, linear and non-   linear peak locations coincide like in Fig. 4, where enhanced SH emission occurs at the maximum of linear transmission. In Fig. 8(a), the SHG peak is due to a localized resonance 64-66 excited at the FH wavelength (for an example, see the inset of Fig. 8(a), which represents the FH transmission spectrum at AR=2). At resonance, EM fields are bound to the structure for time scales long enough to generate EM field accumulation at FH. Resonances at SH might not be discarded, however we have not observed any related effect either for surface waves or in localized resonances. Second, SH emission in the reflection region (air side) is forbidden for the chosen period, as expected from Eq. 2. Third, the FH field is more intense at the output surface than at the input surface because of the presence of the dielectric. A mode tends to concentrate its EM energy in regions of high refractive index 67 as shown in Fig. 8(b) for AR=2, producing both strong SH fields at the output [ Fig. 8(c)] and high SH emission in transmission. Again, we can analyze SH radiation from three different regions separately: transmission side (output), reflection side (input) and hole walls. The radiation in transmission at SH originating from the output interface, J T,out.

SH
, is shown in Fig. 8(a) with empty circular symbols, while the remaining SH radiation (input+walls) is shown with triangles. As expected analyzing the SH near-field J T SH ≈ J T,out.

SH
for all aspect ratios. Finally, the corresponding SH efficiency and φ SH 1 is shown in Fig. 9, panels (a) and (b), as a function of the aspect ratio. Efficiency does not show the abrupt change observed in Fig. 7(a). The change of polarization displays two different regimes, below and beyond AR ≈ 1. For AR < 1 the short side of the hole is perpendicular to the incident electric field. In this case we expect a strong change in the polarization of light (from x to y polarization) even at FH. For AR > 1 the long side of the rectangles is now perpendicular to the incident electric field, so there is less depolarization. Having in mind the results of Sec. III, we realize that the period is of utmost relevance compared to other geometrical parameters in SHG from RHAs.

V. CONCLUSIONS
In conclusion, we have investigated SH radiation in periodic arrays of rectangular holes drilled in a metal film. We have conducted FDTD calculations to compare with existing experimental works and investigated the effect of localized resonances in SHG. Furthermore, we have found that SH efficiency in transmission can be enhanced by coupling the FH beam to dark-SPPs associated to EOT minima. In that case, even though FH radiation is small, the FH local field is intense enough at the the exit surface of the array. We have demonstrated how the far-field distribution of SH radiation changes also when the dark-SPP is excited. All these findings can be explained through a subtle mechanism. The FH near-field, the source of SH induced currents, has opposite parity symmetry that of the SH field. Therefore, SHG provides access to EM modes that cannot be excited by the fundamental field. At the end, it is the character of such modes (absorption, cutoff wavelength, overlap with the lattice diffraction orders) what determines whether SH light can be emitted to the far-field or not. We believe our findings will be useful as tool for checking the quality of holey metal systems and probing the near-field at FH through non-invasive methods.