First-principles calculation of nonlinear optical responses by Wannier interpolation

Various nonlinear optical (NLO) responses, like shift current and second harmonic generation (SHG), are revealed to be closely related to topological quantities involving the Berry connection and Berry curvature. First-principles prediction of NLO responses is of great importance to fundamental research and device design, but efficient computational methods are still lacking. The main challenge is that the calculations require a very dense $k$-point sampling that is computationally expensive and a proper treatment of the gauge problem for topological quantities. Here we present a Wannier interpolation method for first-principles calculation of NLO responses, which overcomes the challenge. This method interpolates physical quantities accurately for any desired $k$ point with little computational cost and constructs a smooth gauge by the perturbation theory. To demonstrate the method, we study shift current of monolayer GeS and ${\mathrm{WS}}_2$ as well as SHG of bulk GaAs, getting good agreements with previous results. We show that the traditional sum rule method converges slowly with the number of bands, whereas the perturbation way does not. Moreover, our method is easily adapted to build tight-binding models for the following theoretical investigations. Last but not least, the method is compatible with most first-principles approaches, including density functional theory and beyond. With these advantages, Wannier interpolation is a promising method for first-principles studies of NLO phenomena.

Figure 1

(a) and (b) Atomic structures of monolayer (a) GeS and (b) ${\mathrm{WS}}_2$; top view (upper) and side view (lower). Yellow and purple balls represent S and Ge or W atoms, respectively. (c) and (d) 2D Shift current conductivities as a function light frequency $\omega$; (c) ${\sigma }_{xxx}$ and ${\sigma }_{xyy}$ together with ${\sigma }_{xxx}$ adapted from Ref. [9] for monolayer GeS and (d) ${\sigma }_{yyy}$ and ${\sigma }_{yxx}$ for monolayer ${\mathrm{WS}}_2$. The calculated band gap for GeS and ${\mathrm{WS}}_2$ is 1.90 eV and 1.82 eV, respectively.

Figure 2

Distribution of shift current related quantities in the whole Brillouin zone for monolayer ${\mathrm{WS}}_2$ and an excitation energy $\left(E_{nm}=\hslash {\omega }_{nm}\right)$ of 2.719 eV: (a) transition rate $r_{nm}^y{r}_{mn}^y\delta ({\omega }_{nm}-\omega )$ $({\AA }^2$/eV), (b) shift vector $R_{mn}^{y,y}$ (Å), and (c) their product $r_{nm}^y{r}_{mn}^y{R}_{mn}^{y,y}\delta ({\omega }_{nm}-\omega )$ $({\AA }^3$/eV).

Figure 3

(a) Band structure of monolayer ${\mathrm{WS}}_2$ by DFT calculations (black lines) and by Wannier interpolation (red dots) using two different sets of MLWFs. The left (right) panel uses a relatively small (large) set of MLWFs, corresponding to results of shift current labeled “small” (“large”) in (b) and (c). (b) and (c) Shift current calculated by (b) the perturbation method and (c) the sum rule method, based on different sets of MLWFs. (d) Shift current calculated by the perturbation method without and with applying the tight-binding approximation (discarding all off-diagonal elements of $\langle \mathbf{0}n\left|\widehat{\mathbf{r}}\right|\mathbf{R}m\rangle )$, which are labeled “full” and “tb”, respectively.

Figure 4

(a) The modulus and (b) imaginary part of the second harmonic generation susceptibility ${\chi }^{xyz}$ of bulk GaAs calculated by the Wannier interpolation method. Results based on different first-principles calculations, the “DFT-PBE + scissors correction” and DFT-HSE methods, are labeled “scissors” and “HSE”, respectively. In (a), the results are compared with those adapted from Ref. [21] that uses the “DFT-PBE + scissors correction” method without doing Wannier interpolation.