Surface second-harmonic generation (SSHG) has been shown to be an effective, non-destructive, and non-invasive probe to study surface and interface properties (Chen et al., 1981; Shen, 1989; McGilp et al., 1994; Bloembergen, 1999; Lüpke, 1999; McGilp, 1999; Downer et al., 2001a,b). SSHG spectroscopy is now very cost-effective and popular because it is an efficient method for characterizing the properties of buried interfaces and nanostructures. The high surface sensitivity of SSHG spectroscopy is due to the fact that within the dipole approximation, the bulk second-harmonic generation (SHG) in centrosymmetric materials is identically zero. The SHG process can occur only at the surface where the inversion symmetry is broken. SSHG has useful applications for studying thick thermal oxides on semiconductor surfaces (Hasselt et al., 1995; Kolthammer et al., 2005) and thin films (Yeganeh et al., 1992). The accurate determination of these studies is highly dependent on multiple reflections of both the SH and fundamental waves in the surface region. These considerations have been taken into account to study thin films (Hase et al., 1992; Buinitskaya et al., 2002, 2003) and, using the Maker-fringe technique (Maker et al., 1962), other materials (Tellier and Boisrobert, 2007; Abe et al., 2008).

Bloembergen and Pershan (1962) were the first to consider multiple reflections in their treatment of SHG in a non-linear slab. However, they only considered the second-harmonic (SH) fields and derived results for a dielectric with a small linear reflectance. They also neglected the multiple reflections of the fundamental waves inside the media. Surface effects were modeled by taking the limit of a thin slab with a thickness much smaller than the wavelength of the incoming light. Dick et al. (1985) used this methodology to determine the components of the non-linear optical susceptibility tensor, χ(−2ω; ω, ω), of a fluorescent dye over fused silica. Later works (Sipe et al., 1987; Mizrahi and Sipe, 1988) developed a simplified method using phenomenological models in which the surface is treated as an infinitesimally thin dipole sheet. The inclusion of multiple reflections is necessary for both the SH radiation and the incoming fundamental fields; this was experimentally verified in the study by Morita et al. (1988), where they show that the lineshape of the SSHG radiation is composed of resonances from both the SH and fundamental waves. Recent studies in amorphous silicon thin (and thick) films (Kessels et al., 2004; Aarts et al., 2006; Lettieri et al., 2007) have had success in characterizing the polarization and angular dependence of the SH signal produced from within the film. A very comprehensive review of SHG from films and substrate systems can be found in the study by Gielis et al. (2008).

As mentioned above, SSHG is particularly useful for studying the surfaces of centrosymmetric materials. From the theoretical point of view, the calculation of χ(−2ω; ω, ω) proceeds as follows. To mimic the semi-infinite system, we construct a supercell consisting of a finite slab of material plus a vacuum region. Both the size of the slab and the vacuum region should be such that the value of χ(−2ω; ω, ω) is well converged. A cut function is used to decouple the two halves of the supercell to obtain the value of χ(−2ω; ω, ω) for either half. If the supercell is itself centrosymmetric, the value χ(−2ω; ω, ω) for the full supercell is identically zero. Therefore, the cut function is of paramount importance to obtain a finite value for χ(−2ω; ω, ω) for either side of the slab (Reining et al., 1994; Anderson et al., 2015, 2016). The cut function can be generalized to one that is capable of obtaining the value of χ(−2ω; ω, ω) for any part of the slab. We can easily obtain the depth within the slab for which χ(−2ω; ω, ω) is non-zero; conversely, we can verify that it goes to zero toward the middle of the slab, where the centrosymmetry of the material is restored (Mejía et al., 2004). Therefore, for the surface of any centrosymmetric material, we can find the thickness of the layer where χ(−2ω; ω, ω) is finite.

On the basis of this approach for the calculation of χ(−2ω; ω, ω), in this article, we generalize the “three-layer model” for the SH radiation from the surface of a centrosymmetric material (Anderson and Mendoza, 2016). This model considers that the SH conversion takes place in a thin layer just below the surface of the material that lies under the vacuum region and above the bulk of the material. It is the three-layer model that allows us to integrate the effects of multiple reflections for both the SH and fundamental fields into the SSHG yield. As we show in this article, this treatment can be generalized to take into account the depth dependence of χ(−2ω; ω, ω) perpendicular to the surface. As shown in the study by Anderson and Mendoza (2016), the inclusion of these effects is necessary to accurately model the SSHG radiation.

This article is organized as follows. In Section 2, we present the relevant equations and theory that describe the SSHG yield that includes the depth dependence of χ(−2ω; ω, ω). In Section 3, we present calculated spectra for the Si(111)(1 × 1):H surface as a test case and contrast against experimental data. Finally, we list our conclusions and final remarks in Section 4.

In this section, we will generalize the results from our previous publication (Anderson and Mendoza, 2016) to allow for the depth dependence of χ(−2ω; ω, ω). The three-layer model proposed in the aforementioned reference considers that the surface is represented by three regions or layers. The first layer is the vacuum region (denoted by v) with a dielectric function ϵ_v(ω) = 1, from where the fundamental electric field E_v(ω) impinges on the material. The second layer is a thin layer (denoted by ℓ) of thickness d characterized by a dielectric function ϵ_ℓ(ω); it is in this layer where the SHG process takes place. The third layer is the bulk region denoted by b and characterized by ϵ_b(ω). Throughout this work, we take μ = 1. Both the vacuum and bulk layers are semi-infinite (see Figure 1). The non-linear polarization responsible for the SHG is immersed in the thin layer ℓ and is given by (1)Pℓa(z;2ω)=ϵ0χabc(z;−2ω;ω,ω)Eℓb(z;ω)Eℓc(z;ω), where χ(z; −2ω; ω, ω) is the dipolar surface non-linear depth-dependent susceptibility tensor, and the Cartesian superscripts (a, b, and c) are summed over if repeated. For ease of notation, we simply use χ(z). Also, χ^abc(z) = χ^acb(z) due to the intrinsic permutation symmetry, since SHG is degenerate in Eℓb(z;ω) and Eℓc(z;ω). We first approximate that the linear field E(z; ω) is independent of the z position. The calculation of the position dependence of the linear field is a complicated problem worth pursuing; Tancogne-Dejean et al. (2015) present a promising method for tackling this issue, which is a rather involved, non-trivial calculation. After following the method outlined in the study by Anderson and Mendoza (2016), we define the linear reflection coefficient riM as (2)riM≡riℓbeiφ1+rivℓriℓbeiφ, i=s,p, with (3)riM±=1±riM, i=s,p.

This coefficient accounts for the multiple (M) reflections of the fundamental field that depends on the thickness d of the layer ℓ included in the phase φ = 4π(d/λ₀)w_ℓ(ω), where λ₀ is the wavelength of the incoming light, w_ℓ(ω) = (ϵ_ℓ(ω)−sin²θ₀)^1/2, θ₀ is the angle of incidence, and n_ℓ = (ϵ_ℓ(ω))^1/2. The Fresnel factors tiij and riij for the vacuum-layer (ij = vℓ) and layer-bulk (ij = ℓb) interfaces are defined in equations (13) and (14) of the same reference.

To better view the effects of the z-dependence of χ(z_n) on the SHG yield, we will apply our formulation on a test surface. We choose the Si(111)1 × 1:H surface, since the (111) symmetry relations has only four non-zero components, and we can directly compare our theoretical calculations with experimental data available in the study by Mejía et al. (2002). We only present results for the p-in P-out (RpP) polarization case since it has the strongest yield and thus the best signal-to-noise ratio for the measured data. The (111) surface has only the following non-zero components of χ(z_n): χ^zzz(z_n), χ^zxx(z_n) = χ^zyy(z_n), χ^xxz(z_n) = χ^yyz(z_n), and χ^xxx(z_n) = −χ^xyy(z_n) = −χ^yyx(z_n), so we can easily work out that (11)ϒpP(zn)=TpvℓNℓtpvℓnℓ2sinθ0rpM+2sin2θ0RpM+(zn)χzzz(zn)+rpM−2wℓ2RpM+(zn)χzxx(zn)−wℓWℓ2rpM+rpM−sinθ0RpM−(zn)χxxz(zn)+rpM−2wℓRpM−(zn)χxxx(zn)cos3ϕ, where the threefold azimuthal symmetry of the SHG signal that is typical of the C_3v symmetry group is seen in the 3ϕ argument of the cosine function.

We consider that the Si(111)(1 × 1):H surface is an excellent case to test the versatility of the three-layer model; in particular, to study the effect that the z-dependence of χ^abc(z_n) and the multiple reflections will have on the SSHG yield. This surface is experimentally well characterized (Mitchell et al., 2001; Mejía et al., 2002; Bergfeld et al., 2004), and we have had success in reproducing these experimental results using the three-layer model with and without multiple reflections in our previous publications (Anderson and Mendoza, 2016; Anderson et al., 2016). The details of the ab initio calculation of χ^abc are discussed in the study by Anderson et al. (2016). We note that we apply a scissors shift of 0.7 eV to the theoretical spectra to include the effects of the electronic many-body interactions within the independent particle approach of our ab initio calculation. This 0.7 eV value allows the SH resonant peaks to acquire their corresponding energy positions and is obtained from a G₀W ₀ calculation (Li and Galli, 2010; Anderson et al., 2016).

The number of layers N for which χ^abc converges is a compromise between accuracy and the expenditure of computational time and resources. We found that N = 50, which includes 2 layers of H and 48 layers of Si, is an excellent compromise. Recall that the slab used in the calculation is centrosymmetric and that only half of the atomic layers of the slab is what actually contributes to χ^abc. In Figure 2, we show the largest component that contributes to RpP, χ^xxz(z_n), for several choices of z_n. z₁ is the layer that corresponds to the H layer. To recover the centrosymmetric environment of the (111) surface, we must add pairs of atomic Si layers, so that they include a vertical and a slanted bond of the tetrahedral unit cell corresponding to this face. This is described in detail in the study by Mejía et al. (2004). As we move from the surface toward the bulk of the system, χ^abc will decrease steadily. In the same figure, we show χ^xxz(z₂) + χ^xxz(z₃), which corresponds to the sum of the responses from the first and second Si layers. Likewise, we also include χ^xxz(z₂₄) + χ^xxz(z₂₅), which corresponds to the last two Si layers of the half-slab. We can see that the contribution from the H layer (z₁) is considerably smaller than that of the first two Si layers, as expected from the fact that the H atoms saturate the dangling bonds of the topmost Si, quenching the response (Mejía et al., 2002). The contribution from the deepest Si layers (z₂₄ and z₂₅) is quite small compared to the topmost layers. It is logical to assume that their contribution should be zero, as the slab regains its intrinsic centrosymmetry toward the center. This is not exactly the case, however, as the relatively small size of the slab yields a correct qualitative result for this regime. Figure 2 also presents two insets with the individual responses for the mentioned Si layers. The left inset depicts the spectra for the first two Si layers, at z₂ and z₃; both have a similar line shape with matching signs and thus interfere constructively creating an overall enhancement of the final intensity. On the other hand, the right inset depicts the equivalent spectra for the Si layers at z₂₄ and z₂₅. These two have a very similar line shape but with opposite signs and therefore interfere destructively and contribute very little to the total response.

For this slab, we find that (12)∑n=123χabc(zn)>>χabc(z24)+χabc(z25), and thus χ^abc is well converged. We present this comparison in Figure 3, contrasted against χhsxxz, which is the total response from all 25 layers of the half-slab. The sum of the layer responses from z₁ to z₂₃ is quite close to the total half-slab response. As before, layers z₂₄ and z₂₅ contribute very little to the overall SH spectrum. From these findings, we can establish that the thickness of the layer ℓ, where the SHG takes place is around d = 3.6 nm for N/2 = 25 active layers of SHG. These results prompt us to propose the following plausible scenario. We could use a larger value for d to achieve χ^abc(z_N/2−1) + χ^abc(z_N/2) = 0, for which we need to go to increasingly larger slabs. But to keep the computational burden reasonable, we choose N = 50 and only change the value of z_n such that d = ∑_nz_n gives the new chosen value of d. In view of equation (12), we can keep the same value for each of the χ^abc(z_n) components already calculated for N = 50. This would be equivalent to say that from equation (9), (13)∑n=123ϒiF(zn)>>ϒiF(z24)+ϒiF(z25), regardless of the actual value of z_n. We will refer to this plausible scenario as “Stretched z_n” below, where a factor of 2.7 is used to stretch z_n such that d = 2.7 × 3.6 ≈ 10 nm.

In Figure 4, we compare the theoretical results for the SSHG yield for p-in P-out, RpP, with the experimental results from the study by Mejía et al. (2002) that were taken over a SH energy range of 2.5–5 eV. We use θ = 65°, ϕ = 30°, and a broadening of σ = 0.075 eV. With ϕ = 30°, the contribution of χ^xxx from equation (11) is completely eliminated. The following scenarios are presented in the figure. First, “Nominal z_n” (solid blue line) is the layer-by-layer calculation for a layer thickness of d = 3.6 nm. This is the thickness of N/2 = 25 atomic layers with the different z_n positions obtained directly from the slab used in the full ab initio calculation. Second, “Stretched z_n” is the scenario proposed in the previous paragraph, where the z_n positions are now stretched by a factor of 2.7, but the same χ^abc(z_n) of the layered ab initio calculation are used. Finally, two “Average” curves are presented for d = 3.6 (dashed green line) and d = 10 nm (dashed yellow line) that use (14)χhsabc=∑n=1N∕2χabc(zn), which is the total response for the complete half-slab of 25 layers, along with the average value of equation (8), as proposed in the study by Anderson and Mendoza (2016), (15)R¯iM≡1d∫0dRiM(z)dz=Riℓbeiδ∕21+RivℓRiℓbeiδsinc(δ∕2).

This choice is very similar to placing χ^abc(z_n) at z_n → d/2 in equation (8), which can be interpreted as placing the non-linear polarization sheet in the middle of the thin layer ℓ. These last curves are neither depth dependent nor ab initio, since the value of d = 10 nm does not come naturally from the calculation of χ^abc(z_n). Indeed, the Average spectrum for d = 10 nm is the primary result from the study by Anderson and Mendoza (2016). On the other hand, the Nominal (d = 3.6 nm) result calculated from the theory developed in this work is completely ab initio.

The experimental spectrum shows two very well-defined resonances, which come from electronic transitions from the valence to the conduction bands around the known E₁ ~ 3.4 eV and E₂ ~ 4.3 eV critical points of bulk Si (Yu and Cardona, 2005). The theoretical results reproduce the features of the spectrum, although we see that the E₂ peak is blueshifted by around 0.3 eV. We postulate that this discrepancy is mainly due to three factors. First, the RpP spectrum tends to redshift as temperature increases (Dadap et al., 1997). Since this experiment was conducted at room temperature, our theoretical results (which consider that T = 0 K) will be blueshifted to some degree. Of course, the experimental temperature at which the spectra are measured should be taken into account in a more complete formulation. Second, equation (11) shows that RpP includes all four non-zero components of χ(z_n). In particular, χ^zzz and χ^xxz include out-of-plane incoming fields that are affected by local field effects (Tancogne-Dejean, 2015); these reveal material inhomogeneities that are far more prevalent perpendicular to the surface than in the surface plane. Therefore, we expect that these out-of-plane components will be more sensitive to the inclusion of these local field effects. As mentioned earlier, these effects are quite challenging to compute (Tancogne-Dejean et al., 2015) and are beyond the scope of this article. Third, GW transition energies are needed to accurately predict the SSHG spectrum. A Bethe–Salpeter calculation will improve the position and the amplitude of the peaks, but is far beyond our current computing capabilities. The included ab initio scissors correction produces an E₁ peak position that is similar to experiment, and we have checked there is no single scissors value that can reproduce the energy positions of both the E₁ and E₂ peaks. We consider that RpP requires the proper treatment of all three of these factors to improve the calculated spectrum.

We can clearly see that RpP for the layered calculation using d = 3.6 nm (the value from the slab) differs from the one with the stretched values of z_n that lead to d = 10 nm. These enhancements are larger for E₂ than for E₁. This can be understood from the fact that the corresponding λ₀ for E₁ is larger than that of E₂. From equations (2) and (8), we see that the phase shifts are larger for E₂ than for E₁, producing a larger enhancement of the SSHG yield at E₂ from the multiple reflections. As the phase shifts grow with d, so does the enhancement caused by the multiple reflections. We have also verified that the effects of the multiple reflections from the linear field are significantly smaller than those of the SH field. This is clear since the phase shift of equation (8) is not only a factor of 2 smaller than that of equation (2) but also w_ℓ <W_ℓ. For larger energies, such as E₂, λ₀ becomes smaller, and the multiple reflection effects become more noticeable. The selected value for d < < λ₀, which comes naturally from the ab initio calculation of χ^abc, is thus very reasonable to model a thin surface layer below the vacuum region where the non-linear SH conversion takes place. Moreover, choosing a larger value d improves the peak ratio E₂/E₁ from 1.8 (d = 3.6 nm) to 2.0 (d = 10 nm), which is closer to the experimental value of 2.8 (Anderson and Mendoza, 2016). Note that the average values obtained by using R¯pM with d = 3.6 and d = 10 nm are very similar to the depth-dependent results for the corresponding value of d. In general, this means that using χhsabc in combination with R¯iM is a reasonable first approximation of the SSHG yield.

We remark that there is a significant dielectric contrast between the thin layer ℓ and the bulk region b. As discussed in Figure 6 of the study by Mendoza et al. (2006), the layer-by-layer ϵ_ℓ(z_n;ω) of a Si(100) surface begins to resemble the bulk dielectric function as we go deeper toward the bulk of the system. However, for the atomic layers of the thin layer, ℓ, ϵ_ℓ(z_n;ω) differs substantially from ϵ_b(ω). We have carried out a similar analysis for the Si(111)1 × 1:H surface used in this work and obtained equivalent results. As ϵℓ(ω)=∑n=1N∕2ϵℓ(zn;ω)≠ϵb(ω), it follows that there are sizeable reflections from the layer–bulk interface. It is tempting to use ϵ_ℓ(z_n;ω) instead of ϵ_ℓ(ω) in the theory developed in Section 2; however, this is related to the calculation of E(z; ω), as mentioned above, is outside the scope of this work.

SSHG is a powerful tool for characterizing film/substrate systems, and we consider that the flexibility of the three-layer model provides a good theoretical foundation for calculating the SH spectra at an ab initio level. On the one hand, this revised model does away with all adjustable parameters; namely, the depth of the film d is now determined exactly from the size of the slab. This approach differs from our previous work (Anderson and Mendoza, 2016), which considered d as an adjustable parameter that can be varied to match the experimental resonance intensity. Our new ab initio formulation allows us to readily simulate the effect of film thickness on the spectrum. Calculation time will increase rapidly with the slab size, but can be computed for any film thickness with sufficient computational resources. On the other hand, we must remember that χ(z_n) is strictly a surface quantity that can be calculated on a layer-by-layer basis. Thus, it is independent of the dielectric properties that conform the bulk material (substrate). These properties obviously play an important role in the calculation of the SSHG yield, in the form of the bulk dielectric function ϵ_b(ω). If that quantity is known, then we can effectively calculate the SHG spectra of any crystalline film over any bulk substrate. Indeed, our formulation properly accounts for both the exterior surface (vacuum layer) and interior interface (layer bulk) contributions to the overall SH signal. Our model is agnostic to any surface symmetries and can effectively model passivated (like the Si(111)1 × 1:H example used here) or clean surfaces with dangling bonds. The layer-by-layer approach will even provide us with the SH contribution from the layer that includes the passivation or the dangling bonds. With these elements, we consider that the depth-dependent three-layer model with multiple reflections is a robust, versatile, and computationally efficient method that is well suited for characterizing thin/thick film systems.

We have derived a formalism to calculate the SSHG yield based on the three-layer model that accurately describes the radiating system. This treatment includes the effects of multiple reflections inside the material from both the SH and fundamental fields and also takes into account the depth variation of the second-order non-linear susceptibility χ^abc(z_n). We applied this theoretical development to calculate the p-in and P-out SSHG yield of the Si(111)(1 × 1):H surface. Our depth-dependent three-layer model reproduces key spectral features and also yields an intensity very close to experiment. We consider it an upgrade over our previous model featured in the study by Anderson and Mendoza (2016), since the thickness d of the layer in which the SHG takes place can be directly determined for the calculation and requires no other free parameters.

SA: literature review, programming, and calculations. BM: mathematical framework and general theory. Both: manuscript preparation.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.