Quo Vadis Nonlinear Optics? An Alternative and Simple Approach to Third Rank Tensors in Semiconductors

Abstract

It is well understood that nonlinear optical (NLO) phenomena are deeply related to the material’s symmetry. Mathematically, the material symmetry can be described in terms of the nonzero parameters in the nonlinear susceptibility tensors. Generally, more complex structures involve more nonzero parameters in the tensor. The number of parameters increases rapidly if higher NLO orders are considered, complicating the physical analysis. Conventionally, these parameters are obtained via abstract symmetry analysis, e.g., group theory (GT). This work presents a novel theoretical analysis to approach the nonlinear tensor using the simplified bond hyperpolarizability model (SBHM) and compare it with GT. Our analysis is based on a light–matter interaction classical phenomenological physical framework. Rather than just looking at the symmetry of the crystal, the model applies physical considerations requiring fewer independent parameters in the tensor than GT. Such a simplification significantly improves the determination of the surface–bulk SHG contribution factors, which cannot be extracted from the experiment alone. We also show for the case of perovskite that the SHG contribution can be addressed solely from their surface dipoles with only one independent component in the tensor. Therefore, this work may open the path for a similar analysis in other complicated semiconductor surfaces and structures in the future, with potential applications to nanoscale surface characterization and real-time surface deposition monitoring.

1. Introduction

After the discovery of second-harmonic generation (SHG) by Franken et al. [1], progress in the field of nonlinear optics (NLO) increased at a rapid pace with valuable applications ranging from thin-film material characterization, nanoscale optical sensing, nanoscale engineering, and various other fields. The nonlinear effects originate from the material response to the incoming field. This response is proportional to the electrical field and its higher powers. In the case of SHG, a single incoming wavelength involving two-photon absorption results in a nonlinear response that produces light with a frequency twice the incoming light frequency. Therefore, the amplitude of the SHG frequency generated by a sample is proportional to the square of the incoming field. It soon became apparent from the early theoretical works by Armstrong [2] and Bloembergen [3] that the analysis of the higher-order susceptibility tensorial components is essential to understand the nature of SHG generation in inversion-symmetric material.

It is now well understood that centrosymmetric structures such as silicon (Si) materials possess the vital property of bulk inversion symmetry, where parity symmetry forbids the nonlinear dipole contribution inside the bulk [3,4]. However, the contribution from SHG at the surface can exist due to symmetry breaking [5,6]. This discovery, together with the fact that NLO enables probing surfaces in situ and nondestructively, makes it a powerful tool to study semiconductor surfaces in detail and at the nanoscale. One of the major theoretical challenges in understanding SHG from semiconductor surfaces by optical probes, e.g., by SHG rotational anisotropy (RA), lies in separating the various nonlinear contributions that may arise. For the case of noncentrosymmetric semiconductor material, the SHG bulk dipolar contribution usually dominates due to the sheer amount of bulk layers inside the material. However, when the semiconductor possesses inversion symmetry, the sources become more complex and detailed knowledge of the material symmetry is required. Because physics takes place in real time and space as well as at the nanoscale (e.g., atomic level), a description of the phenomena using quantum mechanics seems natural and more fundamental, but requires a substantial amount of complex computations and lacks a pictorial insight [7]. Therefore, phenomenological models are sometimes combined with quantum mechanical calculations. Examples include the work of Mendoza and Mochán [8], and later Arzata and Mendoza [9], which exhibit interacting dipole bond model to investigate SHG from surface local-field effects. Other modeling efforts involving SHG from centrosymmetric material are mainly phenomenological, applying Green’s function formalism developed by Sipe et al. [10,11,12].

The challenge, however, was the number of independent parameters required to model SHG in inversion-symmetric structure that was still complicated, even if the SH signal was only generated from the surface/interface. This resulted in the difficulty to interpret experimental data in a simple way. This issue was resolved by the development of the simplified bond hyperpolarizability model (SBHM) proposed in 2002 by Dave Aspnes through a breakthrough paper [13], where direct products of bond vectors replace the tensorial formalism.

In this work, we show how the SBHM can produce a significant improvement over previous models using tensorial analysis. First, the SBHM approaches the nonlinear tensor components from a classical physical model rather than an abstract group theoretical calculation, which often hides the physics, thus providing a better understanding of the nonlinear phenomena at the atomic scale. Second, the SBHM relates some previously assumed to be independent tensor components obtained by GT because of the different approaches in obtaining them (through bond vector products), thus substantially reducing the number of independent parameters. For example, a remarkable achievement of the SBHM was to reproduce rotational anisotropy second-harmonic generation (RASHG) intensity experimental data from [14] with a high precision by only three normalized complex hyperpolarizabilities and two relative phases. The result was a reduction from the previous 14 normalized tensor coefficients in [10]. Third, because of this simplification, previous RASHG experiments can be better understood. Furthermore, the SBHM might resolve several important issues regarding the long-standing problem of determining the actual contribution of surface and bulk effects, which was previously difficult due to a large amount of fitting parameters [10] even in highly symmetric semiconductors such as Si and gallium arsenide (GaAs). The SBHM is thus opening the way for better surface characterization and real-time surface deposition monitoring at the nanoscale. Another attractive future potential application of this work is to study tensorial phase control of nonlinear metasurfaces, which enables SHG analysis at the subwavelength level with potential application in ultrathin, free-space photonic devices for nonlinear imaging [15,16].

To better appreciate the novelty of the third point, we mention the work given by Peng et al. [17], who applied four steps of the Ewald–Oseen dipole model to SHG. His work shows additional SH contributions in the form of spatial dispersion and magnetic terms. Applying the SBHM, we showed very recently that spatial dispersion contribution in Si(111), due to the incoming electric field decay, also contributes to the RASHG intensity data in nonvicinal Si(111) surfaces [18]. This effect can be described using a third rank tensor by redefining spatial dispersion from an earlier work by Peng et al. [17], which described spatial dispersion in the form of a fourth rank tensor. It is well known that the nonlinear polarization of the SH field depends on the material symmetry, which is described either by a third- or fourth-order nonlinear susceptibility tensor. It is interesting to study how this tensor is obtained from the SBHM and compared with the conventional group theoretical approach. Therefore, in this work, we access this tensor from both viewpoints and compare them later on. More specifically, we show that the SHG Si(111) surface part obtained from the SBHM only consists of two independent fitting parameters rather than the four independent parameters if GT is applied. Thus, the SBHM better estimates the surface to bulk SHG contribution.

However, despite the various advantages of this model previously mentioned, we also have to point out several shortcomings: (1) the SBHM is a phenomenological model that requires experimental fitting to obtain the nonlinear hyperpolarizability. Ab initio models using quantum mechanics should be applied to solve this issue but lie beyond the scope of this work. (2) This work only assumes oscillation parallel to the bond and neglects perpendicular oscillation. The latter might be significant for third-harmonic generation studies, but the exclusion of such an assumption for SHG still shows excellent agreement with experiments. (3) For complex semiconductors such as perovskite, the determination of the bond vector can be challenging, and calculations can be consuming. However, we show that many bond vectors cancel out when symmetry is considered, and simplified bond vectors can be applied instead.

An experimental comparison with the SBHM that included quadrupolar contribution from the Si bulk was performed using arbitrary input polarization to reproduce RASHG intensity in nonvicinal Si surfaces [19]. It was later shown that these effects were related to the first-order derivatives of the microscopic response function [20]. The contribution of spatial dispersion to SHG due to the incoming field was also incorporated very recently [18]. Investigation using this model to study the effect of H² exposure on the RASHG spectra for Si surfaces was also performed using the SBHM with success [21]. An SBHM study on the SH spectra from zincblende sample was also performed where the dominant SH contribution was indeed coming from bulk dipole radiation [22].

For the sake of clarity, we proceed with explaining the basics of the SBHM, followed by a description of group theory. Afterwards, we explain the main results of this work, namely, comparing the third rank tensor using both approaches for the case of a Si(111) and Si(001) surface and showing how this method can also be applied to more complex structures, e.g., perovskite. We then explain how this result can be applied to determine the different surface and bulk contributions in Si and show how it can increase our understanding regarding the polar–nonpolar nature of perovskite. Finally, a brief conclusion is given at the end.

2. Simplified Bond Hyperpolarizability Model

The simplified bond hyperpolarizability model (SBHM) is a classical phenomenological theory, whose history can be traced back to Clausius–Mossotti polarizable point models. This equation relates the atomic-scale polarizability $\alpha$ and the volume density n of a collection of points, which in SBHM are renamed as bonds, to the macroscopic complex dielectric function [7]. In the SBHM, the dielectric response function is derived from the Newtonian equation of motion of a charge $q_j$. The charges oscillate harmonically and anharmonically parallel to the covalent atomic bonds. The oscillation is driven by the local incoming electric field as illustrated in Figure 1. For the illustration, we chose a Si(001) surface, meaning that the Si surface of the film is a (001) Miller plane through the Si lattice (e.g., in the cubic lattice system, the direction (hkl) defines a vector direction normal to the surface of a particular plane or facet). For simplicity, let us consider a 1D Lorentzian equation of motion for the case of second-harmonic (SH) and third-harmonic (TH) generation [23]:

$$\begin{array}{c}\hfill m\frac{d^2\mathbf{x}}{d{t}^2}=q_j{\mathbf{E}}_{\mathbf{loc}}\left(\omega \right)·{\widehat{\mathbf{b}}}_j{e}^{-i\omega t}-{\kappa }_1(\mathbf{x}-{\mathbf{x}}_0)-{\kappa }_2{(\mathbf{x}-{\mathbf{x}}_0)}^2-{\kappa }_3{(\mathbf{x}-{\mathbf{x}}_0)}^3-\gamma \frac{d\mathbf{x}}{dt}\end{array}$$

(1)

here $\mathbf{x}$ is the charge position in the local coordinate system along the ${\widehat{\mathbf{b}}}_j$ direction, ${\mathbf{x}}_0$ is the charge equilibrium position, ${\mathbf{E}}_{loc}\left(\omega \right)$ is the local field experienced by the charge at the position of the atom, which is related to the macroscopic field through ${\mathbf{E}}_{loc}\left(\omega \right)=\frac{ϵ\left(\omega \right)+2}{3}\mathbf{E}\left(\omega \right)$. In Equation (1), ${\kappa }_1$, ${\kappa }_2$, and ${\kappa }_3$ are, respectively, the harmonic, first anharmonic, and second anharmonic spring constant. These constants will give rise to the harmonic, SH, and TH radiation. A frictional coefficient $\gamma$ is introduced in the equation to account for the dissipation. Here the direction of the j-th bond is defined by the unit vector ${\widehat{\mathbf{b}}}_j$.

Before we proceed with finding the solution for ${\mathbf{x}}_i$, it is necessary to discuss the specific assumptions in this model. Equation (1) readily shows that each charge oscillates harmonically and parallel to the atomic bonds. The incoming electric field drives these oscillations. Thus, only the charges within the optical penetration depth radiate dipolar and quadrupolar electromagnetic waves. The total amount of dipole and quadrupole far-field radiation is thus given by the summation of all electromagnetic radiation directions [24] that are produced by oscillating charges along each of the bonds. Here, the induced dipole oscillates under the action of an applied electric field and appropriate restoring and dissipation forces. The number of nearest neighbors determines the number of bonds. Therefore, they depend on the symmetry of the crystal. The nonlinear part of the motion is represented in terms of complex hyperpolarizabilities obtained by solving Equation (1). According to classical physics, an accelerated charge radiates electromagnetic waves. The total far field is calculated by summation of all the fields radiated by each dipole/quadrupole. Therefore, the model extends the Ewald–Oseen theorem to nonlinear optics. Figure 2 depicts the coordinate system used in the model.

It should be mentioned that the superpositions of all atomic dipoles driven by the external field are coherently driven dipoles, which means that the fundamental and the second/third harmonic outputs radiate in phase. The coherence is a consequence of the difference between the distance between dipoles, which is in the order of the atomic bond distance (1 nm), and the dimension of the incoming wavelength (100 nm). Hence, Huygens’ principle implies that the coherent dipole field superposition gives a constructive interference only if the outgoing angle (${\theta }_o$) is equal to the incoming angle (${\theta }_i$) as depicted in Figure 2. In transmission, a coherent superposition occurs, but this time ${\theta }_i^{\omega }\ne {\theta }_i^{2\omega }$, since the constructive angle is given by Snell’s law.

We will now derive a solution from Equation (1) with the assumption that $\mathbf{x}$ can be written as an expansion:

$$\begin{array}{c}\hfill \mathbf{x}={\mathbf{x}}_0+\Delta {\mathbf{x}}_1{e}^{-i\omega t}+\Delta {\mathbf{x}}_2{e}^{-i2\omega t}+\Delta {\mathbf{x}}_3{e}^{-i3\omega t}+\dots \end{array}$$

(2)

Inserting Equation (2) in Equation (1) the linear ${\mathbf{p}}_{1j}$, quadratic ${\mathbf{p}}_{2j}$, and cubic ${\mathbf{p}}_{3j}$ dipole polarization of the j-th charge are given by

$$\begin{array}{c}{\mathbf{p}}_{1j}=q_j\Delta {\mathbf{x}}_1{\widehat{\mathbf{b}}}_j=\frac{q_j^2{\mathbf{E}}_{loc}\left(\omega \right)·{\widehat{\mathbf{b}}}_j}{{\kappa }_1-m{\omega }^2-i\gamma \omega }{\widehat{\mathbf{b}}}_j={\alpha }_{1j}{\widehat{\mathbf{b}}}_j({\widehat{\mathbf{b}}}_j·{\mathbf{E}}_{loc}\left(\omega \right))\end{array}$$

(3)

$$\begin{array}{c}{\mathbf{p}}_{2j}=q_j\Delta {\mathbf{x}}_2{\widehat{\mathbf{b}}}_j=\frac{q_j^2{\kappa }_2\Delta {\mathbf{x}}_1^2}{{\kappa }_1-4m{\omega }^2-2i\gamma \omega }{\widehat{\mathbf{b}}}_j={\alpha }_{2j}{\widehat{\mathbf{b}}}_j{({\widehat{\mathbf{b}}}_j·{\mathbf{E}}_{loc}\left(\omega \right))}^2\end{array}$$

(4)

$$\begin{array}{c}{\mathbf{p}}_{3j}=q_j\Delta {\mathbf{x}}_3{\widehat{\mathbf{b}}}_j=\frac{q_j^2{\kappa }_3\Delta {\mathbf{x}}_2^2}{{\kappa }_1-9m{\omega }^2-3i\gamma \omega }{\widehat{\mathbf{b}}}_j={\alpha }_{3j}{\widehat{\mathbf{b}}}_j{({\widehat{\mathbf{b}}}_j·{\mathbf{E}}_{loc}\left(\omega \right))}^3\end{array}$$

(5)

here ${\alpha }_{1j}$, ${\alpha }_{2j}$, and ${\alpha }_{3j}$ are, respectively, the atomic-scale (microscopic) linear polarizability, first-order, and second-order nonlinear hyperpolarizability of the j-th bond. The hyperpolarizability parameters are usually obtained from experimental fitting. Therefore, the model is usually called phenomenological, while future improvements to incorporate ab initio quantum mechanical calculations [8] for the derivation of hyperpolarizability parameters may improve the model accuracy. The incoming s- and p-polarized field are expressed in terms of ${\mathbf{E}}_s=E_s\widehat{\mathbf{y}}$ and ${\mathbf{E}}_p=E_p(-\widehat{\mathbf{x}}cos\theta +\widehat{\mathbf{z}}sin\theta )$.

The total polarization can be readily obtained by summing over all the bonds in the conventional cell, where N is the atomic units of dipoles per unit volume (V)

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{P}=\frac{N}{V}\sum _j{\mathbf{p}}_j& =\frac{N}{V}\sum _j({\alpha }_{1j}{\widehat{\mathbf{b}}}_j\otimes {\widehat{\mathbf{b}}}_j)·{\mathbf{E}}_{loc}\left(\omega \right)\hfill \\ & +\frac{N}{V}\sum _j({\alpha }_{2j}{\widehat{\mathbf{b}}}_j\otimes {\widehat{\mathbf{b}}}_j\otimes {\widehat{\mathbf{b}}}_j)··({\mathbf{E}}_{loc}\left(\omega \right)\otimes {\mathbf{E}}_{loc}\left(\omega \right))\hfill \\ & +\frac{N}{V}\sum _j({\alpha }_{3j}{\widehat{\mathbf{b}}}_j\otimes {\widehat{\mathbf{b}}}_j\otimes {\widehat{\mathbf{b}}}_j\otimes {\widehat{\mathbf{b}}}_j)···({\mathbf{E}}_{loc}\left(\omega \right)\otimes {\mathbf{E}}_{loc}\left(\omega \right)\otimes {\mathbf{E}}_{loc}\left(\omega \right))\hfill \end{array}\end{array}$$

(6)

The relation of the atomic-scale hyperpolarizabilities in Equation (6) with the macroscopic nonlinear susceptibility is straightforward

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbf{P}=& {\overleftrightarrow{\chi }}^{\left(1\right)}·\mathbf{E}\left(\omega \right)+{\overleftrightarrow{\chi }}^{\left(2\right)}··\mathbf{E}\left(\omega \right)\otimes \mathbf{E}\left(\omega \right)+{\overleftrightarrow{\chi }}^{\left(3\right)}···\mathbf{E}\left(\omega \right)\otimes \mathbf{E}\left(\omega \right)\otimes \mathbf{E}\left(\omega \right)\hfill \end{array}\end{array}$$

(7)

where ${\overleftrightarrow{\chi }}^{\left(1\right)}$, ${\overleftrightarrow{\chi }}^{\left(2\right)}$, ${\overleftrightarrow{\chi }}^{\left(3\right)}$ are, respectively, the first, second, and third-order susceptibility tensors of the system.

In RASHG experiments, the sample is often rotated along the z-axis to produce a polar nonlinear intensity plot. In our calculation, the azimuthal angle $\varphi$ is measured from the x-axis and attached to each bond direction by applying the rotation matrix $R_z\left(\varphi \right)$:

$$\begin{array}{c}\hfill R_z\left(\varphi \right)=\left(\begin{array}{ccc}cos\varphi & -sin\varphi & 0\\ sin\varphi & cos\varphi & 0\\ 0& 0& 1\end{array}\right)\end{array}$$

(8)

This readily gives us the formula for the nonlinear SH in terms of the atomic-scale nonlinear hyperpolarizability:

$$\begin{array}{c}\hfill {\overleftrightarrow{\chi }}^{\left(2\right)}=\frac{N}{V}\sum _j{\alpha }_{2j}(R_z\left(\varphi \right)·{\widehat{\mathbf{b}}}_j)\otimes (R_z\left(\varphi \right)·{\widehat{\mathbf{b}}}_j)\otimes (R_z\left(\varphi \right)·{\widehat{\mathbf{b}}}_j)\end{array}$$

(9)

as well as the expression for the TH atomic-scale nonlinear susceptibility:

$$\begin{array}{c}\hfill {\overleftrightarrow{\chi }}^{\left(3\right)}=\frac{N}{V}\sum _j{\alpha }_{3j}(R_z\left(\varphi \right)·{\widehat{\mathbf{b}}}_j)\otimes (R_z\left(\varphi \right)·{\widehat{\mathbf{b}}}_j)\otimes (R_z\left(\varphi \right)·{\widehat{\mathbf{b}}}_j)\otimes (R_z\left(\varphi \right)·{\widehat{\mathbf{b}}}_j)\end{array}$$

(10)

Equations (3)–(5) can be applied to obtain an expression of the far field:

$$\begin{array}{c}\hfill {\mathbf{E}}_{ff}=k^2\frac{e^{ikr}}{r}\left(\sum _j{\mathbf{p}}_j-\widehat{\mathbf{k}}\left[\widehat{\mathbf{k}}·\sum _j{\mathbf{p}}_j\right]\right)=k^2\frac{e^{ikr}}{r}\left(\overleftrightarrow{\mathbf{I}}-\widehat{\mathbf{k}}\widehat{\mathbf{k}}\right)·\sum _j{\mathbf{p}}_j\end{array}$$

(11)

Here $\overleftrightarrow{\mathbf{I}}$ is the unit tensor and $\mathbf{k}=k\widehat{\mathbf{k}}$ is the outgoing wave vector pointing in the direction of the observer/detector. The derivation of Equation (11) from classical electrodynamics by applying the vector potential $\mathbf{A}(\mathbf{r},t)$ which satisfies the Lorentz gauge can be found in [22]. A possible phase difference can exist between the radiated fields of different bonds. However, it can be neglected because the light wavelength $\lambda$ is large compared to the distance between the atomic nucleus and the bond. This work considers the case where the RASHG experiments have been performed at a single wavelength. The sample is being rotated along the z-axis. The obtained nonlinear intensity can readily be obtained by multiplying Equation (11) with its complex conjugate. The considered intensity are obtained for s- and p-polarizations of the incoming fundamental and outgoing nonlinear field resulting in four combinations labeled $p(-in)$-$p(-out)$, $ps$, $sp$, and $ss$.

We focus here on the nonlinear susceptibility of the SH and in some cases TH given in Equations (9)–(10), respectively. The nonlinear susceptibility in the SBHM is obtained via bond vector summation rather than the standard group theoretical approach by assessing the crystal symmetry with regard to reflection and rotation. We describe in the following section the methods to obtain the nonlinear tensor via group theory.

3. Group Theory

Although the derivation of the higher-order tensor from the SBHM in Equations (9)–(10) is straightforward and clear, it is unfortunate that the SBHM is not widely used in the analysis of recent nonlinear optical studies. Earlier phenomenological models such as [10,11,12] apply group-theoretical analysis to derive the higher-order tensor, which as we will see, may result in unnecessary fitting parameters due to a large number of independent components in the tensor. Before we describe how the SBHM can simplify the number of independent components, it is best to briefly recall how the higher-order tensors are usually obtained using group theory (GT).

In this context, group theory (GT) can be understood as a mathematical structure that can be applied to the symmetry operations that are allowed in a crystal. The symmetry operations consist of rotations and reflections and can be represented mathematically by a matrix. Certain crystals have structural properties that may be unchanged upon specific angles of rotation and some specific existing mirror planes. An ensemble of matrices that represents all the rotations and mirror planes that keep the above-mentioned properties unchanged for a particular crystal is called a symmetry group, or more precisely called a point group [25]. Furthermore, these specific symmetry operations must preserve the crystal physical properties. This is in agreement with Neumann’s principle, or the principle of symmetry: “if a crystal is invariant with respect to certain symmetry operations, any of its physical properties must also be invariant with respect to the same symmetry operations or otherwise stated, the symmetry operations of any physical property of a crystal must include the symmetry operations of the point group of the crystal”.

To understand the matter more clearly let us recall the general macroscopic relation between the nonlinear polarization and the incoming driving field in terms of component form, which is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill P_i=& {\chi }_{ij}^{\left(1\right)}{E}_j+{\chi }_{ijk}^{\left(2\right)}{E}_j{E}_k+{\chi }_{ijkl}^{\left(3\right)}{E}_j{E}_k{E}_l\hfill \end{array}\end{array}$$

(12)

The first term includes the linear susceptibility denoted by subscript i and j. This is a 3 × 3 matrix containing nine elements. In contrast, the second is a third rank susceptibility denoted by subscripts i, j, and k containing 9 × 3 = 27 elements. The third term is the fourth rank susceptibility tensor with subscripts i, j, k and l therefore consisting of 9 × 9 = 81 components. The label i, j, k, and l refer to the three Cartesian coordinates x, y, and z. For clarity, an explicit representation of a general third rank tensor ${\chi }_{ijk}^{\left(2\right)}$ is given below:

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}{\chi }_{xxx}& {\chi }_{xyx}& {\chi }_{xzx}\\ {\chi }_{xxy}& {\chi }_{xyy}& {\chi }_{xzy}\\ {\chi }_{xxz}& {\chi }_{xyz}& {\chi }_{xzz}\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{yxx}& {\chi }_{yyx}& {\chi }_{yzx}\\ {\chi }_{yxy}& {\chi }_{yyy}& {\chi }_{yzy}\\ {\chi }_{yxz}& {\chi }_{yyz}& {\chi }_{yzz}\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{zxx}& {\chi }_{zyx}& {\chi }_{zzx}\\ {\chi }_{zxy}& {\chi }_{zyy}& {\chi }_{zzy}\\ {\chi }_{zxz}& {\chi }_{zyz}& {\chi }_{zzz}\end{array}\right)\end{array}\right)\end{array}$$

(13)

As can be seen in Equation (13) the first index “i” in the tensor ${\chi }_{ijk}^{\left(2\right)}$ is related to the rows in the main matrix. Therefore, all the elements in the first row of the inner 3 × 3 matrix have ${\chi }_{xjk}^{\left(2\right)}$ indices and for the second and third row it will be ${\chi }_{yjk}^{\left(2\right)}$ and ${\chi }_{zjk}^{\left(2\right)}$, respectively. Analogously, the indices that follow after “i”, which are “j” and “k”, refer to the usual way of labeling a 3 × 3 matrix, namely the rows and columns in the inner 3 × 3 matrix, respectively.

At first, analyzing nonlinear effects involving a third rank (or worse fourth rank) tensor looks pretty complicated due to the many elements involved. Fortunately, the symmetry of the crystal can reduce the amount of independent tensorial elements. The simplest form of the tensor that still fully represents the crystal symmetry group can be called the irreducible form of the tensor. In GT, an irreducible tensor can be obtained by performing all the symmetry operation that belongs to a certain crystal group on the most general tensor, such as in Equation (13). Therefore, in GT, rather than assessing the nonzero components using bond vectors as in Equations (9)–(10), the latter uses the conventional symmetry analysis of the crystal structure.

Now, if a symmetry operation is applied to the crystal, the third rank tensor can be described mathematically as a linear transform through the relation

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\chi }_{ijk}^{\prime }=R_{im}{R}_{jn}{R}_{ko}{\chi }_{mno}\end{array}\end{array}$$

(14)

or for the case of a fourth rank tensor

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\chi }_{ijkl}^{\prime }=R_{im}{R}_{jn}{R}_{ko}{R}_{lp}{\chi }_{mnop}\end{array}\end{array}$$

(15)

where $R_{ab}$ is a matrix defining a symmetry operation, e.g., it can be in the form of a mirror/inversion operation or general rotation of an arbitrary angle. Neuman’s principle then says that if symmetry operations belonging to the point group of the crystal are performed, then the initial and final tensor of Equations (14)–(15) must be equal, e.g., ${\chi }_{ijk}^{\prime }={\chi }_{mno}$.

For simplicity, let us calculate the resulting third rank susceptibility tensor for a crystal with a symmetry point group C_2v. Certain semiconductor surfaces can possess such symmetry. We apply Equation (14) to calculate the susceptibility third rank tensor for a crystal with C₂ symmetry (inversion symmetric under 180${}^{\circ }$ rotation). Equation (9) can be readily applied where $R_{ab}$ is now the rotational matrix given in Equation (8) and applied to the most general form of a third rank tensor in Equation (13). Contracting the matrices with the tensor to the right side of Equation (14) results in:

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}{\chi }_{xxx}& {\chi }_{xyx}& {\chi }_{xzx}\\ {\chi }_{xxy}& {\chi }_{xyy}& {\chi }_{xzy}\\ {\chi }_{xxz}& {\chi }_{xyz}& {\chi }_{xzz}\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{yxx}& {\chi }_{yyx}& {\chi }_{yzx}\\ {\chi }_{yxy}& {\chi }_{yyy}& {\chi }_{yzy}\\ {\chi }_{yxz}& {\chi }_{yyz}& {\chi }_{yzz}\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{zxx}& {\chi }_{zyx}& {\chi }_{zzx}\\ {\chi }_{zxy}& {\chi }_{zyy}& {\chi }_{zzy}\\ {\chi }_{zxz}& {\chi }_{zyz}& {\chi }_{zzz}\end{array}\right)\end{array}\right)=\left(\begin{array}{c}\left(\begin{array}{ccc}-{\chi }_{xxx}& -{\chi }_{xyx}& {\chi }_{xzx}\\ -{\chi }_{xxy}& -{\chi }_{xyy}& {\chi }_{xzy}\\ {\chi }_{xxz}& {\chi }_{xyz}& -{\chi }_{xzz}\end{array}\right)\\ \left(\begin{array}{ccc}-{\chi }_{yxx}& -{\chi }_{yyx}& {\chi }_{yzx}\\ -{\chi }_{yxy}& -{\chi }_{yyy}& {\chi }_{yzy}\\ {\chi }_{yxz}& {\chi }_{yyz}& -{\chi }_{yzz}\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{zxx}& {\chi }_{zyx}& -{\chi }_{zzx}\\ {\chi }_{zxy}& {\chi }_{zyy}& -{\chi }_{zzy}\\ -{\chi }_{zxz}& -{\chi }_{zyz}& {\chi }_{zzz}\end{array}\right)\end{array}\right)\end{array}$$

(16)

Now, for an element to be the same as its negative, the values must be zero. The final tensor thus takes the form:

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& {\chi }_{xzx}\\ 0& 0& {\chi }_{xzy}\\ {\chi }_{xxz}& {\chi }_{xyz}& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& {\chi }_{yzx}\\ 0& 0& {\chi }_{yzy}\\ {\chi }_{yxz}& {\chi }_{yyz}& 0\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{zxx}& {\chi }_{zyx}& 0\\ {\chi }_{zxy}& {\chi }_{zyy}& 0\\ 0& 0& {\chi }_{zzz}\end{array}\right)\end{array}\right)\end{array}$$

(17)

It can be seen that there are 13 nonzero components in the tensor in Equation (17). Although the form of the tensor is not similar to that in the standard GT literature [26,27], this discrepancy simply stems from the arbitrary choice of the initial coordinate system, where the latter uses the y-coordinate rotation

$$\begin{array}{c}\hfill R_y\left(\varphi \right)=\left(\begin{array}{ccc}-1& 0& 0\\ 0& 1& 0\\ 0& 0& -1\end{array}\right)\end{array}$$

(18)

rather than $R_y\left(\varphi \right)$ in Equation (8). Repeating the step for a rotation in the y-coordinate, this value produces results similar to [26] and its step-by-step derivation can be seen in [28]. In the following section, we focus our analysis on the higher-order susceptibility tensor diamond-like structures such as Si, and show that the SBHM and GT are in perfect agreement. Afterwards, we discuss different semiconductors’ crystal symmetry such as zincblende and wurtzite.

4. Third Rank Tensors of Diamond Semiconductor Structures

This section discusses the susceptibility tensor related to SHG from diamond semiconductor structures such as Si. It is well known that diamond semiconductors are inversion symmetric inside the bulk. Thus, the main contribution should come from surfaces dipoles where the symmetry is broken and from bulk quadrupoles that are always present within the bulk. In addition, very recently, it was discovered that spatial dispersion in the form of incoming field decay also contributes to the SHG signal [18].

4.1. Si(111) Surface

The first step is to determine the number of bonds representing a repetitive structure on the Si(111) surface. As can be seen in Figure 3, the atomic bond structure can be seen in terms of bonds whose arrow represents the direction of the anharmonic electron oscillation. Due to the surface symmetry, only four bonds are required to model the SHG radiation. The bond vector direction follows the coordinate system depicted in the middle of Figure 3 and is given as follows

$$\begin{array}{c}\hfill {\widehat{\mathbf{b}}}_1=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right){\widehat{\mathbf{b}}}_2=\left(\begin{array}{c}sin\beta \\ 0\\ cos\beta \end{array}\right){\widehat{\mathbf{b}}}_3=\left(\begin{array}{c}-\frac{1}{2}sin\beta \\ \frac{\sqrt{3}}{2}sin\beta \\ cos\beta \end{array}\right){\widehat{\mathbf{b}}}_4=\left(\begin{array}{c}-\frac{1}{2}sin\beta \\ -\frac{\sqrt{3}}{2}sin\beta \\ cos\beta \end{array}\right)\end{array}$$

(19)

Here, $\beta$ is, as mentioned previously, the bond angle between the Si bonds, which is 109.47${}^{\circ }$. Two different hyperpolarizabilities are required, namely the up bonds pointing normal to the surface ${\alpha }_u$ or along the z axis and the down bonds pointing downwards ${\alpha }_l$.

Applying the SBHM calculation in Equation (9), we obtain the susceptibility third rank tensor for the Si(111) surface:

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}{c}_1& c_2& c_3\\ c_2& -c_2& 0\\ c_3& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}{c}_2& -c_1& 0\\ -c_1& -c_2& c_3\\ 0& c_3& 0\end{array}\right)\\ \left(\begin{array}{ccc}{c}_3& 0& 0\\ 0& c_3& 0\\ 0& 0& c_4\end{array}\right)\end{array}\right)\end{array}$$

(20)

where

$$\begin{array}{c}{c}_1=\frac{3{\alpha }_l}{4}{sin}^3\beta cos3\varphi \end{array}$$

(21)

$$\begin{array}{c}{c}_2=\frac{3{\alpha }_l}{4}{sin}^3\beta sin3\varphi \end{array}$$

(22)

$$\begin{array}{c}{c}_3=\frac{3{\alpha }_l}{2}{sin}^2\beta cos\beta \end{array}$$

(23)

$$\begin{array}{c}{c}_4={\alpha }_u+3{\alpha }_l{cos}^3\beta \end{array}$$

(24)

From Equations (21)–(24), one can see that there are only two fitting parameters namely ${\alpha }_u$ occurring in Equation (24) and ${\alpha }_l$ occurring in Equations (21)–(24). In other words, if the SHG source is only attributed to the Si(111) dipoles at the surface, we only require two independent fitting parameters to be adjusted when modeling the experimental data.

We now compare the results with the conventional method of obtaining the tensor via the Si(111) structural symmetry as depicted in Figure 3 using group theory (GT). According to GT, the symmetry point group for a Si(111) surface is C_3v. When analyzing surfaces, the symmetry operations do not involve the vertical z-axis. Thus, the mirror planes in the z direction or rotations that transform in some way the z-coordinate are not allowed. However, if the azimuth angle $\varphi$ is fixed to some arbitrary value, the resulting susceptibility tensors in Equation (20) is generally different than those in the standard literature, e.g., [13]. The latter are calculated using mirror planes containing two Cartesian axes where the rotation axis could be arbitrarily chosen, resulting in the final tensor with different element symmetry distribution. In contrast, the SBHM reference system is usually chosen so that the unit vectors describing the bonds have minimum components to simplify the calculation.

In this work, the Si(111) surface tensor elements are fixed to $\varphi =\pi /2$. Thus, the $C_{3v}$ symmetry tensor as reported in [13] is:

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}0& -2{\chi }_{yyy}& {\chi }_{xzx}\\ -2{\chi }_{yyy}& 0& 0\\ {\chi }_{xzx}& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}-2{\chi }_{yyy}& 0& 0\\ 0& {\chi }_{yyy}& {\chi }_{xzx}\\ 0& {\chi }_{xzx}& 0\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{zxx}& 0& 0\\ 0& {\chi }_{zxx}& 0\\ 0& 0& {\chi }_{zzz}\end{array}\right)\end{array}\right)\end{array}$$

(25)

with a small difference in the constant 2 in some of the tensor elements for the $C_{3v}$ symmetry. Here, the coefficient ${\chi }_{xxy}$ satisfies ${\chi }_{xxy}={\chi }_{xyx}=-{\chi }_{yyy}$.

To compare the tensors with the SBHM, they have to be rotated by $\varphi$ using the following relation:

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {\chi }_{ijk}=R_{il}\left(\varphi \right)R_{jm}\left(\varphi \right)R_{kn}\left(\varphi \right){\chi }_{lmn}\end{array}\end{array}$$

(26)

In the case of a $C_{3v}$ symmetry, the tensor is rotated in such a way that the system of reference has a mirror plane perpendicular to the y-axis. To compare with the literature, the rotation angle should be evaluated at $\varphi =\pi /2$ with ${\chi }_{xxx}{\stackrel{\chi }{\to }}_{yyy}$. Afterwards, the general $C_{3v}$ tensor from GT is readily obtained:

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}{\chi }_{xxx}cos3\varphi & {\chi }_{xxx}sin3\varphi & {\chi }_{xzx}\\ {\chi }_{xxx}sin3\varphi & -{\chi }_{xxx}cos3\varphi & 0\\ {\chi }_{xzx}& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{xxx}sin3\varphi & -{\chi }_{xxx}cos3\varphi & 0\\ -{\chi }_{xxx}cos3\varphi & -{\chi }_{xxx}sin3\varphi & {\chi }_{xzx}\\ 0& {\chi }_{xxx}sin3\varphi & 0\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{zxx}& 0& 0\\ 0& {\chi }_{zxx}& 0\\ 0& 0& {\chi }_{zzz}\end{array}\right)\end{array}\right)\end{array}$$

(27)

Comparing Equation (27) with Equation (20), we can solve a system of equations for the unknown tensorial coefficients ${\chi }_{ijk}$ in terms of the physical constants from the SBHM as can be found in [25].

It can be seen that the tensor obtained from GT in Equation (27) is similar to the one obtained by the SBHM in Equation (20), but with an important distinction. The number of independent fitting parameters in the GT tensor is four (${\chi }_{xxx},{\chi }_{xzx},{\chi }_{zxx},{\chi }_{zzz}$) compared to only two when using the SBHM. Moreover, if we include Kleinman symmetry, the susceptibility remains unchanged when the frequencies of the two input signals and the resulting beam are permuted. Therefore the indices in the susceptibility tensor are invariant under all possible permutations. In this case, the original four independent parameters in the GT tensor are now reduced to only three. However, even with Kleinman symmetry considered, the final tensor is still more complex than in the SBHM, which requires only two. The physical interpretation of why SBHM simplifies the tensor is that it only considers bond vectors rather than the general symmetry form of the structure. Another feature is the occurrence of a threefold rotation symmetry 3$\varphi$, meaning that the generated SHG polar plot intensity should produce a three or sixfold-like pattern when analyzed using RASHG. This feature is indeed observed, as can be seen in [18].

4.2. Si(001) Surface

This procedure can be repeated to study other Si surface orientations. Let us consider the Si(001) surface as in Figure 4.

The bond vector orientation for this surface is

$$\begin{array}{c}\hfill {\widehat{\mathbf{b}}}_1=\left(\begin{array}{c}-\sqrt{2}sin\left(\beta /2\right)\\ -\sqrt{2}sin\left(\beta /2\right)\\ -cos\left(\beta /2\right)\end{array}\right){\widehat{\mathbf{b}}}_2=\left(\begin{array}{c}-\sqrt{2}sin\left(\beta /2\right)\\ -\sqrt{2}sin\left(\beta /2\right)\\ -cos\left(\beta /2\right)\end{array}\right)\end{array}$$

(28)

$$\begin{array}{c}\hfill {\widehat{\mathbf{b}}}_3=\left(\begin{array}{c}\sqrt{2}sin\left(\beta /2\right)\\ -\sqrt{2}sin\left(\beta /2\right)\\ cos\left(\beta /2\right)\end{array}\right){\widehat{\mathbf{b}}}_4=\left(\begin{array}{c}-\frac{1}{2}sin\beta \\ -\frac{\sqrt{3}}{2}sin\beta \\ cos\beta \end{array}\right)\end{array}$$

(29)

Performing the same procedure as before the SBHM third rank susceptibility tensor takes the form

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& c_5\\ 0& 0& c_6\\ c_5& c_6& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& c_6\\ 0& 0& c_7\\ c_6& c_7& 0\end{array}\right)\\ \left(\begin{array}{ccc}{c}_5& c_6& 0\\ c_6& c_7& 0\\ 0& 0& c_8\end{array}\right)\end{array}\right)\end{array}$$

(30)

where

$$\begin{array}{c}{c}_5=({\alpha }_u{cos}^2\varphi +{\alpha }_l{sin}^2\varphi )\end{array}$$

(31)

$$\begin{array}{c}{c}_6=\frac{1}{2}({\alpha }_u-{\alpha }_l)csin2\varphi \end{array}$$

(32)

$$\begin{array}{c}{c}_7=({\alpha }_u{sin}^2\varphi +{\alpha }_l{cos}^2\varphi )\end{array}$$

(33)

$$\begin{array}{c}{c}_8=2({\alpha }_u+{\alpha }_l)co{s}^2\left(\beta /2\right)\end{array}$$

(34)

As seen in Equations (31)–(34), the SBHM Si(001) surface third rank tensor consist of only two independent tensor elements namely ${\alpha }_l$ and ${\alpha }_u$. In the same fashion, the analysis using GT shows that such a structure has a $C_{2v}$ symmetry. For a rotating angle of $\varphi =0$ we have:

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& {\chi }_{xzx}\\ 0& 0& 0\\ {\chi }_{xzx}& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& 0\\ 0& 0& {\chi }_{yzy}\\ 0& {\chi }_{yzy}& 0\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{zxx}& 0& 0\\ 0& {\chi }_{zyy}& 0\\ 0& 0& {\chi }_{zzz}\end{array}\right)\end{array}\right)\end{array}$$

(35)

If the tensor is rotated as before for an arbitrary rotation angle with regard to the z-axis, then the tensor takes the form:

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& c_9\\ 0& 0& c_{10}\\ c_9& c_{10}& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& c_{10}\\ 0& 0& c_9\\ c_{10}& c_9& 0\end{array}\right)\\ \left(\begin{array}{ccc}{c}_9& c_{10}& 0\\ c_{10}& c_9& 0\\ 0& 0& c_{11}\end{array}\right)\end{array}\right)\end{array}$$

(36)

where

$$\begin{array}{c}{c}_9={\chi }_{xzx}{cos}^2\varphi +{\chi }_{yzy}{sin}^2\varphi \end{array}$$

(37)

$$\begin{array}{c}{c}_{10}=\frac{1}{2}({\chi }_{xzx}-{\chi }_{yzy}){sin}^{2\varphi }\end{array}$$

(38)

$$\begin{array}{c}{c}_{11}={\chi }_{zzz}\end{array}$$

(39)

It can be seen that the tensor associated to this point group has five independent elements (${\chi }_{xxz}={\chi }_{xzx},{\chi }_{yyz}={\chi }_{yzy},{\chi }_{zxx},{\chi }_{zyy}$, and ${\chi }_{zzz}$) according to [13]. If we now fix the angle $\varphi$ = 0 we readily obtain the tensor elements in [13]. When we apply Kleinman symmetry, these five independent parameters are directly reduced to three.

It has to be stated here that there are SHG sources other than the surface ones responsible for the total SHG intensity in Si. These sources include bulk quadrupole contribution and spatial dispersion, which are not covered in this work but are briefly mentioned here. It is well known that Si is centrosymmetric inside the bulk. Therefore, it cannot produce dipolar SHG due to parity symmetry. However, quadrupole contribution, which is weaker than dipole contribution, can still exist inside the Si bulk and compete with surface dipole due to the sheer amount of bulk layers. Recently, a SBHM model has been developed to account for the bulk quadrupole contribution whose nonlinear tensors also contain fewer independent parameters than GT. These quadrupolar contributions can be differentiated from the surface dipole contributions by applying an arbitrary input polarization in the RASHG experiment. For a more detailed discussion, the interested reader should consult [19]. Another contribution that has recently been discovered is the possibility that an incoming field decays and produces asymmetric driven oscillations within the vertical layers inside the optical penetration depth. This effect produces a unique bulk dipolar contribution called spatial dispersion, and a more elaborate discussion of how to separate this effect with surface dipoles and bulk quadrupole can be found in [18]. The main point here is that the simplification of nonlinear tensors enables a more straightforward estimation of the various SHG contributions within the bulk and surface than if GT is applied to find the nonlinear tensor.

5. Third Rank Tensors of Perovskite Semiconductor Structures

We now proceed with a more complex structure, namely perovskite, which is a semiconductor with ABO₃ symmetry. To the best of our knowledge this is the first time that the SBHM is applied to the perovskite structure. In this work, we focus on the room temperature methyl ammonium lead iodide (MAPbI₃) tetragonal structure. As before, we apply the SBHM to investigate the third rank tensor associated with SHG. The perovskite 3D structure under analysis is given by Figure 5.

As can be seen in Figure 5, the perovskite MAPbI₃ surface and bulk can be described using a different symmetry. The surface structure possesses a C_2v point group symmetry which is noncentrosymmetric whereas the bulk structure is centrosymmetric and is described by the D_4h. However, the centrosymmetry of the bulk structure is still under debate as whether such bulk structure is polar or nonpolar [29]. If the MA ions can be neglected as an SHG source as claimed by [30], then it can be assumed that the bulk structure made of Pb and I covalent bonds are centrosymmetric.

We now show using the tensor formalism of the SBHM that the neglect of MA ions indeed results in a fully zero third rank susceptibility tensor. The main challenge in obtaining the third rank tensor expression using the SBHM is determining how many bond vectors are required and how they are arranged. Afterwards, one can create an effective or reduced structure based on translational analysis of the bond vector. The MAPbI₃ bulk structure and its bond vector creation is given in Figure 6. As can be seen, we apply four Pb atoms surrounded by I atoms forming a tetragonal structure. The total involved bond vectors, in this case, is 36. Inserting all the bond vectors, one readily obtains a zero third rank tensor, confirming the centrosymmetry of the MAPbI₃ bulk structure. However, such a cumbersome calculation can be avoided if we look at the reduced form of the bond vector. Using translation to merge the vectors into a centralized point vector structure, one can see, as depicted in Figure 7 that the 26 bond vectors will exactly cancel each other out, thus producing a zero net vector. The complete expression of the perovskite bulk and surface bond vectors is avaiable as Supplementary Materials.

We now apply the same procedure to investigate the MAPbI₃ surface symmetry. The difference between the surface and bulk lies in the absence of the four up bond vector. This breaks the centrosymmetry into a C_2v structure, as can be seen in Figure 8. As before, we apply four Pb atomic centers for the SBHM model. From the reduced form of the bond vectors, it can be seen that the twofold rotational symmetry is readily obtained due to the four downward bond vectors. The number of bond vectors that are now involved consists of 40 bonds. Each of the bonds represents, as before, the covalent bond between the Pb and I atoms, neglecting the MA ions. Note that the number of bond vectors required is 40 or 4 more bonds than the previous 36 used to model the perovskite bulk structure.

The third rank susceptibility tensor from the SBHM can now be obtained. Using the same procedure by inserting all 40 bond vectors, the tensor takes the following form:

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& c_{12}\\ 0& 0& 0\\ c_{12}& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& 0\\ 0& 0& c_{13}\\ 0& c_{13}& 0\end{array}\right)\\ \left(\begin{array}{ccc}{c}_{12}& 0& 0\\ 0& c_{13}& 0\\ 0& 0& c_{14}\end{array}\right)\end{array}\right)\end{array}$$

(40)

where

$$\begin{array}{c}{c}_{12}={\alpha }_{PbIsurf}\left(2sin\frac{17\beta }{18}sin\frac{17\beta }{9}+2sin\frac{19\beta }{18}sin\frac{19\beta }{9}+\frac{3}{2}sin\frac{35\beta }{18}sin\frac{35\beta }{9}\right)\end{array}$$

(41)

$$\begin{array}{c}{c}_{13}=\frac{1}{2}{\alpha }_{PbIsurf}\left(4sin\frac{17\beta }{18}sin\frac{17\beta }{9}+4sin\frac{19\beta }{18}sin\frac{19\beta }{9}+sin\frac{35\beta }{18}sin\frac{35\beta }{9}\right)\end{array}$$

(42)

$$\begin{array}{c}{c}_{14}=4{\alpha }_{PbIsurf}\left(2cos{\left(\frac{17\beta }{18}\right)}^2+cos{\left(\frac{19\beta }{18}\right)}^2\right)+cos{\left(\frac{35\beta }{18}\right)}^2\end{array}$$

(43)

Based on Equations (41)–(43), we can conclude that the SBHM gives seven nonzero elements with only one fitting parameter, namely the effective Pb-I nonlinear hyperpolarizability (${\alpha }_{PbIsurf}$).

The C_2v third rank tensor from GT can be obtained by performing the point group rotation and mirror symmetry operation as given in the previous example of finding the C_2v tensor, but now adding to the mirror symmetry operation. If this is performed, one readily obtains:

$$\begin{array}{c}\hfill {\chi }_{ijk}^{\left(2\right)}=\left(\begin{array}{c}\left(\begin{array}{ccc}0& 0& {\chi }_{xzx}\\ 0& 0& 0\\ {\chi }_{xxz}& 0& 0\end{array}\right)\\ \left(\begin{array}{ccc}0& 0& 0\\ 0& 0& {\chi }_{yzy}\\ 0& {\chi }_{yyz}& 0\end{array}\right)\\ \left(\begin{array}{ccc}{\chi }_{zxx}& 0& 0\\ 0& {\chi }_{zyy}& 0\\ 0& 0& {\chi }_{zzz}\end{array}\right)\end{array}\right)\end{array}$$

(44)

Therefore, GT requires seven nonzero elements and as in the Si(001) case requires five independent elements (${\chi }_{xxz}={\chi }_{xzx},{\chi }_{yyz}={\chi }_{yzy},{\chi }_{zxx},{\chi }_{zyy}$, and ${\chi }_{zzz}$) but with Kleinman symmetry, this is readily reduced to three elements. Still, the SBHM is far more simple by requiring just one fitting parameter. State-of-the-art RASHG experimental results by Frohna et al. [30] show that there is no SHG contribution from within the bulk of inorganic halide tetragonal perovskite such as MAPbI₃. One of the implications is that the dynamic and large static Rashba effects, which come from within the bulk, cannot be attributed to the SHG source. Thus, the SHG must come from the $C_{2v}$ surface. Using the SBHM, we confirm their results. In addition, we found that the three independent fitting parameters in the GT SHG surface tensor of perovskite can be reduced to only one fitting parameter, namely the nonlinear effective perovskite hyperpolarizability (${\alpha }_{PbIsurf}$). This result is important because it allows the future analysis of more complex perovskite surfaces.

Before we conclude, we would like to briefly discuss the accuracy of the proposed model with experimental results, particularly current results from the SHG intensity in various RASHG experiments. As mentioned earlier, the SBHM is a phenomenological model because the complex nonlinear hyperpolarizabilities must be obtained by experimental fitting. More precisely, after the nonlinear SHG tensor form is obtained using bond vector direct products in Equation (9), the far field can be obtained through Equation (11) and, when multiplied with its complex conjugate, it readily gives the SHG intensity that can be fitted with the experiment to obtain the hyperpolarizability parameter. Because of the complexity of the various SHG sources in nonvicinal Si, an estimate of the residual fitting in vicinal Si(111) where the surface effects dominate has been offered instead by Powel et al. [13]. Applying the SBHM, they found using two independent complex parameters (the SHG up and down hyperpolarizabilities) a perfect fit with the experiment with a residual statistical error of only 0.0039 (e.g., see Figure 2 in the reference). A similar analysis of SBHM on GaAs [22] where the RASHG intensity experimental data are assumed to stem only from within the bulk also has been analyzed using the SBHM with a fitting statistical error of only 0.017, which is a good validation of the model.

6. Conclusions

We have shown that SBHM can simplify the nonlinear tensor compared to group theory. For the case of both Si(111) and Si(001) surfaces, SBHM requires only two fitting parameters namely ${\alpha }_u$ and ${\alpha }_l$ whereas GT with no Kleinman symmetry shows that four (${\chi }_{xxx},{\chi }_{xzx},{\chi }_{zxx},{\chi }_{zzz}$) and five (${\chi }_{xxz}={\chi }_{xzx},{\chi }_{yyz}={\chi }_{yzy},{\chi }_{zxx},{\chi }_{zyy}$, and ${\chi }_{zzz}$) independent parameters are required for the Si(111) and Si(001) surfaces, respectively. For the case of perovskite, the SBHM surprisingly requires only one independent fitting parameter. This parameter is the effective perovskite hyperpolarizability (${\alpha }_{PbIsurf}$). Meanwhile, GT requires five (${\chi }_{xxz}={\chi }_{xzx},{\chi }_{yyz}={\chi }_{yzy},{\chi }_{zxx},{\chi }_{zyy}$, and ${\chi }_{zzz}$) independent elements when modeling the tetragonal perovskite surfaces. Therefore, the SBHM provides a major simplification even if Kleinman symmetry is invoked.

Our results open the possibilities for future research directions. First, with the nonlinear tensor containing fewer independent parameters, one can start to estimate better the various bulk and surface nonlinear contributions in other semiconductors, paving the way for a better understanding of nonlinear processing at the nanoscale. Second, the SBHM calculation of perovskite can be developed further to produce the SHG intensity that can be compared with existing experiments, e.g., Frohna et al. [30]. For the case of perovskite, this approach seems to agree with their result that the SHG comes from the perovskite surface and that the large static and dynamic Rashba effect is not responsible for the SHG intensity. Third, real-time atomic/molecular deposition can be monitored and modeled more easily with the SBHM because we can attribute various features of the SHG intensity to single atomic bonds, including the bonds connected to the deposited atoms/molecules.