We investigate second harmonic generation (SHG) from hexagonal periodic arrays of triangular nano-holes of aluminum using a self-consistent methodology based on the hydrodynamics-Maxwell–Bloch approach. It is shown that angular polarization patterns of the far-field second harmonic response abide to threefold symmetry constraints on tensors. When a molecular layer is added to the system and its parameters are adjusted to achieve the strong coupling regime between a localized plasmon mode and molecular excitons, Rabi splitting is observed from the occurrence of both single- and two-photon transition peaks within the SHG power spectrum. It is argued that the splitting observed for both transitions results from direct two-photon transitions between lower and upper polaritonic states of the strongly coupled system. This interpretation can be accounted by a tailored three-level quantum model, with results in agreement with the unbiased numerical approach. Our results suggest that the hybrid states formed in strongly coupled systems directly contribute to the nonlinear dynamics. This opens new directions in designing THz sources and nonlinear frequency converters.

## INTRODUCTION

The research in plasmonics found various applications in chemistry,^{1,2} biology,^{3} and engineering.^{4,5} All are due to the ability of plasmonic materials to sustain plasmon–polariton resonances with a very small mode volume.^{6} Even though characteristic Q-factors of such modes are relatively small, large local field enhancement^{7} makes such systems very attractive for numerous applications, including nonlinear nano-optics^{8,9} and plasmon–polariton chemistry.^{10} The latter combines small plasmon mode volumes with molecular excitons.^{11} When the coupling strength of a local field and an exciton surpasses all the damping rates in a system, the strong coupling regime is reached. In such a regime, new hybrid states called upper and lower polaritons are formed.^{12} Under strong coupling conditions, these states have a mixed plasmon–exciton character.

While the linear optics of such systems is being actively explored, the nonlinear phenomena such as second harmonic generation (SHG) recently began attracting considerable attention.^{13–17} What would make nonlinear plasmonics an exciting playground in terms of fundamental as well as applied objectives? From a fundamental viewpoint, interactions of organic and inorganic systems, including molecules of various complexity with metals, challenges our abilities to properly describe such systems using conventional semi-classical methodology. At the application side, advanced sensing and bio-sensing devices relying on the sensitivity of nonlinear effects to perturbations and on the unlimited possibilities for targeted surface sensitization open-up new horizons in diagnosis and therapy.

In this manuscript, we report theoretical investigations of the nonlinear properties of a strongly coupled exciton–plasmon system. It is comprised of a planar array of metallic nanoscale cavities with a contacting molecular layer of a general nature. Whereas SHG has been chosen as the central target at this stage, other phenomena can be predicted along the same theoretical lines, such as other three- and four-wave mixing schemes, among which are third-harmonic generation, the intensity dependent index of refraction, and parametric effects including near IR and THz generation.

Our approach is based on the numerical exploitation of a non-perturbative hydrodynamic-Maxwell model. It can be therefore viewed to some extent as a “numerical experiment,” then calling, very much as in the case of an actual experiment, for simple physical interpretations of numerical findings. Interpretations of the numerical data are based on a semi-classical quantum picture relying on a minimal three-level system that appears sufficient to account for the main features obtained from simulations. The proposed analytical model allows for simple expressions that are prone to straightforward interpretations in terms of photon mediated coupling between the basis states. The purpose of our studies is to theoretically explore in a relatively generic way the potential for nonlinear optics of mixed metal–molecule structures where strong interactions between both species are expected to display new features with the potential to enhance their efficiency based on well identified parameters.

Early theoretical works in nonlinear plasmonics combined either a nonlinear material with otherwise linear plasmonic interface,^{18} employed *ab initio* calculations limited due to obvious reasons to small sizes,^{19} or utilized a semi-classical hydrodynamic model for conduction electrons.^{20} The latter has recently proved to be highly effective even on a quantitative level.^{21} It was recently demonstrated^{15} that by combining the hydrodynamic model for metal with Bloch equations for molecules in the time domain, one can truly explore nonlinear optics of strongly coupled systems. The major question that we attempt to address in this work is as follows: in a strongly coupled system with two well-defined polaritonic states, can we expect to observe the hybrid states directly participating in nonlinear dynamics? And if so, are there direct lower-to-upper polariton transitions? The positive answer to this question would lead to a large body of exciting applications ranging from THz generation to electromagnetic (EM) induced transparency and slow light.^{22}

In the section titled Theory, we overview the main features of the hydrodynamic-Maxwell–Bloch model. Next, we discuss our main results based on this model in terms of linear and nonlinear spectroscopic as well as polarization properties. Finally, we propose a simple analytical three-level quantum model, confront it to the former “numerical experiment,” and then draw conclusions. Additionally, we provide the Appendix, which discusses the polarization properties of SHG based on abstract tensor and symmetry considerations as well as a more visual geometric argumentation.

## THEORY

We consider periodic arrays of triangular nanoholes arranged in a hexagonal lattice, as schematically depicted in Fig. 1(a). The optical response of this system is evaluated using Maxwell’s equations in the space–time domain,

where the macroscopic polarization **P**_{pl} is determined by the hydrodynamic model of conduction electrons for a given plasmonic system [our notation with the subscript “pl” (plasmon) is to distinguish between polarizations in the plasmonic system and the molecular layer as discussed below].^{23} The equation of motion describing **P**_{pl} combines both the continuity equation relating the polarization current and the number density, *n*, and the equation determining dynamics of the electron velocity field.^{24} Such an equation can be written in a concise form as^{25,26}

where $me*$ is the effective mass of a conduction electron, *n*_{0} is the equilibrium number density of electrons, and *γ*_{e} is the phenomenological damping rate. In the simulations, we parameterize the plasma frequency, $\Omega p=n0e2/\epsilon 0me*$, the number density, and the damping rate. In this paper, we consider aluminum with the following parameters: Ω_{p} = 10.85 eV, *γ*_{e} = 0.4641 eV, and *n*_{0} = 1.81 × 10^{29} m^{−3}. The geometrical parameters defining the array are shown in the caption of Fig. 1. A few general remarks regarding the hydrodynamic model employed are in order. We first note that this approach is not perturbative and thus, in principle, contains all possible nonlinear phenomena. The terms leading to a nonlinear dependence of the optical response on the field intensity are the Lorentz magnetic term, $P\u0307pl\xd7B$, the quadrupole term, $E\u2207\u22c5Ppl$, and the convective terms [last two terms in Eq. (2)]. All the nonlinear terms were shown to be important in SHG simulations of periodic arrays of nanoholes.^{15} It was also demonstrated^{24} that the second order perturbation theory applied to Eqs. (1) and (2) can be used to quantitatively describe the SHG process in plasmonic systems.^{25} We also note that the hydrodynamic model coupled to Maxwell’s equations results in the SHG signal at flat metal interfaces. This is because the incident k-vector breaks the symmetry leading to SHG. All these features make the hydrodynamic model superior compared to approaches that rely on considering the second-order nonlinear polarization using phenomenological second-order susceptibility tensor.^{21}

In addition to the metal, we also consider two-level molecules forming a thin layer on the input side of the array. The molecular layer is simulated using rate equations^{12} that determine the dynamics of the number density of molecules in the excited state, *N*_{2}, and the corresponding macroscopic polarization of the molecular layer, **P**_{mol},

where *γ*_{21} is the population relaxation rate, Γ = *γ*_{21} + 2*γ*_{d}, *γ*_{d} is the pure dephasing rate, *ω*_{12} and *μ*_{12} are the transition frequency and the transition dipole, respectively, and *N*_{0} is the total number density of molecules. The initial conditions correspond to all molecules in the ground state. The parameters used in this work are *γ*_{21} = 4.14 meV, *γ*_{d} = 7.03 meV, *μ*_{12} = 10 D, and *N*_{0} = 2 × 10^{26} m^{−3}, and the thickness of the layer is 20 nm. In addition, in Eq. (3), we assumed that the total number of molecules is conserved, *N*_{0} = *N*_{1} + *N*_{2}.

Equations (1)–(3) are discretized in space and time following the finite-difference time-domain method.^{15} When evaluating both linear and nonlinear responses, the system is excited by a planewave at normal incidence. We follow the conventional total field/scattered field approach^{27} to generate plane waves. Periodic boundary conditions are applied along the X and Y axes [see Fig. 1(a)], and convolution perfectly matched layers^{28} are set on top and bottom of the grid to absorb outgoing waves. The EM energy flux is detected on the input (reflection) and output (transmission) sides of the array in the far-field zone as the corresponding Poynting vector integrated over a given detection plane. In order to speed the simulations up, we employ a three-dimensional domain decomposition technique and message-passing interface subroutines. Our home-build codes are parallelized for 12 × 12 × 8 = 1152 processors with a spatial resolution of 1.5 nm and a time step of 2.5 attoseconds. Linear spectra shown in Figs. 1 and 3 are usually converged after propagating equations for 300 fs. The average execution times vary between 15 and 30 min. The nonlinear simulations are performed using a 500 fs long pump. To eliminate spurious nonlinear contributions that always occur due to spatial discretization, we apply the Blackman–Harris time window to all detected fields.^{19} Such a procedure results in clear power spectra with well-resolved harmonics. The execution times of our nonlinear codes are between 2 and 3.5 h. The nonlinear response is calculated the following way. The system is driven by a 500 fs strong laser pulse at a fixed pump frequency and the peak amplitude of 5 × 10^{7} V/m. The latter is chosen such that the intensity of the second harmonic scales quadratically with the pump intensity. The signal is calculated as an outgoing EM energy flux integrated over detection planes on the input (reflection) and output (transmission) sides of the system. At each grid point along the detection planes, we compute corresponding EM field amplitudes as functions of time and space. Once the time propagation is complete, we perform Fast Fourier Transform (FFT) of all detected EM field components. The spatial integration of the corresponding Poynting vector is performed after FFT. The obtained signal is obtained as a function of the outgoing frequency. The SHG signal is extracted from the total signal by integrating over outgoing frequency in the vicinity of the second harmonic of the pump.

## RESULTS AND DISCUSSION

Figures 1(b) and 1(c) present linear transmission and reflection spectra of the nanohole array in vacuum for different thicknesses of the metal film. We note that both transmission and reflection exhibit multiple resonances with the broad mode clearly demonstrating the splitting, which is more pronounced in the reflection. The splitting is gradually disappearing with increasing thickness indicating that the physics of this is due to the coupling between surface modes on the input and output sides of the array, forming symmetric and antisymmetric states with respect to the center of the film (such modes were first discussed in Ref. 29 and then later in Ref. 30, and experimental implications of these modes and their properties were examined in Ref. 31). The broad mode with the splitting centered at 2.45 eV corresponds to a localized plasmon mode with an EM field primarily localized near the vertices of holes. The other modes observed near 3 eV are Bragg-plasmon modes. In this paper, we consider the localized plasmon mode, its coupling to molecules, and its influence on SHG. We note that the localized plasmon enhances SHG efficiency notably higher compared to Bragg plasmons, as demonstrated in Refs. 15 and 17.

Figure 2 shows angular diagrams of SHG response detected in the transmitted signal. We perform two types of detection: (1) calculating the SHG signal polarized along the pump polarization and perpendicular to it [Figs. 2(a) and 2(b)] and (2) calculating the SHG signal polarized along the X and Y axes in the fixed (laboratory) system of coordinates [see Fig. 1(a) where the X and Y axes are defined].

Figure 2 displays two strikingly different SHG polarization patterns. Panels (a) and (b) exhibit six lobes in a sixfold symmetry configuration, whereas panels (c) and (d) exhibit four lobes in a fourfold symmetry configuration. These angular symmetry features are discussed in great detail in the Appendix based on the tensor symmetry considerations along with a pictorial discussion of the non-intuitive fourfold emission from a threefold symmetric nonlinear response.

It is informative to examine the SHG signal as a function of the pump frequency. In general, any dispersive material with a second order nonlinearity exhibits the second-order susceptibility, which can be written as a product of the first-order susceptibility at the fundamental and second harmonic frequencies following Miller’s rule,^{32}

We thus should expect to observe peaks in the SHG signal corresponding to one- and two-photon emission. Figure 3 presents the conversion efficiency in transmission and reflection as a function of the pump frequency. The conversion efficiency defined as the ratio of the intensity of the second harmonic signal, $I2\omega $, and the intensity at the fundamental frequency, $I\omega $.

The results presented in Fig. 3 are compatible with Miller’s rule for plasmonic systems. There are clear signatures of both one- and two-photon emission enhanced by corresponding plasmon modes of the array. As expected, the largest enhancement of the conversion efficiency is due to the localized plasmon at 2.45 eV, which is higher than that due to the Bragg plasmon at 3.1 eV by more than an order of magnitude as was reported elsewhere.^{15} The two-photon emission peaks also have an obvious correspondence to the plasmon modes of the array when compared with the red dashed line. Interestingly, Wood’s anomaly, i.e., the peak associated with the period of the array,^{33} influences the conversion efficiency as seen near 2 eV as a knee shape of the conversion efficiency. Another observation worth emphasizing is the splitting of the plasmon peak discussed in Figs. 1(b) and 1(c), which is mostly pronounced in the reflection near 2.4 eV. One can note that the two-photon emission peak at 1.2 eV also evidences a clear signature of the splitting. Interestingly, the splitting due to the hybridization of the surface modes on the input and output sides of the array has a different structure in the one-photon and two-photon peaks. Note that the low frequency mode has a higher conversion efficiency at the one-photon peak compared to the higher frequency mode. The two-photon peak has this reversed. This again reinforces the hypothesis that the local symmetry plays an important role in the SHG process.

We now add a thin molecular layer placing it on top of the array and examine its influence on the SHG. Figure 4 explores the linear optical response of such a system. The transmission with molecules being resonant with the localized plasmon mode exhibits a clear Rabi splitting indicative of strong coupling between molecular excitons and the corresponding plasmon mode. Figure 4(b) shows frequencies of the upper and lower polaritons as functions of the periodicity of the array exhibiting expected avoided crossing. The Rabi splitting observed for the array with a period of 395 nm is measured to be 179 meV. Next, we explore the SHG process in such a system under strong coupling conditions.

The main result of this manuscript is shown in Fig. 5. Here, we plot the SHG efficiency for the molecule/metal system under strong coupling conditions and compare it directly with the results obtained without molecules. Unexpectedly, the two-photon peak exhibits a clear signature of the Rabi splitting. Previously, it was argued that the Rabi splitting observed in the second harmonic signal for strongly coupled systems can be explained by expanding the two coupled oscillator model.^{14,16} *With molecules treated as simple two-level emitters with no diagonal elements of the transition dipole, the only source of the second harmonic is the metal*. The model based on a strongly driven oscillator with a second order nonlinearity (plasmons) coupled to a linear oscillator (molecular excitons) results in the modified second-order susceptibility (5).^{14} Only the second term in (5) carries the signature of the strong coupling,

where *g* is the coupling strength between the plasmon (SPP) and the corresponding molecular exciton (exc). Equation (5) evidently predicts that the Rabi splitting does not appear in the two-photon peak since only plasmon enhancement is observed.

In contrast with the coupled oscillator model, our results shown in Fig. 5(a) indicate that the molecules also contribute to the two-photon emission. One may argue that the splitting seen near 1.2 eV may correspond to the absorption of two 1.2 eV photons by the molecules. However, the splitting observed at 1.2 eV equals exactly the half of the Rabi splitting, while the absorption linewidth for molecules is significantly smaller. This is also illustrated in Figs. 5(c) and 5(d), where we plot the linear absorption of the molecular layer right next to the corresponding one- and two-photon peaks. One can see that the linear absorption is significantly narrower and notably red-shifted compared to the observed resonances. Another interpretation of the observed splitting at 1.2 eV is thus needed to account for the 2*ω* splitting.

We propose a three-level quantum model in the semi-classical approximation to describe a two band plasmon–polariton spectrum close to the avoided crossing point. 0 designates the ground state, 1 denotes the lower polaritonic branch with energy *ℏω*_{1}, and 2 denotes the upper one with energy *ℏω*_{2}. We limit ourselves to the component $\chi yyy2$ for threefold symmetry as the only other component $\chi yxx2$ can be obtained via $\chi yyy2+\chi yxx2=0$. We also assume that the transitions are only weakly populating the excited states, making the population dominantly in the ground state. In the density matrix formulation,^{32} the only remaining diagonal term is $\rho 000=1$, thus neglecting excited state populations. In the following, we shall omit the Cartesian indices, which will all be implicitly taken as *yyy*.

We can now proceed to adopt the general expressions given [page 174, Eq. (3.6.16) in Ref. 32] in the case of a non-dipolar three-level system. For an *N*-level system, we have

where each term in *a* and *b* corresponds to a three-photon transition from the ground state and back, with any pair of intermediate excited states chosen among the *N* states.

Reduction to a three-level system limits the summation to *n* = 1, 2, likewise for *m*. Moreover, threefold symmetry and the absence of dipolar properties in the molecular layer lead to the cancellation of all participating dipoles in the summation, resulting for the case of a three-level system in *μ*_{00} = *μ*_{11} = *μ*_{22} = 0. The *a* and *b* terms take then the following simplified expressions:

where *ω*_{12} = *ω*_{1} − *ω*_{2} < 0 and *ω*_{21} = *ω*_{2} − *ω*_{1} > 0.

The two denominators of the term clearly exhibit two- and one-photon resonances, respectively, at *ω*_{1} and *ω*_{2} in agreement with the *a priori* non-perturbative fully vectorial hydrodynamic-Maxwell–Bloch model. This behavior corresponds to one- and two-photon resonant second harmonic generation processes for a pair of two-level systems, namely, $0,1$ and $0,2$. Using simplified Eq. (7), we performed fitting of the second-order susceptibility to numerical results. The result is shown in Fig. 5(b). Although we note that the fitting is obviously not perfect, the model adequately describes the splitting observed at one- and two-photon peaks.

The four possible corresponding schemes are sketched in Fig. 6 under (*a*_{1,2}) and (*b*_{1,2}). Following conventions in Fig. 6, $\omega a1$ and $\omega b1$ are in one-photon resonance with levels $1$ and $2$, whereas $2\omega a2$ and $2\omega b2$ are in two-photon resonance with the same levels. In contrast with expression *a*, the two denominators in *b* exhibit only one-photon resonances, respectively, at both *ω*_{1} and *ω*_{2}.

The b term does not exhibit resonance in the optical regime since *ω*_{21} ≪ *ω*, 2*ω*. However, a modified b term could be of interest in the context of a possible difference-frequency generation (*DFG*) process $\omega \u2032\u2192\omega ,\Omega =\omega \u2032\u2212\omega $, inverse to the corresponding $\omega ,\Omega \u2192\omega +\Omega $ sum-frequency generation (*SFG*) process, whereby a coherent beam at frequency Ω in the THz regime could be generated in resonance with the *ω*_{21} transition. The corresponding susceptibility can be evaluated in two steps: the first one by way of generalizing the *SHG* susceptibility as detailed below onto *SFG* $\omega p,\omega q\u2192\omega p+\omega q$ and then replacing *ω*_{p} by −*ω*_{p} corresponding to *DFG*.

For *SFG*, we have

For *DFG*, we obtain

We further chose the incoming frequencies close to each other such that the frequency difference *ω*_{q} − *ω*_{p} be much smaller than optical frequencies, namely, *ω*_{q} − *ω*_{p} ≪ *ω*_{1}, *ω*_{2} and of the order of the frequency shift between the higher and lower polaritonic branches, i.e., *ω*_{q} − *ω*_{p} ≃ *ω*_{21}, which can span the THz domain.

Under such conditions, the *α*_{1} expression cannot support any resonances at optical frequencies as it would otherwise contradict the above assumptions. The *α*_{2} expression can be resonant, respectively, at *ω*_{q} = *ω*_{1} for its first term or at *ω*_{q} = *ω*_{2} for the second one.

More interesting contributions are arising from the *β*_{1} and *β*_{2} expressions, i.e., resonances at *ω*_{q} − *ω*_{p} = *ω*_{21}. In addition, *ω*_{q} = *ω*_{2} in the second term of *β*_{2} can lead to a second resonance. Such a doubly resonant scheme is displayed under (c) in Fig. 6. Whereas such a scheme would be probably of limited interest in a bulk medium with the benefit of a double resonance offset by the cost of corresponding absorptions for each resonance, the situation appears much more favorable for surface driven nonlinearities at the nanoscale where propagative absorption is negligible at the scale of the structure.

We note that the doubly resonant process described above could be of major interest for the efficient generation as well as detection of THz radiation, the structure studied here playing the role of a nanoscale emitter–receptor antenna in the THz range.

At a more fundamental level, such a process would also be of major spectroscopic interest as it could permit to resonantly access the otherwise elusive *μ*_{12} transition dipole that couples levels 1 and 2. In the generic octupolar molecular system described in Ref. 34, levels $1$ and $2$ then referred to as $1\u0304$ are degenerate for essential symmetry reasons resulting from threefold symmetry. Our present plasmonic system is mimicking such a system with the advantage of providing a handle for tuning a strong field induced split between $1\u0304$ and $1$.

## CONCLUSION

Utilizing the non-perturbative hydrodynamic-Maxwell–Bloch approach, we examined the properties of the SHG process in a system comprised of two-level emitters strongly coupled to a localized plasmon mode supported by a periodic array of nano-holes in a thin metal film. It is shown that the effective second-order susceptibility exhibits the Rabi splitting at both one- and two-photon peaks. We proposed a simple analytical three-level quantum model, which allows setting on a firm basis the following main findings:

The Rabi splitting observed for both one- and two-photon transitions provides the evidence for strong field interaction between plasmonic and molecular parts of the system. Whereas the occurrence of such splitting and its sensitivity to the lattice period have been observed and reported for linear properties, it is to our knowledge a first prediction in the context of nonlinear properties. Such features are recognized in both approaches at comparable spectral positioning and intensity levels.

A set of four resonances corresponding to one- and two-photon driven transitions from the ground level to the two plasmon–polariton branches is evidenced in the simulated SHG spectrum for both models.

Polarization properties from the simulation properties show well-defined threefold and fourfold patterns depending on the experimental consideration. These are accounted for in simple algebraic and geometric terms in the Appendix.

Our theoretical studies are surely not to be seen by any means as replacement for actual experiments currently under way but can be considered of a high predictive value in terms of interesting and possibly unexpected phenomena well worthy of a subsequent experimental endeavor. Our results are to be seen as a preliminary step that paves the way onto the design, nanofabrication, and optical testing of original configurations of considerable fundamental and applied interests.

## ACKNOWLEDGMENTS

This work was sponsored by the Air Force Office of Scientific Research under Grant No. FA9550-19-1-0009. The authors also acknowledge computational support from the Department of Defense High Performance Computing Modernization Program. A.S. acknowledges the ISF (Grant No. 2525/17).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

### APPENDIX: THE FOURFOLD SYMMETRY PATTERNS IN SHG ANGULAR DIAGRAMS

We proceed to clarify here to account for the patterns of SHG angular polarization plots (APPs) from nonlinear entities with threefold symmetry.^{35} Such patterns may exhibit fourfold and sixfold, depending on the experimental configuration used, namely, the respective orientations of the incoming *ω* polarizer vs the 2*ω* outgoing analyzer, together with the sample rotation scheme, either rotating in a fixed polarizer–analyzer configuration or conversely with a fixed sample in a rotating polarizer–analyzer configuration. As expected from adequate overall symmetry considerations that encompass both the intrinsic symmetry of the microscopic nonlinear entity at stake together with that of the measurement configuration at the macroscopic laboratory scale, the former microscopic symmetry can be lowered or modified by its convolution with the macroscopic level. Such a response has been already observed in many instances from optically poled polymers^{36} to molecular crystals^{37} or semiconductor nano-crystals with tetrahedral symmetry^{38} onto metallic holes or metallic nano-particles with triangular shapes.^{14,39,40}

We shall further explain here the counter-intuitive difference between the threefold symmetry shape of the nonlinear entities and the fourfold symmetry of the resulting SHG angular polarization plot in relation with the chosen polarizer–sample–analyzer (PSA) experimental configuration. A different setup configuration in which the sample is rotating leads to a more intuitive threefold symmetric APP polar plot that directly reflects the symmetry of the sample.

We shall consider in the following the two configurations depicted in Fig. 7 [note that Fig. 7(b) has been demonstrated elsewhere^{41,42}] and provide a tensor-based justification for their APPs as well as a more physically intuitive rationale. Our ultra-thin sample can be considered to a good approximation as quasi two-dimensional thereby canceling all tensor coefficients with a *z* out-of-plane index, thereby limiting the number of non-vanishing tensor coefficients to that with coefficients limited to combinations of the sole *x* and *y* indices. Under this approximation, only eight non-zero tensor coefficients survive, out of which four are independent coefficients of $\chi ijk2$ as a result of index permutation symmetry

In the case of isosceles triangular symmetry, an additional mirror plane symmetry perpendicular to the symmetry axis of the triangle leads to the cancellation of additional tensor elements that are anti-symmetric with respect to the mirror plane, e.g., featuring an odd number of *x* or *y* coefficients, namely,

This further reduces the number of non-zero coefficients from eight to four coefficients, out of which only two are independent,

A further increase in symmetry elements from that of isosceles triangles to equilateral triangles, that is, from *D*_{2h} to *D*_{3h} imposes an additional relation that results from the very definition of octupolar symmetry, namely, the cancellation of all dipolar tensors attached to the entity of interest. Such symmetry induced cancellation applies to proper physical vectors such as their permanent dipole moment or adequate combinations of higher order tensor coefficients behaving as vectors under point-wise rotations. In the case of a three rank-order tensor such as $\chi ijk2$ of interest here, the resulting $\chi \u20d7i2$ irreducible vector component of the full $\chi ijk2$ tensor is given by $\chi \u20d7i2=\chi ijj2$ (adopting Einstein’s convention for repeated index summation). This leads here to $\chi ijj2=0$ = 0 for all cartesian indices *I* = *x*, *y*, *z*. Implementation of this condition in the case of a planar threefold object leads to

This leaves only four mutually dependent non-zero coefficients, which all relate to a single independent $\chi \u22252$ coefficient, the ensuing relations between the transverse incoming field, $E\u20d7$, and the induced polarization components being

The projection operation along the unit vectors $e\u20d7x$ and $e\u20d7y$ is the mathematical expression of our experimental configuration (a) in Fig. 7, whereby the output harmonic polarization analyzers are set along those same $e\u20d7x$ and $e\u20d7y$ orientations. This is an important remark that will further account for the fourfold symmetry of the measured APPs that convolves the symmetry of the nonlinear entity proper, that is, threefold, together with that for the chosen instrumental PSA configuration.

Taking a linearly polarized input beam with a transverse fundamental field making an angle *α* with the X axis, with components *E*_{x} = *E*_{0} cos *α*, where *E*_{0} stands for the field modulus, the second order harmonic polarization is then given by

The far-field SHG intensity is, to a multiplicative constant, the above squared expressions of the polarization components

These harmonic intensities reported here in Fig. 7 correspond to experimental configurations with a rotating input linear polarizer on a rotated mount allowing it to cover a full circular rotation in the 0 ≤ *α* ≤ 2*π* interval and two fixed linear harmonic analyzers along the X and Y axes, with the sample itself being set on a fixed mount. These expressions are in agreement with the experimentally observed as well as calculated fourfold symmetry pattern in Figs. 2(c) and 2(d) as from their respective $sin22\alpha $ and $cos22\alpha $ angular dependences. It also reflects the observed $\pi /4$ angular shift inbetween the two *x*- and *y*-analyzed APPs, namely,

is a simple consequence. Figure 8 shows specific configurations that illustrate the ensuing fourfold pattern of the APP for this configuration.

One can further illustrate the $\pi /4$ shift when the analyzer is set along the X axis, with respect to the former configuration, as shown in Fig. 9.

For the sake of illustrating the contribution of the PSA configuration onto the symmetry pattern of the APP in that configuration, let us now investigate an alternative one where the linearly polarized input polarizer and output analyzer are parallel and co-rotating in the 0 ≤ *α* ≤ 2*π* interval, as shown in Fig. 10, whereby the $\alpha =0,2\pi 3,4\pi 3$ angular positions are being singled-out. It comes clearly out from inspection of Fig. 10 that those three positions (as well as any other three related to this one by an arbitrary in-plane rotation) are equivalent, thus suggesting that the APP must exhibit a threefold symmetry pattern or rather a sixfold one when squaring the outgoing field into the corresponding considering intensity. Such a configuration has been used, for example, in Refs. 41 and 42.

The above projections along the $e\u20d7x$ and $e\u20d7y$ unit vectors become then irrelevant, and we turn now to the more appropriate projection onto the rotating unit vector $u\u20d7=e\u20d7x\u2061cos\u2061\alpha +e\u20d7y\u2061sin\u2061\alpha $,

This PSA configuration leads to the more “natural” sixfold [due to the squaring from the threefold $Pu\alpha $ to a then sixfold $Iu2\omega \alpha $ angular dependence].

## REFERENCES

_{2}monolayer

_{2}and h-BN by optical second-harmonic generation

_{2}Se

_{3}