Toy model of second harmonic generation due to structuring of centrosymmetric films

: We show how structuring of matter can lead to second order optical nonlinearity. Coulomb interactions involving bound electrons cause a nonlinear optical response at boundaries. We demonstrate that, in a planar structure cut from a centrosymmetric lattice of harmonic oscillators, second order nonlinearity is proportional to the perimeter of the structure. This proportionality and our model can instruct the design of dielectric nonlinear particles, surfaces and metamaterials for optical second harmonic generation.


Introduction
The second-order nonlinear interaction of two waves in a medium can generate a third wave at a combinational frequency, underpinning the phenomena of second harmonic generation, optical rectification, parametric down conversion, sum-and difference frequency generation. Conventionally, second harmonic is generated in crystals with symmetry lacking an inversion centre, however, centrosymmetric media can generate second harmonic waves in the presence of inhomogeneity associated with interfaces [1][2][3][4][5][6][7][8][9], optical field gradients [9][10][11] or chirality [12]. A considerable effort is now focused on artificial nonlinear materials that can be manufactured by top-down nanofabrication processes: second harmonic can be efficiently generated in structured metallic and dielectric films, as well as metasurfaces with asymmetric patterns even if the bulk of the constituent materials is centrosymmetric [13][14][15][16][17][18][19][20][21][22].
In this paper we report a classical oscillator model that gives rise to second-order nonlinearity in a structured dielectric film made of a centrosymmetric dielectric material. The film is structured into flakes (particles). We describe the dielectric material of the film as a twodimensional lattice of harmonic oscillators with optical electrons. By design, the oscillators do not exhibit any second order nonlinear response in isolation. However, we show that flakes of such structured films can generate second harmonic where nonlinearity emerges from the Coulomb interactions of charges of neighbouring optical electrons and nuclei in the confined anisotropic environment of the flake. The model leads to the following scaling rules: for a particle of a given shape and orientation (relative to the driving field), the first order polarizability is proportional to the particle surface while the second order polarizability is proportional to the particle perimeter.

Optical nonlinearity of structured dielectric films
In the model we consider a flake consisting of N atoms and describe atom k as a classical Lorentz oscillator [23] with an optical electron with coordinate rk(t) bound to a stationary nucleus at Rk. As is standard Lorentz model approximation, we assume that the electron and the nucleus are elastically bound which gives rise to a linear restoring force and the atom's linear optical response. In order to account for the influence of the N-1 other atoms on the optical electron k, we consider the Coulomb interactions of the optical electron with remaining electrons and nuclei. The resulting potential energy Uk is Where q and m are the electron charge and mass, and ω0 is the angular resonance frequency of the isolated harmonic oscillator. The resulting equation of motion in an optical field E(t) becomes Where γ is the damping frequency. A similar approach has been previously employed to study linear macroscopic properties of dielectrics [24] and for treating surface nonlinearity in uniform planar surfaces [7,[25][26][27][28]. Here, our focus is on highly structured materials. Second harmonic generation is strongly shape-dependent and our model is applicable to planar dielectric structures of any shape. Therefore, it complements existing models that describe second harmonic generation in specific materials and/or shapes [29][30][31][32], and symmetry-based selection rules [33,34] that only identify cases where second harmonic generation is forbidden. In order to predict the edge nonlinearity of a structured dielectric film, such as a metamaterial, 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 (see Fig. 1a). The electromagnetic response of the single atom is strictly linear and arises from a harmonic potential as described above. 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 a nonlinear response of optical electrons in collections of coupled atoms, which is illustrated by curved optical electron trajectories in Fig. 1b. We note, that similar nonlinearity arising from electrostatic interactions has been considered as 'multipole nonlinearity' in the context of nonlinear plasmonic metamaterials [35]. 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 hexagonal lattice. In order to increase the speed and stability of numerical modelling, a perturbative scheme is used. In the first stage, the linear response of N coupled atoms is found. This linear response is then used, in the second stage, as a source for determining the nonlinear response. In order to approximate a typical atomic response in dielectrics such as SiN, ITO and TiOx, we choose ω0≈9.4×10 15 s -1 and γ=0.01ω0 (This ω0 corresponds to a UV resonance at 200 nm wavelength). The amplitude of the driving electric field (E) is chosen to be 8.68×10 7 V/m, which corresponds to pump light of 1 GW/cm 2 intensity. The wavelength of the pump is that of the Nd:YAG laser (1064 nm, angular frequency ω≈0.19×ω0), chosen purely because it is a common wavelength for the nonlinear optics community and because it lies far outside the resonant region defined by the choice of ω0 and γ. The lattice constant, i.e. the smallest inter-atomic spacing, is 0.5 nm in all cases. We will calculate the total linear and non-linear electric dipole induced in collections of atoms -2D particles carved out of a hexagonal lattice by imposing a closed border. Given a displacement of the k th electron relative to its nucleus, rk -Rk, the electric dipole due to the k th atom is dk=q(rk -Rk), and the total electric dipole is p=∑k dk, which can be separated into linear and nonlinear components, p (1) , p (2) , …, oscillating at the driving frequency ω and its harmonics, 2ω, … .
Calculation times for the structures considered below -containing up to 3691 atoms -are up to 32 seconds in Matlab R2018b on a Windows 10 computer with Intel ® Core TM i5-6600 3.3 GHz CPU and 32 GB RAM.

Triangular particles
Irrespective of the origin of nonlinear response, the vectorial nature of the electromagnetic fields implies that a total second-order nonlinear response of a 2D (two-dimensional) particle driven by in-plane optical fields, i.e. a total induced non-linear electric dipole, is allowed only for structures with three-fold rotational symmetry, or no rotational symmetry at all [33,34]. We will therefore focus on two-dimensional particles with three-fold overall rotational symmetry, the simplest case that permits second-order nonlinearity, represented by equilateral triangular arrangements of atoms cut out of a hexagonal lattice.
We use the notation y xx d to denote per-atom second-order nonlinear electric dipole along the y-axis induced due to driving field polarized along the x-axis. The other components follow in the same way. Fig. 2 shows the second-order nonlinear electric dipole induced in each atom of a triangular 2D particle, cut out of a hexagonal lattice of atoms. Several important phenomena are readily observable: (1) the quadratic nonlinear dipole is always weakest in the central region of the triangle -as one would expect for a centrosymmetric arrangement of atoms (hexagonal lattice); (2) a nonlinear response at the edge seems to be induced in all cases with similar amplitude, but different phases of second harmonic response at different edges can yield vanishing (components of) total second harmonic response for a whole particle. x-polarized field (Fig. 2a) as well as y-polarized field (Fig. 2c) create a second-order nonlinear response in the x-direction along the edges of the triangle, but the nonlinear response of opposite edges cancels out, thus the x-component of the total nonlinear electric dipole is zero, as it should be for a particle with mirror symmetry, x ↔ -x. In contrast, the y-polarized second-order nonlinear response resulting from x-polarized (Fig. 2b) and y-polarized (Fig. 2d) illumination does not cancel.

Particle symmetry
Basic symmetry considerations [34] show that total second harmonic response is forbidden for 2-fold rotationally symmetric shapes such as a rectangle, but is allowed for 1-fold symmetric shapes such as isosceles triangle. An interesting question to ask is how this selection rule appears when a triangle is gradually converted via a trapezoid into a rectangle. How does the total nonlinear electric dipole vanish? How does it depend on changes in geometry (for a given driving field)? What effect does the finite lattice size have on the changes in the induced nonlinear electric dipole? To address this, a series of calculations has been carried out for symmetric triangular, trapezoidal and rectangular particles driven by the same electric field along the symmetry axis (y-axis). The results are shown in Fig. 3. We fix the height and bottom side length of an isosceles trapezoid, and increase the slope of its left and right sides by extending its top side from 0 nm to 30 nm, see Fig. 3a. The effect of such changes in particle symmetry, from three-fold rotational symmetry via absence of rotational symmetry to two-fold rotational symmetry, on the magnitude of the total linear dipole (p (1) ) and the total second-order dipole (p (2) ) is shown in Fig. 3b. The magnitude of the total linear dipole (p (1) ) grows linearly as the triangle evolves into a rectangle. This is consistent with the linear polarizability of each atom being roughly independent of its neighbourhood -so more atoms translate into correspondingly higher total linear electric dipole. Indeed, we find that the linear dipole per atom remains constant. Clearly, such a simple response is a consequence of operating in the off-resonance regime (the driving frequency to atomic resonant frequency ratio is ω/ω0≈0. 19). In contrast, the magnitude of the total quadratic electric dipole (p (2) ) decreases linearly with the increase in the top side length. This can be explained by considering the plots of the per-atom  (Fig. 3c,d). In the rectangle, the second-order electric dipole excitation at the top and bottom sides is of equal magnitude but opposite phase, giving zero total effect. In case of the triangle, this cancellation does not occur. The gradual drop in the magnitude of the total second-order dipole with increasing top side length, shown in Fig. 3b, is consistent with the cancellation of second-order dipole contributions from atoms in corresponding sections of the top and bottom sides, which leads to complete cancellation of the particle's second-order nonlinear response as it becomes rectangular. Thus, while local edge nonlinearity is also present for particle shapes with inversion symmetry, anti-phase nonlinear response of opposite edges causes cancellation of the overall second-order nonlinear response of inversion-symmetric particles.  (1) and p (2) at the fundamental and second harmonic frequencies induced by a y-polarized fundamental wave in a particle evolving from lack of inversion symmetry to inversion symmetry. Colour maps (c) and (d) show the y-component of the atomic dipole y yy d at frequency 2ω for the initial and final particle shapes.
We note that second-order nonlinear response with same magnitude and opposite phase for opposite edges (one rotated 180° relative to the other) has an interesting implication for inverse structures (a particle and the corresponding particle-shaped hole). An inverse structure results from interchanging which side of the edge the material is on, which is equivalent to a 180° rotation of each segment of particle edge. Therefore, inverse structures must have second-order dipoles of same magnitude and opposite sign.

Particle size
The dipoles in Fig. 3 show remarkably smooth dependencies on particle geometry, despite the fact that the length of the longest side of the rectangle is just 61 atoms. One can therefore conjecture that whilst the nonlinear response at each point on the edge may be strongly influenced by the local orientation of the edge-cut relative to the driving field and atomic lattice, the total second-order nonlinear response can be approximated by treating edges as 'smooth' even in the presence of a coarse lattice. This conjecture can be tested by choosing a particle of some shape, changing its overall size (whilst keeping the lattice size fixed), driving the particle with electric field of fixed magnitude and orientation, and calculating the total induced electric dipole. The results of this test, for the case of triangular particles, are shown in Fig. 4. The linear and second-order nonlinear induced electric dipoles have been calculated for triangles with side lengths from 1 nm (6 atoms in total) to 17 nm (630 atoms). Despite the coarse lattice (lattice constant of 0.5 nm), the total linear electric dipole is proportional to the number of atoms, i.e. the linear electric dipole per atom does not change much, see inset. In contrast, the total secondorder nonlinear electric dipole is proportional to the size S of the triangle's edges, i.e. to the particle's perimeter. Thus, the overall quadratic nonlinear response in mesoscopic twodimensional particles, with centrosymmetric structure of the bulk, is proportional to the perimeter of the particle. We note that size effects of (surface) second-order nonlinear response of metallic particles have been investigated experimentally in the past, where an increase in second-order nonlinear response with particle size has been observed for very small particles [36]. The noteworthy effect detected here is that such behaviour persists to levels of just few atoms in dielectric particles.

Optimization strategy
The proportionality of the total induced second-order nonlinear electric dipole to the perimeter of the particle suggests a simple strategy of optimizing the non-linear response of mesoscopic structures. As illustrated by Fig. 5, given a 2D particle with nonlinear response, such as an equilateral triangle, an increased second-order nonlinear response can be achieved by substituting the particle (triangle) with smaller particles of the same shape (triangles), which have a larger combined perimeter and the same combined area. This process will increase the total length of edges, thus increasing the second-order nonlinear response, but will leave the overall number of atoms (approximately) unchanged, thus preserving the linear response. We note that this method for maximizing the nanostructure's nonlinear response will naturally lead to a metasurface, i.e. a periodically structured film (Fig. 5). While a particle that is small in comparison to the wavelength will radiate over a large solid angle, second harmonic generation from a 2D array of nanostructures will become increasingly directed with increasing overall size of the array. The direction(s) of second harmonic generation will be determined by phase matching, i.e. by constructive interference of the fields radiated by all unit cells of the array.
While we consider 2D structures, we expect that our method can be extended to cover structured films, that are thin in comparison to the wavelength under consideration, by treating interactions between columns of atoms rather than individual atoms. For such structures, the overall nonlinear response will increase with increasing thickness due to the increasing number of atoms on the structure's perimeter. However, as the thickness increases further, phase matching along the thickness direction would need to be included to describe the overall nonlinear response. The magnitude of total linear dipole and total second harmonic dipole for the cases considered in (b). The total linear dipole is almost the same in all cases, since the number of atoms in the three arrangements is almost the same. The total SH dipole is proportional to the total edge length, i.e. it doubles with every step.

Conclusion
In summary, we have shown that second-order nonlinear optical response can arise from the nonlinearity of inter-atomic Coulomb interactions -a universal mechanism, applicable to all media. Coulomb interactions between neighbouring atoms cause nonlinear oscillations of charges at the edges of 2D particles in response to light. Patterning enables second-order nonlinear response even in case of normal incidence onto centrosymmetric dielectric films, where the second-order nonlinear response is proportional to the perimeter of the patterned shape. The proportionality is extremely robust and persists down to the level of just few atoms. Edge nonlinearity, as described here, should be expected in a wide variety of metamaterials, photonic crystals and other planar nonlinear nanophotonic structures driven by light at normal incidence. It is of direct relevance to the rapidly growing range of applications of twodimensional dielectric metamaterials (metasurfaces) in nonlinear photonics and quantum optics.

Funding
This work is supported by the UK's Engineering and Physical Sciences Research Council (Grant No. EP/M009122/1).