Toy model of harmonic and sum frequency generation in 2D dielectric nanostructures

Optical nonlinearities of matter are often associated with the response of individual atoms. Here, using a toy oscillator model, we show that in the confined geometry of a two-dimensional dielectric nanoparticle a collective nonlinear response of the atomic array can arise from the Coulomb interactions of the bound optical electrons, even if the individual atoms exhibit no nonlinearity. We determine the multipole contributions to the nonlinear response of nanoparticles and demonstrate that the odd order and even order nonlinear electric dipole moments scale with the area and perimeter of the nanoparticle, respectively.

In nonlinear media, electromagnetic waves can interact, generating new waves at combinatorial frequencies.
Depending on the generated frequency, the nonlinear phenomena can be classified as harmonic generation, parametric down conversion, sum or difference frequency generation, which are important for optical signal processing, bio-imaging and spectroscopy. Even order harmonics are usually forbidden in materials with inversion symmetry, but with the presence of inhomogeneity associated with interfaces [1][2][3][4][5][6][7][8][9] , optical field gradients 10,11 or chirality 12 , centrosymmetric media can support even harmonic response. In particular, nanostructuring can introduce inhomogeneity and asymmetry in metallic and dielectric films, resulting in even order sum-frequency and harmonic generation in structures made from centrosymmetric media [13][14][15][16][17][18][19][20][21][22] . Such nonlinearities depend strongly on geometry, i.e. the shape of particles or the unit cell of periodic structures. The advance of nanofabrication to smaller and smaller scales provides an opportunity to engineer the nonlinear properties of particles and (meta)surfaces on the nanoscale. Such engineering depends on understanding the dependence of nonlinear optical properties on the arrangement of atoms. Full ab-initio quantum treatments of second- [23][24][25][26][27][28][29] and higherorder 30,31 nonlinearity have been demonstrated, but only for specific materials and under specific conditions, whilst experimental measurements, of surface nonlinearity, require pristine surfaces and strict control over the bulk effects 3,5,[32][33][34] . Sophisticated classical/quantum test models have been previously developed for metals 35,36 , 'dipolium' dielectrics 29,37 , and dielectric particles of quite specific shapes 38 . Meanwhile, symmetry-based selection rules 39,40 identify cases where even harmonic generation is forbidden, but cannot predict its magnitude where it is allowed.
Therefore, the motivation of this work is to develop a simple classical model that can predict how harmonic and sum frequency generation depend on the size and geometry of dielectric nanostructures. We consider a 2D nanoparticle described as a lattice of linear oscillators coupled by Coulomb interactions. The model reveals the different origins of the even and odd order nonlinear processes and predicts the dependence of the efficiency of nonlinear processes on the incident polarization state and the nanoparticle's geometry, symmetry and size. We show that the even order and odd order nonlinear electric dipole moments scale linearly with the perimeter and area of the nanoparticle, respectively, and quantify magnetic dipole and electric quadrupole contributions to nonlinearity. The model may be applied to optimize and understand harmonic and sum frequency generation in dielectric nanostructures of various shapes. www.nature.com/scientificreports/ is described by a linear restoring force, which gives rise to a linear optical response. In this model, there is no nonlinearity within a single atom (see "Discussion"). However, if we consider a particle, which contains N atoms, the Coulomb interaction of electrons with electrons, and electrons with nuclei of other atoms will introduce nonlinearity. The equation of motion for each optical electron can be written as where r k and R k is the positions of electron k and its nucleus in the particle; γ is the damping frequency; q and m are the electron charge and mass; and ω 0 is the angular resonance frequency of an isolated atomic harmonic oscillator. Second harmonic generation (SHG) due to Coulomb interactions has been observed experimentally in dimers of metallic nanoparticles of deeply sub-wavelength size 43 and a similar nonlinearity arising from electrostatic interactions has been considered in the context of nonlinear plasmonic metamaterials 44 . Similar approaches have been previously employed to study linear 45 and nonlinear properties of dielectrics 46 , harmonic generation by metallic structures 47 and anharmonic oscillators 48,49 , and for treating surface nonlinearity of uniform planar surfaces 7,[27][28][29]37 . A range of models of sum frequency generation has been reported, e.g. for harmonic oscillators 50 and liquid interfaces [51][52][53] . Nonlinear oscillator models predicting the nonlinear response of metacrystals have been reported 54,55 and applied to split ring resonator magnetic media. Here, our focus is on highly structured materials and our model is applicable to planar dielectric nanostructures of any shape. Therefore, it complements existing models that describe harmonic or sum frequency generation in specific materials and/or shapes 16,23,35,36,38 , and symmetry-based selection rules 39,40,56 that predict cases where such processes are forbidden.
In order to predict the nonlinear response of a structured dielectric film, such as a nanoparticle or metamaterial unit cell, we consider plane wave illumination (along z) of a 2D (two-dimensional) lattice of 'atoms' consisting of charges that are confined to the xy-plane. The electromagnetic response of the single atom is strictly linear and arises from a harmonic potential 41 . The nonlinear response arises exclusively due to inter-atomic Coulomb interactions. The ∝ 1/r 2 dependence of the Coulomb force, where r is the inter-charge separation, gives rise to a nonlinear response of optical electrons in collections of coupled atoms, which is illustrated by curved optical electron trajectories in Fig. 1.
We apply this basic principle of one-to-one coupling via the Coulomb force to modelling of the response of large collections of atoms that form particles with a square or hexagonal lattice. The resonance frequency of an isolated atom is ω 0 = 9.4 × 10 15 rad/s, corresponding to a UV resonance at 200 nm wavelength, and the damping frequency is γ = 0.01 ω 0 to approximate a typical atomic response in dielectrics such as ITO, TiO x and SiN. In order to access a highly nonlinear regime, the amplitude of the pumping electric field (E 0 ) is set to 8.68 × 10 10 V/m, corresponding to 10 15 W/cm 2 pump intensity. While this extreme level of intensity may not be feasible experimentally, it allows us to explore physical trends and higher-order optical nonlinearities in a regime where they www.nature.com/scientificreports/ are of sufficient magnitude to make numerical errors of differential equation solvers irrelevant. The pump wavelength is chosen to be 1064 nm, a common choice for nonlinear optics that is far from the atomic resonance (angular pump frequency ω p = 0.19 ω 0 ). In modelling of sum frequency generation, we pump the particle at two frequencies, ω p and 0.1ω p , simultaneously. The ratio of 10 of the two pump frequencies is chosen to make the sum frequency order easily recognizable. The atomic spacing is chosen to be 0.5 nm in all cases. By solving the coupled equations of motion [Eq. (1)] of all optical electrons within 2D nanoparticles numerically, we calculate the displacement of each optical electron in the time domain. We solve the coupled system of differential equations in Matlab starting without optical electron displacement and allowing transitional effects to pass before analysing the oscillation of all electrons over 1000 periods of the driving field. This allows the calculation of the particle's net multipole moments. In the interest of simplicity, our initial discussion will be limited to the particle's electric dipole moment in order to describe radiation emitted normal to the plane of the particle. The electric dipole moment is proportional to the displacement of electrons. Therefore, a displacement of the k th electron relative to its nucleus, r k − R k , results in an electric dipole due to the k th atom of d k = q(r k − R k ), and the total electric dipole of the particle is the sum over all atoms, P = ∑ k d k . This gives the electric dipole moment d(t) of each atom and the net electric dipole moment P(t) of the nanostructure in the time domain. The net electric dipole moment is separated into linear and nonlinear components, P (1) , P (2) , P (3) …, oscillating at the driving frequency ω p and its harmonics 2ω p , 3ω p …, by Fourier transformation. In the presented Fourier series, the peak values are the dipole amplitudes in Coulomb meters at harmonic and sum frequencies. For pumping of a 2D particle at normal incidence, as considered here, radiation emitted normal to the plane of the particle is determined by its net electric dipole moment alone. Towards the end of this paper, we will consider other multipoles, which contribute to radiation in other directions.
Harmonic and sum-frequency generation in nanostructures of D4 symmetry. Here we consider how a structure in the xy-plane radiates along z in response to pumping along z. According to symmetry analysis, radiation of even order harmonics along z is forbidden in structures with even-fold rotational symmetry and allowed in structures of no or threefold rotational symmetry 39,40 . In order to reveal how nonlinear frequency generation appears/vanishes in cases where it is allowed/forbidden, we study harmonic and sum frequency generation in simple 2D particles of either square shape (D4 symmetry, with square lattice) or triangular shape (D3 symmetry, with hexagonal lattice). Figure 2a illustrates harmonic generation in a square particle of 25 atoms. The frequency-dependence of the total electric dipole moment P generated by electric pump field E 0 cos(ω p t) is shown, where P i,j refers to the i-component of the dipole moment caused by a j-polarized pump field, with i,j being x or y. As required by symmetry, x-polarized excitation generates a dipole only along x, while y-polarized excitation generates a dipole only along y with the same magnitude. i.e. P x,x = P y,y are observed, while P x,y and P y,x are forbidden. As expected for a particle with even-fold rotational symmetry 39,40 , we observe only odd harmonics-P (1) , P (3) , P (5) , P (7) and P (9) -which appear as peaks in Fig. 2a.
Harmonic generation is a special case of sum frequency generation. The latter becomes apparent when pumping the square particle with a combination of two frequencies, ω p and 0.1ω p . This is illustrated by Fig. 2b, which considers a pump field E 0 /2[cos(ω p t) + cos(0.1ω p t)]. As for harmonic generation, also mixed frequency pumping results in an electric dipole moment of the particle that is parallel to the pump polarization, P x,x and P y,y . We observe only odd order sum frequencies, i.e. odd harmonics of the higher and lower pump frequencies (ω p , 3ω p , 5ω p , … and 0.1ω p , 0.3ω p , …) and sums of odd combinations of both pump frequencies (1.2ω p , 2.1ω p , 2.3ω p …), indicating that the model can predict wave mixing in 2D particles. We observe that different dipole moments of the same order tend to be of similar magnitude, but resonant enhancement near the atomic resonance at 5.26 ω p and the background level also influence the magnitude of the calculated dipole moments. (Similar results are obtained for different choices of pump frequency, see Supplementary Fig. S1).
The origin of harmonic and sum frequency generation is revealed by the distribution of atomic dipole moments at the relevant frequencies. Figure 2c illustrates this for the square particle pumped at frequency ω p . The top row shows the distributions of dipole components oriented normal to the pump polarization at odd (left) and even (right) harmonic frequencies. Notably, all harmonics are generated locally, but the symmetry of the particle implies that the atomic electric dipole component normal to the pump field has the same magnitude and a π phase difference on opposite sides of the structure, resulting in overall cancellation. Therefore, the structure of D4 symmetry has no net electric dipole moment in the direction normal to the pump polarization. The net electric dipole component parallel to the pump polarization also cancels for even harmonics (bottom row, right). In contrast, the dipole component parallel to the pump polarization does not vanish for odd harmonics (bottom row, left). In-phase oscillation of all atoms causes the atomic dipoles to add up to a net electric dipole moment of the particle at odd harmonic frequencies.
Indeed, the symmetry of the experimental arrangement, including structure and pump beam, means that no light can be emitted at even harmonics along z and therefore the net electric dipole moment must vanish. In these cases, the fields "radiated" from different parts of the particle cancel by destructive interference in the normal direction. Radiation in other directions may not vanish, can be described in terms of other multipoles, and will be discussed towards the end of this paper. Here, we describe radiation emitted normal to the plane of the particle.
While the allowed component of odd harmonic generation originates from the entire particle, we notice that the even harmonics-which do not add up to any net electric dipole moment-are associated with edges of the particle. In the interest of simplicity, we have shown dipole distributions for harmonic frequency generation in Fig. 2c. Dipole distributions for sum frequency generation with different pump frequencies (see Supplementary  Fig. S2) have the same characteristics as those for harmonic generation of the same order. www.nature.com/scientificreports/ The dependence of the nanoparticle's odd harmonic response on its size is shown by Fig. 2d. We find that the total nonlinear electric dipole at the generated odd harmonics is proportional to the total number of atoms (2D particle area). The pump power dependence of the generated harmonics is illustrated by Fig. 2e, which shows the electric dipole moment generated in the square particle for different pump field amplitudes at different harmonic frequencies. The magnitude of the electric dipole moments at the first, third and fifth harmonic frequencies-P (1) , P (3) , P (5) -is proportional to the first, third and fifth power of the pump amplitude, as expected for first, third and fifth order (non)linear optical processes.
Harmonic and sum-frequency generation in nanostructures of D3 symmetry. While the evenfold rotational symmetry discussed above only permits odd order nonlinear electric dipole moments, three-fold rotational symmetry allows a more complex nonlinear response. Now we will discuss harmonic generation in a particle of D3 symmetry, a triangular particle consisting of 28 atoms arranged in a hexagonal lattice as shown by Fig. 3. The spectral dependence of the triangular particle's electric dipole moment in response to pumping at frequency ω p is shown by Fig. 3a-d. Also here, the electric dipole moment at the kth harmonic is proportional to the kth power of the pump field amplitude, causing higher harmonics to vanish quickly with reducing pump amplitude, as it should be (see Supplementary Fig. S3). An x-polarized pump field generates an electric dipole moment of the particle along x at odd harmonic frequencies (Fig. 3a) and along y at even harmonic frequencies (Fig. 3b). A y-polarized pump field generates a dipole moment along y only at all harmonic frequencies (Fig. 3c,d). The distribution of the atomic dipole moments within the nanoparticle is shown for the second and third harmonic frequencies (insets), which are representative of the even and odd harmonic cases. While the particle exhibits local charge oscillations along x and y at all harmonics, overall cancellation of components of  . (a, b) Frequency dependence of the electric dipole moment of a square particle (blue outline) consisting of 25 atoms (red crosses) in response to optical pumping at (a) a single frequency ω p and (b) a combination of two frequencies, ω p and 0.1ω p . P y,x indicates the y-component of the particle's electric dipole moment caused by x-polarized pumping. For either x-or y-polarized pumping, the orientations of pump polarization and generated dipole are parallel, and the generated dipoles have the same magnitude for both cases. (c) Magnitude (colours) and sign ("+" and "−") of the dipole moment per atom, d generated at harmonic frequencies (multiples of ω p ) in response to pumping at frequency ω p . The top (bottom) row shows the dipole component orthogonal (parallel) to the pump polarization. Stacked images for different-either even or odd-harmonics show the same qualitative behaviour. (d) Scaling of the particle's electric dipole moment with its size at the first three odd harmonic frequencies [corresponding to peaks in (a)]. These odd order dipole moments scale with the particle's number of atoms (area). (e) Scaling of the particle's electric dipole moment with the pump field amplitude for the same cases. www.nature.com/scientificreports/ the particle's electric dipole moment at even (Fig. 3a), odd (Fig. 3b) or all (Fig. 3d) harmonics results from charge oscillations with equal amplitude and opposite phase at opposite edges of the triangle. It has been reported that the total linear and second-order nonlinear responses of a particle scale differently with its size. For a 3D particle made from centrosymmetric material, the total linear electric dipole moment is proportional to its volume (number of atoms), while the second-order dipole is proportional to its surface area 32 . For a 2D particle, as considered here, the linear electric dipole moment is proportional to its area (number of atoms), while the second-order electric dipole is proportional to its perimeter 41 . However, how does particle size affect higher order harmonics? Fig. 3e-h shows how the electric dipole moment of a triangular 2D particle changes with particle size. The magnitudes of the total dipole moments, P (1) , P (2) , P (3) , P (4) , at the first four harmonic frequencies are shown for triangular particles with perimeter from 1.5 nm (3 atoms in total) to 9 nm (28 atoms). We find that the total nonlinear electric dipole at odd harmonics is proportional to the total number of atoms (2D particle area), while the dipole at even harmonics is proportional to the particle's perimeter. At the first three harmonic frequencies, x-polarized and y-polarized pump light induce dipole moments of the same magnitude, while at the fourth harmonic P (4) y,y = 4P (4) y,x . Sum frequency generation in the triangular particle is consistent with harmonic generation. An x-polarized pump electric field at a combination of two different frequencies generates a nonlinear electric dipole moment along x for odd sum frequencies and along y for even sum frequencies, while a y-polarized pump field generates a nonlinear dipole moment along y only for all sum frequencies, where more and more sum frequencies are generated with increasing pump amplitude (Supplementary Fig. S4). For sum frequency generation resulting from cross-polarized pump fields of different frequencies see Supplementary Fig. S5. In this case, x-polarized electric dipole moments are generated at sum frequencies containing an odd multiple (1, 3, …) of the x-polarized pump frequency and generated y-polarized sum frequencies contain an even multiple (0, 2, …) of the x-polarized pump frequency. The inset shows the outline (blue line) and the 28 atoms (red crosses) of the triangular particle. Colourmaps show the distribution of the atomic dipole moments at the 2nd and 3rd harmonic frequencies, which are representative of even and odd harmonics. Colour represents the dipole magnitude and regions oscillating in anti-phase are indicated by opposite signs as in Fig. 2c. (e-h) Scaling of the particle's dipole moment with its size at the first four harmonic frequencies. The particle's electric dipole moment scales with its number of atoms (area) at odd harmonics and with its perimeter length at even harmonics. www.nature.com/scientificreports/ Multipole contributions to harmonic generation. While radiation normal to the plane of the 2D particle pumped at normal incidence is determined by its electric dipole moment, other multipoles 57 contribute to radiation in other directions. The α-component of the particle's net electric dipole moment is given by where ρ(r) is the charge density as a function of position r . In our 2D case, only P x and P y are allowed and the radiated power is given by the following expression, where Z 0 is the impedance of free space, c is the speed of light and k = ω/c The α-component of the particle's magnetic dipole moment is where J is the current density and ṙ i is the electron velocity. Here, only m z is allowed and the radiated power is The α, β-component of the particle's electric quadrupole is given by where δ α,β is the Kronecker delta. In our geometry, Q x,x , Q y,y and Q x,y = Q y,x are allowed and the radiated power is Figure 4 shows the power radiated by different multipoles for the square and triangular particles, where the coordinate origin is placed at the centre of each particle. For the square particle, electric dipole radiation is emitted at odd harmonics and considerably weaker electric quadrupole radiation ( Q x,x , Q y,y ) is emitted at even harmonics, while magnetic dipole contributions are negligible, Fig. 4a. Therefore, only odd harmonics will be radiated along z. In other directions, all harmonics will be radiated, with even harmonics being several orders of magnitude weaker than odd harmonics of similar order. For the triangular particle, electric dipole radiation dominates at odd harmonics, while electric dipole and quadrupole radiation have a comparable magnitude at the second harmonic frequency, Fig. 4b,c.
The multipole contributions are also revealed by the polarization surface densities (electric dipole distributions) in Figs. 2c and 3a-d. For example, for both the square and triangular particles, the net quadrupole at even harmonics originates from electrons on opposite edges of the particle, that oscillate towards and away from each other (→ ←, ← →). The net magnetic dipole at the third harmonic frequency for x-polarized pumping of the triangular particle arises from anti-parallel electron oscillations (↑↓,↓↑), see Fig. 3b.

Discussion
The presented model describes the contribution of Coulomb interactions between different atoms to the nonlinear response of dielectric nanostructures. To focus on the universal collective nonlinear response of any dielectric, we purposefully remove all possible atom-specific nonlinearities by describing atoms with the harmonic potential of a classical Lorentz oscillator model (and pumping far away from resonance). However, real atoms can be nonlinear. Considering the inversion symmetry of atoms in the ground-state, opposite directions of electron displacement yield opposite restoring forces. Therefore, the simplest approximation of a nonlinear atom would involve a symmetric anharmonic potential. Within this approximation, an individual atom can only have odd orders of optical response, i.e. contribute to 3rd, 5th, … order nonlinearity. With respect to the particle's electric dipole, such atomic nonlinearity would contribute to the odd-order components that are parallel to the excitation field. This atomic contribution to a particle's odd order nonlinear dipole moment will be proportional to the particle's number of atoms, just like the contribution from Coulomb interactions between atoms (Figs. 2d,  3g). Therefore, the atomic nonlinearity will affect the magnitude of the particle's odd order nonlinear electric dipole response parallel to the pump field, but not its qualitative behaviour. It has no bearing on the particle's even order nonlinearity, which arises entirely from interactions between atoms.

Conclusion
In summary, we have shown that harmonic and sum frequency generation in 2D dielectric nanostructures can be modelled in terms of harmonic oscillators (atoms) coupled by the Coulomb force. This approach allows modelling of large collections of atoms and mapping of the origin of harmonic and sum frequency generation in nanostructures. Odd order sum frequency generation has been found to originate from the whole structure with corresponding net electric dipole moments scaling with the number of atoms. Even order sum frequency generation has been traced to the boundaries of 2D nanostructures, with corresponding net electric dipole moments scaling with the structure's perimeter. The model allows the design and optimization of nonlinear dielectric nanostructures for nanophotonics, such as structures and particles cut from 2D materials and dielectric layers that are thin compared to the pump wavelength. 3D geometries may also be considered by solving Eq. (1) in three dimensions.

Data availability
Data from this paper is available from the University of Southampton ePrints research repository: https:// doi. org/ 10. 5258/ SOTON/ D1792.