Full ab initio band structure analysis of interband and intraband contributions for third harmonic generation coefficient of bulk silicon: implementation and application of the sum-over-states

We fully implement the Aversa and Sipe sum-over-states formulism and make a full ab initio band structure analysis of interband and intraband contributions for the third-order nonlinear optical susceptibilities of bulk silicon. The band structure and momentum matrix elements were calculated by using the highly accurate all-electron full potential linearized augmented plane wave method within the local density approximation. The convergence tests including the scissor correction with different k-points meshes and empty states were performed. Both real and imaginary parts of susceptibility were directly calculated and checked by the Kramers-Kronig relation. The converged results are compared with other theoretical and experimental ones and in agreement with the recent ab initio real-time-based calculation. The nonlinear optical coefficient comes from three parts: the pure interband contribution (Pinter), the modulation of interband terms by intraband terms (Pmod), and the intraband contribution (Jintra). For each part, the origin of enhanced peaks is explored by tracing the sum-over-states process. The interband contribution is found to be dramatically modulated by the intraband contribution.


Introduction
In the last three decades, the perturbation-theory-based sum-over-states (SOS) method 1 has been widely used to calculate the optical polarizability of isolated atoms or molecules [2][3][4] . Using the SOS method, one can not only carry out the frequency-dependent optical response calculation but also investigate the electronic origin of the optical response. T he form er produces the theoretical results that more directly compared to experiments, as the nonlinear optical measurements are performed at different optical frequencies. The latter is helpful in identifying which excited states play a significant role in the optical response, then analyzing the charge transition contribution to these selective excited states, and ultimately identifying which functional groups dominate the optical response of the whole molecule. The corresponding findings will guide us to design and synthesize the material with a large nonlinear optical response 2,3 .
For condensed semiconductor materials, the SOS method has also been developed to determine their linear and nonlinear optical susceptibilities [5][6][7][8][9][10][11] . For the linear optical susceptibility χ (1) , we can easily implement the SOS calculations providing t h a t t h e b a n d s t r u c t u r e a n d momentum matrix elements are obtained. For the nonlinear optical susceptibilities such as χ (2) and χ (3) , however, the difficulty increases rapidly owing to the complexity of the equations 6,7,[9][10][11][12][13][14][15] . For χ (2) , the theoretic technique has been developed to a rather sophisticated level for both static and dynamic cases. At the outset, the static and dynamic calculations were separately developed because merging the static (ω→0) calculation into the dynamic one was hindered by apparently diverging terms in the SOS equations 13,14,16,17 (i.e., factors of ω -1 and ω -2 ). To overcome the difficulty of unphysical divergences, Sipe and Ghahramani (SG) 10 developed the formalism for calculating the nonlinear optical coefficients within the independent-particle approximation. SG eliminated the unphysical divergences by a carefully separate treatment of interband and intraband motion and provided the detailed expressions for the calculation of second harmonic generation.
For χ (3) , application of the SOS technique is much less than for χ (2) because of the complexity of equations. The connection between static and dynamic calculation was also plagued by the apparently diverging term 12,15 . SG has pointed that the divergence-free expressions of χ (3) can be developed in a similar way to χ (2) but the derivation will be a formidable task. On the other hand, also within the independent-particle approximation, by applying the perturbation theory to the dynamical equation of the electronic density o p e r ator and using a so-called length-gauge formulation, Aversa and Sipe (AS) 5 presented well-behaved, general expressions for χ (2) and χ (3) for arbitrary frequency mixings in a simpler way than SG. The expressions for χ (2) are in well consistent with those shown in SG. Significantly, the derivation of expressions for χ (3) from AS requires much fewer efforts than from SG. The expressions for χ (3) were given explicitly in the divergence-free form. Clearly, these expressions are also helpful in understanding the third-order nonlinear optical response of semiconductors. However, to our best knowledge there has been no implementation of χ (3) based on these expressions so far. Although there are a few reports 12,15,[22][23][24][25] about calculations of χ (3) based on the SOS method, these works are based on the two-band or empirical tight-banding and semi-ab initio b a n d m o d e l s . The more accurate band model calculations are required for the calculation of χ (3) which is highly sensitive to the details of the band structures and wave functions 15 .
The purpose of this work is two-fold. One is to fully implement the AS SOS formulism 5 to calculate the χ (3) . The other is to perform a full ab initio band structure analysis of interband and intraband contributions for third-order nonlinear optical susceptibilities of bulk silicon. The band structure is calculated by using the highly accurate all-electron full potential linearized augmented plane wave (FP-LAPW) method [26][27][28] within the local density approximation. The calculated susceptibilities are in agreement with the recent ab initio real-time-based computational approach combined with the Berry-phase formulation of the dynamical polarization 29 . By tracing the SOS process, we easily e x t r a c t t h e transition term w i t h a s i g n i f i c a n t c o n t r i b ution to nonlinear coefficients and identify the ω, 2ω, or 3ω r e s o n a n t c o n t r i b u t i o n s t o the peak in the frequency-dependent nonlinear optical spectra. Finally, we summarize our results in Sec. 5.

Formulism and implementation
The equations used to calculate the third-order response function were originally obtained by AS 5 w ho used the len gth-gauge formalism based on the position operator r (Ref. 30 ). For convenience of reading, we inherit the AS's notations. The third-order susceptibility tensor represented by χ (3) dcba (-ω 3 ;ω γ ,ω β ,ω α ), where ω 3 =ω γ +ω β +ω α , is decomposed to χ (3) χ and χ (3) σ based on the decomposition of the physical contributions to the polarization 5 , namely, dP/dt = dP χ /dt + J σ . The final expressions for χ (3) χ and χ (3) σ are as follows: where C = e 4 K/ћ 3 with the K factor depending on the particular combination 1 of ω γ , ω β , and ω α , for example, K is 1/4 for the third harmonic generation (THG) polarizability χ (3) dcba (-3ω;ω,ω,ω), ω 1 and ω 2 are defined by ω 1 = ω α and ω 2 = ω α + ω β , respectively. ω mn = ω m -ω n is the energy difference between the bands m and n, f mn = f m -f n is the difference of the Fermi distribution functions, the indices of a, b, and c are Cartesian directions, and all four band indices l, m, n, p are different (one exception is shown below) because r mn [= p mn /(imω mn )] is defined to be zero unless n ≠ m (Ref. 5 ). So, to calculate the χ (3) dcba by using Eqs.1 and 2, we first have to obtain the band structure of periodic system and the momentum matrix elements.
There is no obvious divergence in Eqs.1 and 2 except that the first summation of Eq.1 shows an apparent divergence arising from both a lack of r lm and r np elements in numerators and a factor of 1/ω 2 for χ (3) dcba (0;0,0,0) when l = m and n = p. Introducing intrinsic permutation symmetry and relabeling indices, we reduced these two troublesome terms (TTT) to In evaluation of the matrix elements r mn as p mn /(imω mn ), a numerical problem occurs when the bands m and n are nearly degenerate. As mentioned by SG, one can always choose the appropriate wave functions for the bands m and n such that the matrix elements r mn (or p mn ) vanishes. Therefore, for nearly degenerate bands m and n decided by a small cutoff value, e.g., ωmn ≤ 0.001 a.u., we set r mn to be zero, in consistent with the definition of r mn above. This strategy was smoothly used by Rashkeev et al. 7 in their calculations of frequency-dependent second-order optical response of semiconductors.

Computational details
We applied our implementation to the calculation of THG of cubic silicon crystal (Si) with a lattice parameter of 5.43Å. The wave functions and momentum matrix elements were computed with the highly accurate all-electron FP-LAPW method [26][27][28] within the local density approximation as implemented in the ELK code 31 . Since the LDA calculation underestimates the band gap of Si, we applied the widely used scissor correction 32 in the optical calculation. Both the real and imaginary parts of χ (3) were directly calculated and checked by the Kramers-Kronig relation (KKR) 33 . The maximum angular momentum used for APW functions is l max = 8. Since the nonlinear optical calculation possibly requires much denser k-points mesh and more empty states than the linear optical calculation 6,9 , we performed the convergence tests of χ (3) on the 10×10×10, 20×20×20, 30×30×30, and 40×40×40 k-points meshes and the number of empty states (10, 14, and 18) per atom. In terms of the limit of cubic symmetry, we only calculated two nonzero independent elements of χ (3) , namely, χ (3) 1111 and χ (3) 1212 .

Scissor correction
It is well known that the LDA underestimates the band gap in semiconductors. Since the denominators of Eqs.1 and 2 depend on the 1/ω 4 nm like factors, the underestimation of the band gaps definitely leads to an error in calculation of χ (3) . The simple and effective way is to introduce a so-called scissor correction, in which the band energies are shifted by a factor of Δω and the momentum matrix elements are corrected by p mn = p mn (1+f mn Δω/ω nm ) (Ref. 32 ). The scissor correction is derived by adding to the LDA Hamiltonian the scissor operator that is an effective self-energy 34,35 . The f mn Δω/ω nm is considered as a nonlocal contribution to the matrix element 34,35 , in consistent with a point of view of the nonlocal exact density functional 36 . This nonlocal correction may have a significant contribution to the matrix element, as shown by Nastol et al. 32 .
Their estimation for GaAs had shown that the Δω/ω nm value is about 4.4 for a lowest conduction band n, and a highest valence band m near the Γ point. As shown in calculations of χ (2) of GaAs and GaP, the magnitude of χ (2) was dramatically improved by applying the scissor correction 32 .
Using the FP-LAPW/LDA method, we obtain for Si an indirect gap of 0.46 eV which is lower by 0.69 eV than the experimental value of 1.15 eV. In the following calculations of χ (3) we used a scissor correction of 0.69 eV.

Convergence test
The calculations of χ (2) 7,9 have shown that a denser k-points mesh is required for nonlinear  Fig. 1b and 1c the convergence tests on the real parts of χ (3) 1111χ (ω) (Eq.1) and χ (3) 1111σ (ω) (Eq.2), and observe that a poor convergence of χ (3) 1111 (ω) in the low applied frequency region should be attributed to the poor convergence of χ (3) 1111σ (ω). In section 4.5, we will show that the existence of (f n /ω mn )-like terms in χ (3) 1111σ (ω) leads to a possible resonance in the low applied frequencies. In addition to the k-points convergence test, we have also run the tests on the number of empty states included in the SOS calculations. Figure 1d shows a good convergence behavior for the number of empty states (i.e., 10, 14, and 18) per atom. For , we obtain a similar behavior for convergence tests. In the following section, we will use the 30×30×30 results to make further discussions.

KKR test
The KKR 33 As an example, we shows in Fig. 2

Comparisons with other theoretical results and experiments
To get a better evaluation for our results, we compare our results with other theoretical reports 15 Figure 4 shows the dispersions of interband and intraband contributions of χ (3) (ω). In terms of the decomposition of position operator (r = r i + r e , r i and r e are intraband and interband parts of r) 5 , Eq.1 related to the interband part r e is referred to as the interband contribution while Eq.2 related to the intraband part r i is the intraband contribution (Jintra in Fig.4). Note that the later three summations of Eq.1 contain diagonal momentum matrix elements (p nn of Eq.4). So the later three summations of Eq.1 is considered as a modulation of interband terms by intraband terms (Pmod in Fig.4), similar to the decomposition of χ (2) (ω) 6,8,9 , and the first summation of Eq.1 that does not contains diagonal momentum matrix elements is considered as the pure interband contribution (Pinter in Fig.4).

Interband and Intraband contributions
Firstly, as shown in Fig.4, Pinter presents the enhancement at some applied frequencies. For example, the Re[χ (3),Pinter 1111 (ω)] presents two peaks at ω = 1.55 and 2.16 eV. To understand the origin of these two peaks, we traced the SOS calculations in terms of Eq.1 for these two applied frequencies. Figure 5 shows are given for clarity. From insets, we can see that a few k-points (labeled by a, b, c1, c2, d1, d2, e1,   e2, f1, and f2) have a relative large contribution to Re[χ (3),Pinter 1111 (ω)]. To check whether these relative large contributions are attributed to the possible ω, 2ω, or 3ω resonance in Pinter of Eq.1, we separately traced the SOS calculations for these k-points.

Conclusions
We        f This is a relative value measured relative to LiF at the same applied field (ω = 1.16 eV). Table 2. Selected contributions labeled in insets of Fig. 6 to Re[χ (3),Pinter 1111 (ω)] (a.u.). The single-particle channels with significant contributions to each k-points a r e g i v e n . 2 ( 0 .