Ab initio calculation of surface nonlinear optical response

We review recent advances in ab initio calculations of surface optical second harmonic (SH) responses of several systems: Si(001)-adsorbate surfaces, Si(111) surface, GaN(0001) surface, and Ni/Cu(001) surface. Optical second harmonic generation (SHG) has recently become a tool for studying surface electronic states. However, the interpretation of experimental SHG data has been rather difficult due to the lack of guiding theoretical predictions. For well understanding of the results of the measured SH intensity, a number of ab initio theoretical calculations have been competitively developed in these several years, and they are reviewed in this contribution. The results indicate that for simple covalent systems such as the Si surface, quantitative agreement between experiment and calculation has been obtained. [DOI: 10.1380/ejssnt.2003.57]


I. INTRODUCTION
Surface optical spectroscopy is of importance as useful probes for analyzing various surface phenomena because it is sensitive to the surface electronic structure. Optical probes generally have attractive advantages: high spectral and temporal resolution, non-destructive character, and insensitivity to ambient gas (if it is sufficiently transparent in the relevant spectral range). Thus, they have received widely practical application in the field of surface science. Linear optical reflectance methods such as surface differential reflectance (SDR) and reflectance anisotropy spectroscopy (RAS) were often used to analyze surface electronic states of semiconductor of Si and GaAs [1,2]. Surface sensitivity of these linear optical probes is based on the technique to eliminate bulk optical response. On the other hand, second order nonlinear optical processes such as second harmonic generation (SHG) have intrinsic surface sensitivity. This is because the SHG is forbidden inside the bulk of centrosymmetric materials within the dipole approximation, but it is allowed at the surface, where inversion symmetry is broken. The surface specific character of SHG is an attractive feature. From the 1970's, a large number of experimental SHG studies on various metal and semiconductor surfaces have been performed, and their experimental usefulness has been sufficiently shown [3][4][5]. However, the interpretation of the observed SHG spectra is -without any further knowledge -rather difficult. Therefore, microscopic theoretical understanding of nonlinear optical processes is needed in order to understand the results of measured SH intensity.
Let us briefly survey the theoretical studies of surface SHG. So far, there have been many theoretical frameworks or models according to the nature of each material of metal, semiconductor, and molecular crystal. In the case of the metal surface, the hydrodynamic model of electron collective motion has been used to describe the nonlinear optical response of a free-electron-like metal surface from the old time [6]. Liebsch and co-workers have performed microscopic calculations due to the timedependent local density approximation (TDLDA) for jellium metal surface [7,8], and then performed calcula- * Corresponding author: h-sano@jaist.ac.jp tions considering bulk lattice potential beyond the jellium model in order to understand the anisotropic nonlinear optical response of Al surface [9,10]. Recently, first-principle calculations of surface SH response of metals have been performed: Cu(001) [11], Fe [12], Ni(110) [13], Fe/Cu(001) [14], and Ni/Cu(001) [15].
Turning to semiconductor surfaces, nonlinear optical properties of the Si surface was successfully described by the polarizable bond model that a crystal is treated as a lattice of bonds and each bond is modeled via a point harmonic oscillator [16]. From the 1990's, more realistic calculation of the surface SH response has been performed by semi-empirical tight binding (TB) approach for several Si surfaces, including clean reconstructed, Hterminated, and As-covered ones [17][18][19][20][21][22][23]. Some of these studies referred to Cini's formalism [24]. The TB calculations correctly predicted qualitative behavior of surface SH response for different surfaces, but it gave not so good quantitative account of SHG spectra. One of the reasons for this shortcoming is that it is difficult to obtain real surface electronic states strongly depending on the change in environment by the TB calculation based only on a few parameters. Because the SH signal is generated in the very thin surface region limited to a few atomic layers, only a small error in the calculation of the surface states probably give rise to significant change of calculated SH response. Most recently, more precise and reliable ab initio calculations based on density functional theory (DFT) have been performed to obtain SH response of several semiconductor surfaces: Si(001) clean and adsorbate surface systems [25][26][27][28], a Si(111) surface [29,30], and a GaN(0001) surface [31].
There are several immature points left in the calculation method of the surface SHG as compared with that of the bulk linear optical response. However, the microscopic calculation of the surface SHG has considerably progressed due to the ab initio approach, and it is nearly reaching the stage that the surface SH response can be quantitatively predicted. For the development of the surface SHG study, it is very beneficial to make sure what current status of the surface SHG calculation is.
In this paper, we review recent theoretical study of the surface SHG, or the state-of-the-art surface SHG calculations based on the modern DFT. Excellent previous reviews on SHG calculation have been written [32,33]. Therefore, the present paper specifically focuses on ab initio calculation of the surface SHG. In section II, the computational aspect including formulation of nonlinear optical susceptibility is given. In section III, recent results of surface SHG calculation of semiconductor and metal are reviewed.

A. Nonlinear polarization and nonlinear susceptibility
The second order nonlinear polarization P NL (2ω) induced in a material by the incident electric field E(ω) is given by (1) where optical susceptibility χ (2) represents dipole allowed nonlinear optical response. Γ describes bulk nonlocal response depending on the gradient of the incident electric field, and this second term gives the bulk electric quadrupole contribution to P NL . Each term of Eq. (1) consists of surface and bulk contributions. The bulk contribution of the first term χ (2) b vanishes in the centrosymmetric materials such as Si and Ge. In this case, the surface contribution of the first term χ (2) S is dominant, because the inversion symmetry is broken at the surface. The bulk quadrupole contribution of the second term with Γ b is quite small compared with the contribution of the first term with χ (2) b of non-centrosymmetric materials, and its amplitude is of the order of that with χ (2) S . In the surface region, the amplitude of the electric field normal to the surface drastically changes. Therefore, surface contribution of the second term with Γ S is expected to be important. Because two surface contributions from χ (2) S and Γ S are indistinguishable from the viewpoint of experiment, surface nonlinear polarization P NL S is often expressed by using effective surface nonlinear susceptibility χ In general, the second order nonlinear susceptibility has 3 3 = 27 components. Some of them are zero according to the symmetry of the crystal structure. To give an example, C 3v -symmetry surface like Si(111)-(1×1) or (7×7) surfaces has only four independent non-zero components as follows: S,xyy = −χ (2) S,yxy , χ (2) S,xxz = χ (2) S,yyz , χ (2) S,zxx = χ (2) S,zyy , χ ( (2) S component possesses its own anisotropy in the surface plane. Therefore, in a certain case the amplitude ratio between the susceptibility components can be determined by measuring anisotropy of the SH response, e.g. the SH intensity as a function of the azimuthal angle of a sample surface or that as a function of the polarization angle of incident light beam [34,35].

B. Screening of the electric field
The nonlinear susceptibility is not so a convenient value to compare with the raw experimental SHG data. In order to enable direct comparison with the measured SHG data, SH light intensity should be calculated from the nonlinear susceptibility considering the dielectric screening effect on the electric field. The screening of the electric field is an important issue in the surface SHG. For calculation of the SH light intensity radiated from surface nonlinear polarization, the three-layer model shown in Fig. 1 has often been used [36]. In the three-layer model, the very thin surface layer has the nonlinear susceptibility χ (2) S . In the surface layer the surface SH polarization is generated by incident light, and consequently it emits SH light into the vacuum layer. Both input and output fields are corrected due to the dielectric screening effect. This correction coefficient, i.e. Fresnel factors, can be obtained by conventional electromagnetic calculation for light propagation in the stacked dielectric layers. In this electromagnetic calculation, dielectric functions of the surface and bulk layers are needed. Because it is usually difficult to find the surface dielectric function ε S , one often assumes that the surface dielectric function is set to be the same as the bulk dielectric function: ε S = ε b . For the more proper calculation of the screening, ε S ( = ε b ) is obtained by surface linear optical calculation [37]. Previous studies [18,29] indicate that there is no remarkably quantitative difference in SHG spectra calculated by the pure three-layer model (ε S = ε b ) and the effective two-layer model (ε S = ε b ).

C. Microscopic description of nonlinear optical susceptibility
The formulation of nonlinear optical susceptibility χ (2) is obtained by solving the density matrix equation of motion using a perturbative technique. The detailed derivation is shown in several references [38,39,32]. As the in-teraction Hamiltonian describing the interaction between the matter and the radiation, one can use two types of the form: In the present paper, representative expressions of χ (2) corresponding to each gauge are provided below.
The formulation of the surface nonlinear susceptibility in the velocity gauge is give by [28,30] Im[χ where the notations r, s, n, and l indicate electronic states of the valence (V ) or conduction (C) bands. E ns denotes the direct energy gap between the one-electron energy levels n and s. The symbol p i sn denotes the matrix element of the i-Cartesian component of the momentum operator p between states s and n: . The symbolp i sn marks the matrix element of the modified momentum operator described later. This expression must be symmetrized in the last two indices (jk) in order to comply with the intrinsic permutation symmetry of χ (2) . The real part of χ (2) is obtained by employing the Kramers-Kronig transformation.
The nonlinear susceptibility in the length gauge [40,41] has been used in the surface SHG calculation [27,31]. The latest formalism (atomic units) is given by and ω nm is the energy difference between levels n and m, and f nm ≡ f n − f m , with f i the Fermi occupation factor at zero temperature. Matrix elements r nm are evaluated through momentum matrix elements p nm as r nm = p nm /imω nm . Equation (3b) represents the purely interband contribution, and Eqs. (3c) and (3d) represent the contribution associated with the intraband motion of the electrons [40,41]. When energy eigenvalues and eigenvectors of a surface system are calculated, for modeling the surface a repeated slab construction, or a supercell approach, is usually applied. This repeated slab model consists of vacuum layer and slab layer with two surfaces. In the SHG calculation using the repeated slab model, one must use the modified momentum operator defined byp = 1 2 [S(z)p + pS(z)], where S(z) is a function which decays inside the slab: it is unity at the top layer of the slab and zero at the bottom layer. By introducing the modified momentum operator into the calculation one can avoid destructive interference of the SH polarizations possibly induced by the two surfaces of the finite slab, and thus the SH response of just one surface is obtained selectively. As the function of S(z), a continuous function or a step function centered at the middle of the slab has been used.

D. Ab initio calculation of surface electronic structure
For the calculation of the nonlinear optical susceptibility of Eqs. (2) and (3), the electron energy eigenvalues and eigenvectors are used as an input, and the most modern method to obtain them is the first-principle calculation based on DFT such as ab initio pseudopotential approach and the self-consistent full potential linearized augmented plane-wave (FLAPW) method. It is known that the FLAPW method provides better precision for optical properties than the pseudopotential approach. In DFT, exchange correlation (XC) effects have been often described by the local density approximation (LDA), and recently the generalized gradient approximation (GGA) is also used as an improved way to describing XC interaction.
Despite the limitation of Kohn-Sham theory to the ground-state properties, the band topology of valence and conduction bands are for most cases rather well repre-sented by standard DFT calculation. However, band gaps are typically smaller than experimental data. This is an important problem for calculation of optical properties involving excited states. In order to improve agreement with experimental data, quasi-particle (QP) corrections are required. The GW approximation is one of the QP corrections [42]. In Hedin's GW approximation, the spatially nonlocal and energy dependent exchangecorrelation self-energy Σ is approximated by a convolution of single particle Green's function G and dynamically screened Coulomb interaction W . The GW approximation was applied mostly to the calculation of small bulk system [42], because its calculation requires enormous computational complexity. Recently it has been used in the calculation of surface systems, and has provided good results of surface linear optical properties [43,44]. The GW approximation will be applied to the surface SHG calculation in the near future.
Instead of the GW approximation, a scissors-like operator scheme has been often employed as the QP corrections. The scissors-like operator works all the conduction bands to shift upward rigidly with the constant energy shift ∆ QP , leaving the valence ones unchanged. The value of ∆ QP is usually chosen to match the experimental gap energy. This simple treatment requires little computational complexity. However, it provides not so good optical spectra especially in the high energy region, because the scissors-like approach does not improve the bandwidth of the conduction band. Another shortcoming is that this approach does not apply for cases with half-filled metallic bands [33]. Furthermore, G. Onida et al. [45] pointed out that the scissors-like approach is successful for many bulk semiconductors, but in the case of surfaces the validity of the scissors-like operator can be questioned, since selfenergy corrections can be different for surface and bulk states, which experience a different screening. Despite of several shortcomings, the scissors-like approach will be used as expedient method.
In the calculations of surface systems, blue shift of the conduction bands sometimes occurs because of quantum size effects in the thin slab structure employed in the calculation. In this particular case, the calculated spectra can be compared with the experimental data without any further QP corrections [27,30].
In addition to the GW approximation, excitonic and local-field effects have been taken into account for bulk optical calculations, and their results precisely reproduced the experimental data [45]. However, these effects have not yet been applied to the surface SHG calculations.

A. SH responses of clean, H-, and Ge-covered Si(001) surfaces
V. I. Gavrilenko et al. performed ab initio calculations of SHG spectra of Si(001)-adsorbate systems, and indicated that the calculated results explain all essential features of the measured SHG spectra in the range 3.0 < 2 ω < 3.5 eV [25,27]. Surface structures were optimized by the Car-Parinello molecular-dynamics method based on DFT with ab initio pseudopotentials using the GGA. Their calculated clean Si(001)-(2×1) and Ge/Si(001) structures have asymmetric dimers with buckling of 0.76Å, 0.79Å(1ML Ge), and 0.805Å(2ML Ge), respectively, and corresponding monohydride surfaces have symmetric dimers. Calculations of optical properties were also based on the GGA. They modeled the surface by a slab of N Sl = 12 atomic layers, and used up to 48 k points in the irreducible 2D Brillouin zone, with energy cut off E cut = 15Ry. No QP corrections were used because quantum size effects in the slab coincidentally shifted the E 1 feature very close to the measured position.
In Fig. 2(a), the calculated SHG spectra of clean Si(001)-(2×1) (solid curve) and monohydride Si(001)-(2×1) (dashed curve), and the corresponding measured spectra (open and filled circles, respectively) are presented. For the clean surface, a pronounced E 1 peak at 3.35 eV and a weak shoulder at 3.06 eV in the calculated spectra agree well with the observed features. In addition, the quenching of these peaks by hydrogen termination is quantitatively reproduced by the ab initio calculation. This quenching feature was not explained by semi-empirical tight-binding method [19].
In Fig. 3, calculated SHG spectra over the expanded range 0 < 2 ω < 6 eV are presented. Intense peaks are seen in lower-frequency region (2 ω < 3 eV) for clean Si(001) surface. These peaks result from surface states associated with dimers and dangling bonds. This assignment is supported by the calculated result that the lowerfrequency SH peaks drastically change by H and Ge adsorption. Namely, these peaks are mostly eliminated by monohydride adsorption [Fig. 3(a), lower solid curve], because the monohydride adsorption symmetrizes and depolarizes the surface dimer and passivates dangling bonds. On the other hand, they are strengthened and blueshifted by Ge adsorption [ Fig. 3(b), upper two solid curves], because the Ge adsorption increases the asymmetry and polarization of the dimer.
A SHG peak near 3.3 eV in spectra of clean Si surface is related to E 1 bulk resonance. Since the bulk SH response of a centrosymmetric Si crystal vanishes, the peak near 3.3 eV arises not directly from the bulk bonds, but from the deformed back-bonds located near the surface. The intensity of this bulk-like peak is strongly sensitive to H and Ge adsorption as shown in Figs. 2 and 3. The most consistent explanation of the origin of this high sensitivity is that adsorbates chemically influence underlying bonds via orbital rehybridization, or charge transfer. The calculated result of Fig. 9 of Ref. [27] (not shown here) supports this mechanism, because it indicates that charge density around back-bonds near the surface is significantly changed by H adsorption.
V. I. Gavrilenko et al. [27] reported interesting facts in the SHG calculation of Si(001) surfaces as follows. They carefully analyzed convergence of SHG and reflectance difference spectroscopy (RDS) spectra with respect to slab thickness N Sl , and found that the SHG spectra converged very well with N Sl = 8, while the low-frequency ( ω ≥ 3 eV) RDS converged within 10-15% only for N Sl ≥ 12, and the higher-frequency RDS peaks at 4.3 and 5.3 eV continued to change significantly with increasing N Sl . The faster convergence of SHG with N Sl results from its stronger surface localization in covalent cubic materials.

B. SH responses of H-terminated and clean 2×1-reconstructed Si(111) surfaces
The monohydride Si(111)-(1×1) surface is the most simple and characterized of all semiconductor surfaces. Therefore, it represents an ideal reference system for SHG studies of semiconductor surfaces. J. E. Mejia et al. have performed ab initio SHG calculation of the H-terminated Si(111)-(1×1) surface by using the DFT-LDA with pseudopotentials, and obtained a semi-quantitative agreement between theory and experiment [29].
The surface was modeled by a slab of N Sl = 36 atomic layers. They used a plane wave basis set with an energy cutoff E cut = 15Ry, and the QP correction within the scissors-like approach with ∆ QP = 0.5 eV. In order to match the experimental condition, an angle of incidence θ = 65 • and an azimuthal angle φ = 30 • were taken in the calculation, where φ is defined as an angle between the plane of incidence and the [112] direction.
In Fig. 4, calculated SHG spectra (solid curves) and measured data (crosses) are presented. The experimental spectra in Fig. 4 show two peaks originating from two photon resonances of the E 1 (3.4 eV) and E 2 (4.3 eV) transitions of bulk Si. Below E 1 the SHG yield is zero, since there are no surface related states. The theoretical curves have shapes similar to the experimental spectra: especially in Fig. 4(a), energy separation and intensity ratio between E 1 and E 2 peaks agree well with the exper-iment. However, significant difference of peak positions between theory and experiment is seen in Fig. 4. This result indicates that QP corrections within the scissors-like approach do not give an adequate reproduction of the experimental SHG spectra, and that these corrections affect differently the nonzero components of χ (2) S . Gavrilenko et al. [29] also showed that the intensity and shape of the E 1 peak significantly change by variation of 30% of a vertical position of the second-layer Si atoms, which is equivalent to 3-4% change in bond length of the Si-Si back bond, but changing the Si-H bond length by as much as 5% of its equilibrium value does not appreciably alter the SHG spectra (not shown). Judging from this calculated result, the bulk-like E 1 peak of the monohydride Si(111) surface probably results from back-bonds located near the surface.
Another ab initio SHG calculation of Si(111) surface by using the self-consistent FLAPW method within the LDA was performed by H. Sano et al. [30] They calculated the SH light intensity from the reconstructed Si(111)-(2×1) surface as a function of the polarization angle of the incident light, and have found a good agreement between the calculated and the experimental results [34].
In their calculation, the number of basis functions was restricted by momentum cutoff of 3.85 a.u. −1 for the (2×1) surface. Density, potential, and basis functions inside the atomic spheres were expanded into spherical harmonics up to l max = 8. For modeling the (2×1) surface, they used a slab of 14 Si atomic layers, and adopted the geometrical structure of the π-bonded chain model [46]. The QP correction was not used. The reconstructed Si(111)-(2×1) surface is well known to have optical absorption at ω = 0.5 eV due to optical transitions between occupied and unoccupied surface states. In order to show the reliability of the calculations, they confirmed that the calculated linear optical response reproduces the above feature well (not shown).
The SH light intensity for the Si(111)-(2×1) surface as a function of the polarization angle of the normally incident pump beam having a photon energy of 1.17 eV measured by Heinz et al. [34] is shown in Fig. 5(a). I x and I y are the SHG signal polarized along x: [211] and y: [011] directions, respectively, and I total is the total SHG signal. Under their experimental conditions [34], only three surface nonlinear susceptibility components χ (2) S,xxx , χ (2) S,xyy , and χ (2) S,yxy are non-vanishing. Using these susceptibility components, the SH intensity is written as where θ is the polarization angle of the incident light with respect to the x direction, and A is a coefficient including Fresnel factors for the incident and output fields. Heinz et al. [34] carried out a least-square fit of the calculated val- ues of Eqs. (4a)-(4c) to the measured SH patterns in Fig.  5(a), and obtained the ratio of the three susceptibility components as | χ (2) xxx |:| χ (2) xyy |:| χ These calculated patterns indicate that the anisotropy of the SH response in the surface plane depend drastically on the incident photon energy. The angular pattern of the SH response reflects the shapes of wave functions involved in the relevant optical transitions. Physically, it reflects the angular dependence of nonlinear mobility of surface electrons by the incident electric field. Since calculated results based on DFT often show energy shifts, the experimental data of Fig. 5(a) at ω = 1.17 eV was compared with the calculated ones in a sufficiently broad energy range including 1.17 eV. In the energy range from 0.15 to 4.0 eV the calculated patterns agreed well with the measured ones only at ω = 1.01 eV according to Fig. 5(c). At this energy, the ratio of the calculated susceptibility components is | χ (2) xxx |:| χ (2) xyy |:| χ (2) yxy |= 1 : 0.66 : 0.87, and this is also in reasonable agreement with the measured ratio.
Since the incident photon energy of 1.01 eV in Fig. 5(c) is close to the experiment value of 1.17 eV, their calculation reproduces the measured SH response of Si(111)-(2×1) surface very well. The shift of 0.16 eV in the incident photon energy is probably due to the fundamental error of ground-state DFT.
According to Eqs. (4a)-(4c) the phase of the complex nonlinear susceptibility components should affect the shapes of the polarization angle pattern enormously. Thus the agreement between the calculated and measured polarization angle patterns indicates the correctness of the phase of their calculated complex susceptibility tensor components.

C. SH response of GaN(0001) surface
GaN is a promising semiconductor material for various technological applications such as optical devices. Since it is a polar wurtzite crystal, the bulk SH response is allowed. In this case, it is an interesting issue to predict the surface specific features of SHG spectra of GaN crystal. V. I. Gavrilenko and R. Q. Wu calculated SHG spectra of the Ga-and N-terminated GaN(0001) surfaces by using the FLAPW method [31]. The FLAPW calculations have provided very accurate results of linear and nonlinear optical responses for A III B V compound semiconductors [47].
The surface was modeled by a slab containing eight GaN monolayers. The calculations for SH response were based on the GGA. No QP correction was adopted. They calculated electronic band structures of N-and Gaterminated GaN(0001)-(1×1) surfaces. The results (not shown) indicate that the surface states of the N-and Gaterminated surfaces are significantly different from each other. The N-terminated surface has a well pronounced low dispersive metallic surface state in the gap, while the Ga-terminated surface has three strongly dispersive surface states, indicating the strong interaction between the neighboring surface atoms.
Absolute value of the nonlinear susceptibility component | χ (2) zxx | of the GaN(0001) surface is shown in Fig.  6(a), and the SHG efficiencies for the p-in/p-out and the s-in/p-out polarization configurations are shown in Fig.  6(b) and (c), respectively. The solid curves and the gray dashed curves in Fig. 6 indicate the calculated results of Ga-terminated GaN(0001) and N-terminated GaN(0001) surfaces, respectively. The thin dashed curves in Fig. 6 indicate the rescaled spectra (by factors of 0.01 and 1/25 for the susceptibility component and the SHG efficiencies, respectively) of the N-terminated surface. The bold gray curve in Fig. 6(a) indicates the calculated susceptibility component of the bulk wurtzite GaN [47]. Let us separate the SHG spectra into two spectral regions: region 1 (0 ≤ ω ≤ 5.0 eV) and region 2 ( ω ≥ 5.0 eV). In region 2, SHG spectra of both surfaces are roughly similar to each other. Furthermore in region 2 calculated susceptibility components of both surfaces and the bulk GaN crystal are also similar to one another in magnitude as shown in Fig. 6(a). Thus, SHG spectra in region 2 probably result from the bulk SH response. On the other hand, in region 1, calculated SHG spectra appear to be dominated by contributions from the surface states because the surface susceptibility components are significantly larger than the bulk component in region 1. Fig. 6 indicates that surface specific SH responses are much stronger than the bulk SH responses. Especially N-terminated surface has highly intensive SH peaks in the low-energy region. These N-related peaks are reduced by adsorption of a Ga monolayer. From this result, SH intensity in the low-energy region can be used as a sensitive in situ probe to monitor the changes in chemical composition on the surface during the growth processes.

D.
SH response of ferromagnetic Ni/Cu(001) surface M. N. Dadoenkova et al. [15] calculated the surface SH response of the magnetized Ni bilayer on the Cu(001) surface by using the FLAPW method implemented in the package WIEN97 [48] on the basis of DFT within the GGA. For modeling the surface, a repeated slab construction (supercell) with the vacuum thickness equivalent to 10 Cu-lattice spacing was used. In the calculation of the nonlinear susceptibility, only the muffin-tin spheres of the Ni bilayer and the interstitial region containing vacuum and vacuum/slab interface were taken into account [14]. The incident angle and electric field strength of the incident light beam were taken to be 45 • and 10 8 V/m, respectively. The plane of incidence was set to be parallel to x ([100]) direction.
The SHG spectra of Ni/Cu(001) bilayer for the pinput/p-output and p-input/s-output polarization configurations are shown in Figs. 7(a) and (b), respectively. The solid and dashed curves correspond to the in-plane (M//x) and the out-of-plane (M//z) magnetizations, respectively. Fig. 7 indicates that the change of the magnetization direction from the out-of-plane to in-plane direction considerably modifies the SHG spectra: for the pin/s-out configuration, SH intensity grows and new peak appears in the energy range from ω = 7 to 7.5 eV. In the SHG spectra low-energy edge of the SH signal is seen around 2.5 eV, though the calculated nonlinear susceptibility components (not shown) have large amplitude in the energy range from ω = 1 to 3 eV. It is suggested that this edge mostly results from the onset of the d-electron absorption of Cu substrate.
The insets of Fig. 7 are the azimuthal angle dependence of the SH intensity of the Ni/Cu(001) bilayer with in-plane magnetization at the incident photon energy of ω = 1.5 eV. For p-polarized SH signal, the azimuthal angle patterns of M//x (solid curve) and M//-x (dashed curve) are obviously different in their shapes. In contrast, s-polarized SH signal shows no change of the pattern by changing the magnetization from M//x to M//-x.
Their calculated results demonstrate that SH response has high sensitivity to the surface magnetism. However, their theoretical prediction is expected to be verified by corresponding experiments.
Theoretical studies reviewed in Chapter III indicate that the surface SH response is considerably sensitive to the change of the surface electronic states localized within a few atomic layers. Therefore, the most important factor of reliable surface SHG calculation is the accuracy of the calculation of surface electronic states. From this viewpoint, ab initio calculation is very promising because it has been proved to provide the most correct surface electronic states among various methods for various surface systems. Success in the SHG calculation of Si surfaces means that the DFT-pseudopotential method can sufficiently describe Si surface electronic structure. However, for dealing with compound semiconductor surfaces and transition metal oxide surfaces, the FLAPW method is more suitable, since it can calculate their surface electronic states with a higher accuracy.

IV. SUMMARY
Recent studies on ab initio calculations of surface SH response are reviewed. For the H and Ge adsorbed Si(001) surfaces and H-terminated Si(111) surface, ab initio calculations well reproduced the features of the bulk-like E 1 and E 2 transition peaks in the measured SH spectra. For the clean reconstructed Si(111)-(2×1) surface, the calculated SH intensity as a function of the polarization angle of the incident light agreed with the measured data. These results indicate that for simple covalent systems such as Si, surface SH response can be semi-quantitatively predicted by the ab initio calculations.
The SHG calculation of the Ga-and N-terminated GaN(0001) surfaces indicated that two surfaces have different SH responses in the energy region below ω = 5 eV. For the ferromagnetic Ni/Cu(001) surface, the calculation indicated that surface SH response is very sensitive to the magnetization direction. Though these two theoretical studies present very interesting results, the validity of the calculations has not been sufficiently confirmed by comparing them with corresponding experiments.
Because the ab initio SHG calculations for various surface systems rapidly progress due to an increase of the computer power and an improvement of approximations in the calculations, the combined study of experiment and calculation will become a more powerful method to ana-lyze the surface electronic states. Therefore, accurate ab initio SHG calculations are getting more important, and are expected to be used more often in the future.