Non-linear optics from an ab-initio approach by means of the dynamical Berry-Phase:application to second and third harmonic generation in semiconductors

Non-linear optics from an ab-initio approach by means of the dynamical Berry-Phase:
application to second and third harmonic generation in semiconductors

C. Attaccalite CNRS/Univ. Grenoble Alpes, Institut Néel, F-38042 Grenoble, France    M. Grüning School of Mathematics and Physics, Queen’s University Belfast, Belfast BT7 1NN, Northern Ireland, UK Centre for Computational Physics and Physics Department, University of Coimbra, Rua Larga, 3004-516 Coimbra, Portugal

We present an ab-initio real-time based computational approach to study nonlinear optical properties in Condensed Matter systems that is especially suitable for crystalline solids and periodic nanostructures. The equation of motions and the coupling of the electrons with the external electric field are derived from the Berry phase formulation of the dynamical polarization [Souza et al. Phys. Rev. B 69, 085106 (2004)]. Many-body effects are introduced by adding single-particle operators to the independent-particle Hamiltonian. We add a Hartree operator to account for crystal local effects and a scissor operator to correct the independent particle band structure for quasiparticle effects. We also discuss the possibility of accurately treating excitonic effects by adding a screened Hartree-Fock self-energy operator. The approach is validated by calculating the second-harmonic generation of SiC and AlAs bulk semiconductors: an excellent agreement is obtained with existing ab-initio calculations from response theory in frequency domain [Luppi et al. Phys. Rev. B 82, 235201 (2010)]. We finally show applications to the second-harmonic generation of CdTe and the third-harmonic generation of Si.

I Introduction

Ab-initio approaches based on Green’s function theory became a standard tool for quantitative and predictive calculations of linear response optical properties in Condensed Matter. In particular, the state-of-the-art approach combines the approximation for the quasi-particle band structure Aryasetiawan and Gunnarsson (1998) with the Bethe-Salpeter equation in static ladder approximation for the response function. Strinati (1988) This approach proved to effectively and accurately account for the essential effects beyond independent particle approximation (IPA) in a wide range of electronic systems, including extended systems with strong excitonic effects. Onida et al. (2002)

In contrast, for nonlinear optics ab-initio calculations of extended systems rely in large part on the IPASipe and Ghahramani (1993) with correlation effects entering at most as a rigid shift of the conduction energy levels. Cabellos et al. (2009) Within time-dependent density-functional theory (TDDFT), it has been recently proposed Luppi et al. (2010) an approach to calculate the second-harmonic generation (SHG) in semiconductors that takes into account as well crystal local-field and excitonic effects. However, this promising approach Cazzanelli et al. (2012) is limited by the treatment of the electron correlation to systems with weakly bound excitons. Botti et al. ()

Within Green’s function theory the inclusion of many-body effects into the expression for the nonlinear optical susceptibilities is extremely difficult. Furthermore the complexity of these expressions grows with the perturbation order. Therefore it is not surprising that there have been only few isolated attempts to calculate second-order optical susceptibility using the Bethe-Salpeter equation Leitsmann et al. (2005); Chang et al. (2001) and no attempt to calculate higher-order optical susceptibilities. Virk and Sipe (2009)

Alternatively to the frequency-domain response-based approach, one can obtain the nonlinear optical susceptibility in time-domain from the dynamical polarization of the system by using the expansion of in power of the applied field


This strategy is followed in several real-time implementations of TD-DFT Yabana and Bertsch (1996). In these approaches the dynamical polarization is obtained by numerical integration of the equations of motion (EOMs) for the Kohn-Sham system. Takimoto et al. (2007); Castro et al. (2004); Meng and Kaxiras (2008) So far applications regard mostly nonlinear optical properties in molecules.

The time-domain approach presents three major advantages with respect to frequency-domain response-based approaches. First, many-body effects are included easily by adding the corresponding operator to the effective Hamiltonian. Second, it is not perturbative in the external fields and therefore it treats optical susceptibilities at any order without increasing the computational cost and with the only limitation dictated by the machine precision. Third, several non-linear phenomena and thus spectroscopic techniques are described by the same EOMs. For instance, by the superposition of several laser fields one can simulate sum- and difference-frequency harmonic generation, or four-waves mixing. Boyd (2008)

In a recent work, Attaccalite et al. (2011) we proposed a real-time implementation of the Bethe-Salpeter equation, based on nonequilibrium Green’s function formalism. However, due to the problems in defining the position operator and thus , it is not trivial to apply Eq. (1) to systems in which periodic boundary conditions (PBC) are imposed. As it was recognised for example in Ref. Aversa and Sipe, 1995, the same problem appears in the direct evaluation of the nonlinear optical susceptibility in frequency-response based approaches. In particular the dipole matrix elements between the periodic part of the Bloch functions are ill-defined when using the standard definition of the position operator. In that case, it is possible to obtain correct expressions for the dipole matrix elements from perturbation theory Aversa and Sipe (1995); Sipe and Ghahramani (1993); Luppi et al. (2010) at a given order in the external field. Instead, in the real-time approach one needs an expression valid at each order of the perturbation.

A correct definition of the polarization operator in systems with PBC has been introduced by means of the geometric Berry phase in the Modern theory of polarization.Resta (1994) To our knowledge real-time schemes for calculating the electron-field coupling consistently with PBC have been proposed in Refs. Springborg and Kirtman, 2008; Virk and Sipe, 2007; Souza et al., 2004. In those works the dipole matrix elements are evaluated numerically from the derivative in the crystal-momentum () space. The latter cannot be carried out trivially because of the freedom in the gauge of the periodic part of the Bloch functions. In fact, the gauge freedom leads to spurious phase differences in the Bloch functions at two neighbouring points and ultimately to spurious contributions to the numerical derivative. Then, basically the three schemes Springborg and Kirtman (2008); Virk and Sipe (2007); Souza et al. (2004) differ in how the gauge is fixed to eliminate the spurious phase.

This work presents a real-time ab-initio approach to nonlinear optical properties for extended systems with PBC in which the nonlinear optical susceptibility are obtained through Eq. (1). To derive the EOMs we follow the scheme of Souza et al.Souza et al. (2004) based on the generalization of Berry phase to the dynamical polarization (Sec. II.1). Originally applied to a simple tight-binding Hamiltonian, this approach is valid for any single-particle Hamiltonian and, as we discuss in Sec. II.2, it can be applied in an ab-initio context with inclusion of the relevant many-body effects. After detailing how nonlinear optical susceptibility is extracted from the dynamical polarization (Sec. III), we show results for the second- and third-harmonic generation (THG) in semiconductors (Sec. IV) and successfully validate them against existing results from the literature obtained by response theory in frequency domain.

Ii Theoretical background

We consider a system of electrons in a crystalline solid of volume (where is the number of the equivalent cells and the cell volume) coupled with a time-dependent electric field


where is the zero-field Hamiltonian, and describes the coupling with the electric field. Here, we consider a generic single-particle Hamiltonian . In Sec. II.2 we specify the form of and show how many-body effects are included by means of effective single-particle operators. Of course, the choice of a single-particle Hamiltonian prevents applications to systems with strong static correlation such as Mott insulators or frustrated magnetic materials. We assume the ground state of to be non-degenerate and a spin-singlet so that the ground-state wavefunction can be expressed as a single Slater determinant. We also assume, as usual in treating cell-periodic systems, Born-von Kármán PBC and define a regular grid of -points in the Brillouin zone. With such assumptions, the single-particle solutions of are Bloch-functions.

Regarding the electron-field coupling we assume classic fields and use the dipole approximation, ( is the electronic charge). However, because of the PBC the position operator is ill-defined. In order to obtain a form for the field coupling operator compatible with Born-von Kármán PBC, in this paper we use the Berry phase formulation of the position operator and consequently the polarization. As proved in Ref. Souza et al., 2004, in this formulation the solutions of are also in a Bloch function form: , with being the periodic part and being the band index. Notice that, even in the Berry phase formulation, for very strong fields and with the number of -points that goes to infinity the Hamiltonian Eq. 2 is unbounded from below due to the Zener tunnelling.Springborg and Kirtman (2008) Nevertheless the strength of the fields used in non-linear optics is well below this limit.Springborg and Kirtman (2008); Souza et al. (2004)
In Sec. II.1 we detail how, by starting from the Berry phase formulation of polarization, we obtain the EOMs in presence of an external electric field within PBC.

ii.1 Treatment of the field coupling term

ii.1.1 Berry phase polarization

Developed in the mid-90s the Modern Theory of Polarization Resta (1994) provides a correct definition for the macroscopic bulk polarization, not limited to the perturbative regime, in terms of the many-body geometric phase


In Eq. (3) is the macroscopic polarization along the primitive lattice vector , , with the primitive reciprocal lattice vector such that , and the number of -points along , corresponding to the number of equivalent cells in that direction. Note that in this formulation the polarization operator is a genuine many-body operator that cannot be split as a sum of single-particle operators.

By using the assumption that the wave-function can be written as a single Slater determinant, the expectation value of the many-body geometric phase in Eq. (3) can be seen as the overlap between two single Slater determinants. The latter is equal to the determinant of the overlap matrix built out of , the occupied Bloch functions


Then we can rewrite Eq. (3) as


where is the spin degeneracy, equal to since we consider here only spin-unpolarized systems, and is the number of -points in the plane perpendicular to reciprocal lattice vector , with .

The overlap has dimensions , where is the number of doubly occupied bands. However, from the properties of the Bloch functions and by imposing they satisfy the so-called “periodic gauge” , it follows that the integrals in Eq. (4) are different from zero only if . Therefore the determinant of reduces to the product of determinants of overlaps built out of , the periodic part of the occupied Bloch functions:


This leads to the formula by which we compute the polarization of the system


Using matrix properties, Petersen and Pedersen (2012) the logarithm of the matrix determinant can be rewritten as the trace of matrix logarithm, and so Eq. (7) can be transformed as


more suitable to derive the EOMs. By taking the thermodynamic limit ( and ) of the latter expression one arrives at the King-Smith and Vanderbilt formula for polarization. King-Smith and Vanderbilt (1993) Since in a numerical implementation we deal with a finite number of and finite , we stick here to Eq. (8) with to derive the EOMs.

ii.1.2 Equations-of-motion

Following Ref. Souza et al., 2004 we start from the Lagrangian of the system in presence of an external electric field :


where is the energy functional corresponding to the zero-field Hamiltonian:


with . Notice that does not connect wave-functions with different vectors. To simplify the notation we do not explicit the time dependence of the , but they should be considered time-dependent in the rest of the paper.

We derive the dynamical equations from the Euler-Lagrange equations


To obtain the functional derivative of the polarization expression in Eq. (8) we use that Souza et al. (2004); Nunes and Gonze (2001)


and that exchanging arguments () in [Eq. (6)] brings a minus sign in Eq. (8). This leads to (see Ref. Souza et al., 2004 for details):


where , and from which we can define the field coupling operator


Notice that the field coupling operator in Eq. (16) is nonhermitian. In order to have well defined Hermitian operators in the EOMs we replace with . This is possible because at any time .Souza et al. (2004) Finally, by using Eqs. (14)-(16) in Eq. (12) and the Hermitian field coupling operator we obtain the EOMs:


Note that Eq. (16) contains


that has the form of the two-points central finite difference approximation of , but for the fact that are used instead of . As explained in Ref. Souza et al., 2004, the are built from the [Eq. (15)] in such a way that they transform as under a unitary transformation .

In fact, there is a gauge freedom in the definition of , that is , and since the Hamiltonian is diagonalized independently at each , the gauge is fixed independently and randomly at each . Then, standard (numerical) differentiation will be affected by the different gauge choices at two neighbouring -points. Instead the (numerical) derivative in Eq. (18) is gauge-invariant, or more specifically is performed in a locally flat coordinate system with respect to . In fact, in the thermodynamical/continuum limit, Eq. (18) corresponds to the covariant derivative. The problem of differentiating with respect to has been addressed also in Refs. Virk and Sipe, 2007; Springborg and Kirtman, 2008; Andrew and Tan, 1998 that use alternative approaches to ensure the gauge-invariance. In the here-discussed approach the definition of a numerical covariant differentiation originates directly from the definition of the polarization as a Berry phase.

ii.2 Treatment of electron correlation

Correlation effects play a crucial role in both linearOnida et al. (2002) and non linearLuppi et al. (2010); Cabellos et al. (2009) response of solids. Since we assumed that in Eq. (3) can be written as a single Slater determinant, effects beyond the IPA can be introduced in through an effective time-independent one-particle operator that can be either spatially local as in time-dependent density functional theory, or spatially nonlocal as in time-dependent Hartree-Fock.

However, both time-dependent density functional theory and time-dependent Hartree-Fock are not suitable approaches to optical properties of semiconductors: the former, within standard approximations for the exchange-correlation approximations, underestimates the optical gap and misses the excitonic resonances; the latter largely overestimates the band-gap and excitonic effects.

In the framework of Green’s function theory a very successful way to deal with electron-electron interaction in semiconductors is the combination of the approximation for the quasi-particle band structure Strinati et al. (1982) with the Bethe-Salpeter equation in static ladder approximation for the response function. Strinati (1988)

We recently extended this approach to the real-time domain Attaccalite et al. (2011) by mean of non-equilibrium Green’s function theory. In practice, the latter approach corresponds to a time-dependent static screened Hartree-Fock operator that satisfies the above-mentioned restrictions on the choice of and thus can be used within the here proposed framework. In what follows, we reformulate the approach in Ref. Attaccalite et al., 2011 as time-dependent Schrödinger like equations [Eq. (17)].

As starting point we choose the Kohn-Sham Harmiltonian at fixed density as a system of indipendent particles, Kohn and Sham (1965)


where is the electron-ion interaction, the Hartree potential and the exchange-correlation potential. The advantage of such a choice is that the Kohn-Sham system is the independent-particle system that reproduces the electronic density of the unperturbed many-body interacting system , thus by virtue of the Hohenberg-Kohn theorem Hohenberg and Kohn (1964) the ground-state properties of the system. Furthermore, no material dependent parameters need to be input, but for the atomic structure and composition.

As first step beyond the IPA, we introduce the corrections to the independent-particle energy levels by the electron-electron interaction through a (state-dependent) scissor operator


The latter can be calculated ab-initio e.g., via the approach , or can be determined empirically from the experimental band gap . We refer to this approximation as the independent quasiparticle approximation (QPA):


Notice that in our approach the inclusion of a non-local operator in the Hamiltonian does not present more difficulties than a local one, while in the response theory in frequency domain this is not a trivial task.Luppi et al. (2010) As a second step we consider the effects originating from the response of the effective potential to density fluctuations. By considering the change of the Hartree plus the exchange-correlation potential in Eq. 19 we will obtain the TD-DFT response. Here we include just “classic electrostatic” effects via the Hartree part. We refer to this level of approximation as the time-dependent Hartree (TDH)


In the linear response limit the TDH is usually referred as Random-Phase approximation and is responsible of the so-called crystal local field effects.Adler (1962)

Beyond the TDH approximation one has the TD-Hartree-Fock that includes the response of the exchange term to fluctuations of the density matrix . As discussed above this level of approximation is insufficient for optical properties of semiconductors, normally worsening over TDH results. The next step is thus to consider a screened exchange term in which the relevant electron correlation is introduced as a static screening term. Strinati (1988) The latter is calculated for the unperturbed KS system and is fixed to its initial value. We refer to this level of approximation as TD screened exchange or TD screened Hartree-Fock (TD-SHF),


We want to emphasize again that within this approach many-body effects are easily implemented by adding terms to the unperturbed independent-particle Hamiltonian in the EOMs [Eq. (17)]. Limitations may arise because of the computational cost of calculating those addition terms. In the specific the large number of -points needed to converge the SHG and THG spectra makes TD-SHF calculations impracticable. However, much less -points are needed for converging the screened-exchange self-energy itself and currently we are investigating how to exploit this property and devise “double grid” strategy similar to the one proposed in Ref. Kammerlander et al., 2012. In this work effects beyond IPA are limited to the QPA and TDH.

Finally, when the wave-function cannot be approximated anymore with a single Slater determinant (as in strong-correlated systems) the evaluation of the polarization operator [Eq. 3] becomes quite cumbersome.Stella et al. (2011) Also we are not aware of any successful attempt to combine Berry’s phase polarization with Green’s function theory or density matrix kinetic equations (including for example scattering terms), even if some appealing approaches have been proposed.Massidda et al. (1995); Chen and Lee (2011)

Iii Computational scheme and numerical parameters

Figure 1: [Color online] Proposed real-time ab-initio scheme to compute SHG and THG spectra in the energy range for extended systems with PBC: Results from KS-DFT and are input to determine the zero-field Hamiltonian. The EOMs [Eq. (17)] are then integrated to obtain the overlaps from which the polarization is computed as in Eq. (7). In the post-processing step the nonlinear susceptibilities are obtained by inversion of the Fourier matrix [Eq. (28)], see Sec. III for details.

Figure 1 illustrates the computational scheme we use to calculate the SHG and THG spectra. It consists in KS-DFT and calculations to determine the density and the KS eigenvalues, the quasiparticle corrections and eigenfunctions entering the zero-field Hamiltonian; the integration of the equation of motions [Eq. (17)] with a monochromatic electric field to obtain the from Eq. (7) and the post-processing of to extract the nonlinear susceptibilities. The latter two steps are repeated varying the laser frequency within the energy range for which we calculate the spectra.

The scheme in Fig. 1 has been implemented in the development version of the Yambo code. Marini et al. (2009) Kohn-Sham calculations have been performed using the Abinit code, Gonze et al. (2002) and the relevant numerical parameters are summarized in Table 1. All the operators appearing in the EOMs[Eqs. (17),(22),(21)] have been expanded in the Kohn-Sham basis set and the number of bands employed in the expansion is reported in Table 1.

Rigorously to have a fully ab-initio scheme, the scissor operator has to be calculated using e.g., . Here we use an empirical values for the scissor operator (reported in Table. 1) since the scope is to validate the computational scheme, and to facilitate the comparison with other works in the literature.

System PP (Ha) -grid (Å) Bands (eV)
SiC Si: 30 8/16 4.36 1-8 0.8
AlAs Al: 20 8/18 5.66 2-10 0.9
CdTe Cd: 40 8/14 6.48 7-13 1.0
Si Si: 14 8/24 5.39 1-7 0.6
Table 1: Parameters used in the Kohn-Sham density-functional theory and in the real-time simulations: the pseudopotential components for each atom (PP); energy cutoff for the planewaves (); number of points of the Monkhorst-Pack grid in each of the three dimensions for calculating the density and in the RT-simulations (-grid); the lattice parameter (); the range of bands for which the single-particle wave-function is evolved during the RT-simulation; the QP scissor () used within the QPA.

The EOMs [Eq. (17)] have been integrated using the following algorithm Koonin and C. (2008)


valid for both Hermitian and non-Hermitian Hamiltonians, and strictly unitary for any value of the time-step in the Hermitian case. In all real-time simulations we used a time-step of 0.01 fs. The number of states (number of bands number of points) evolved during the simulations is reported in Table 1.

Figure 2: [Color online] Pictorial representation of the signal analysis in the post-processing step. The signal (red line) can be divided into two regions: an initial convergence region (up to ) in which the eigenfrequencies of the systems are “filtered out” by dephasing and a second region where Eq. (26) holds. In this second region the signal is sampled within a period to extract the coefficients of Eq. 28. Note that is not a realistic one: for illustration purposes we enhanced the second-harmonic signal that otherwise would not be visible on this scale.

In our simulations we switch on the monochromatic field at . This sudden switch excites the eigenfrequencies of the system introducing spurious contributions to the non-linear response. We thus add an imaginary term into the Hamiltonian to simulate a finite dephasing:


where are the valence bands of the unperturbed system and is the dephasing rate. Then we run the simulations for a time much larger than and sample close to the end of the simulation, see Figure 2. Since determines also the spectral broadening, we cannot choose it arbitrary small. For example in the present calculations we have chosen of 6 fs that corresponds to a broadening of approximately 0.2 eV (comparable with the experimental one) and thus we run the simulations for 50-55 fs.
Once all the eigenfrequencies of the system are filtered out, the remaining polarization is a periodic function of period , where is the frequency of the external perturbation and can be expanded in a Fourier series


with , and complex coefficients:


To obtain the optical susceptibilities of order at frequency one needs to calculate the of Eq. (26), proportional to by the -th power of the . However, the expression in Eq. (27) is not the most computationally convenient since one needs a very short time step—significantly shorter than the one needed to integrate the EOMs—to perform the integration with sufficient accuracy. As an alternative we use directly Eq. (26): we truncate the Fourier series to an order larger than the one of the response function we are interested in. We sample values within a period , as illustrated in Figure 2. Then Eq. (26) reads as a system of linear equations


from which the component of in the direction is found by inversion of the Fourier matrix . We found that the second harmonic generation converges with S equal to 4 while the third harmonic requires S equal to 6. Finally we noticed that averaging averaging the results on more periods can slightly reduce the numerical error in the signal analysis.
Alternatively one can opt for a slow switch on of the electric field as in Takimoto et al.,Takimoto et al. (2007) so that no eigenfrequencies of the system are excited, and avoid to introduce imaginary terms in the Hamiltonian. We found, however, that the latter approach also requires long simulations, and on the other hand, it is less straightforward to extract the .

Iv Results

The main objective of this section is to validate the computational approach described in Secs. II and III against results in the literature for SHG obtained by the response theory in frequency domain. In particular we chose to validate against results from Refs. Luppi et al., 2010; Hübener et al., 2010 on bulk SiC and AlAs in which the electronic structures is obtained—as in our case—from a pseudopotential plane-wave implementation of Kohn-Sham DFT with the local density approximation, which makes the comparison easier. In the following we considered the zinc-blende structure of SiC and AlAs for which the tensor has only one independent nonzero component, (or its equivalent by permutation).

Figures 3 and 4 show results for the magnitude of SHG in SiC at the IPA, QPA and TDH level of theory. At all levels of approximation we obtained an excellent agreement with the results in Ref. Luppi et al., 2010. The minor discrepancies between the curves are due to the different choice for the -grid used for integration in momentum space: we used a -centered uniform grid (for which we can implement the numerical derivative) whereas Ref. Luppi et al., 2010 used a shifted grid. Figures 5 and 6 show results for the magnitude of SHG in AlAs at the IPA, QPA and TDH level of theory. Also in this case results obtained from our real-time simulations agree very well with the reference results and again the small differences between the spectra can be ascribed mostly to the different grid for -integrations.

Figure 3: [Color online] Magnitude of for bulk SiC calculated within the IPA (black triangles) and QPA (red circles). Each point corresponds to a real-time simulation at the given laser frequency (see Sec. III). Comparison is made with results obtained ab-initio by direct evaluation of the in Ref. Luppi et al., 2010 in IPA (grey solid line) and QPA (brown dashed line).
Figure 4: [Color online] Magnitude of for bulk SiC calculated within the IPA (black triangles) and TDH (red circles). Each point corresponds to a real-time simulation at the given laser frequency (see Sec. III). Comparison is made with results obtained ab-initio by direct evaluation of the in Ref. Luppi et al., 2010 in IPA (grey solid line) and TDH (brown dashed line).

As side results we can also observe the effects of different levels of approximation for the Hamiltonian on the SHG spectrum. In order to interpret those spectra note that SHG resonances occur when either or equals the difference between two single-particle energies. Then one can distinguish two energy region: below the single-particle minimum direct gap where only resonances at can occur, and above where both or resonances can occur.

Regarding the quasiparticle corrections to the IPA energy levels by a scissor operator, below the minimum Kohn-Sham direct band gap the IPA spectrum is shifted by half of the value of the scissor shift (0.4 eV for SiC and 0.45 eV for AlAs) and the spectral intensity reduced by a factor 1.18 (SiC) and 1.25 (AlAs). Above the minimum Kohn-Sham direct band gap instead the QPA spectrum cannot be simply obtained by shifting and renormalizing the IPA one because of the occurrence of resonances at , that are shifted and renormalized differently.

Regarding the crystal local field, their global effect is to reduce the intensity with respect to the IPA. For SiC, the intensity is reduced by about 15% below the gap, while above the band gap TDH and IPA have similar intensities. For AlAs we observe a reduction of about 30% in intensity for the whole range of considered frequencies, but for frequencies larger than 4 eV (that is where the resonances with the main optical transition occur) for which again the TDH and IPA have similar intensities.

Figure 5: [Color online] Magnitude of for bulk AlAs calculated within the IPA (black triangles) and QPA (red circles). Each point corresponds to a real-time simulation at the given laser frequency (see Sec. III). Comparison is made with results obtained ab-initio by direct evaluation of the in Refs. Luppi et al., 2010; Hübener et al., 2010 in IPA (grey solid and dot-dashed line) and QPA (brown dashed and dotted line).
Figure 6: [Color online] Magnitude of for bulk AlAs calculated within the IPA (black triangles) and TDH (red circles). Each point corresponds to a real-time simulation at the given laser frequency (see Sec. III). Comparison is made with results obtained ab-initio by direct evaluation of the in Ref. Luppi et al., 2010 in IPA (grey solid line) and TDH (brown dashed line).

We also computed the SHG of bulk CdTe (zincblende structure) in Fig. 7 and we compared with theoretical results Ghahramani et al. (1991) obtained by a minimal-basis semi-ab-initio approach (linear combination of Gaussian orbitals in conjunction with an Slater potential where is tuned to fit the gap) and with experimental results. Shoji et al. (1997); Jang et al. (2013) At the QPA level (scissor operator of 1.0 eV) the calculated spectrum differs noticeably from the theoretical results in Ref. Ghahramani et al., 1991, and largely overestimates the experimental intensities. Interestingly crystal local fields are strong, reducing the intensity of about 50% with respect to the IPA. The intensity of the TDH spectrum is then consistent with the experimental measurements.

Figure 7: [Color online] Magnitude of for bulk CdTe calculated within the QPA (black triangles) and TDH (red circles). Each point corresponds to a real-time simulation at the given laser frequency (see Sec. III). Comparison is made with results obtained ab-initio by direct evaluation of the in Ref. Ghahramani et al., 1991 in IPA (grey solid line) and available experimental results in Refs. Shoji et al., 1997; Jang et al., 2013 (blue stars and diamonds).

Finally, our approach can also compute third order susceptibilities as shown in Figs. 8 and 9 for the THG of Si. Within the dipole approximation bulk Si does not have SHG since is centrosymmetric. The first nonlinear effects we can extract from our simulation are at the third order. The THG for Si (diamond structure) has two independent components, and . In the expression for the TH polarization along the direction ,

is the isotropic contribution, while the anisotropic contribution. Combination of measurements of THG at different field polarizations provide the anisotropy and the phase of . Fig. 8 shows our results within the QPA and RPA (scissor operator of 0.6 eV) for and compared with theoretical calculations at the tight-binding level with either semi-ab-initio or empirical parameters. Apparently the empirical tight-binding results are the closer to our QPA spectra; however, the semi-ab-initio tight-binding spectra show the same peak structure and similar ratio between and intensities. Both spectra at the RPA level are very similar to the QPA ones: the isotropic contribution is practically identical, while slightly more pronounced differences can be observed for the anisotropic contribution. Fig. 9 shows the anisotropy and the phase compared with both other theoretical results and experiments. For energies below 1 eV our QPA spectra is in good agreement with results obtained from semi-ab-initio tight-binding and with the experimental measurement. For higher energies our spectra are less structured with respect both the semi-ab-initio tight-binding and the experiment. In particular the peak at 1-1.1 eV is missing. On the other hand, the intensities of the spectra are in better agreement with the experiment than the previous theoretical results. It is interesting to observe how small changes in the spectrum of the anisotropic contribution due to crystal local field effects induce quite important changes in the phase and anisotropy spectra. In particular, as previously observed they increase the anisotropy. Then, it seems important to include crystal local field effects even if they are weak.

Figure 8: [Color online] Magnitude of the two independent components of the THG in bulk Si: (left panel) and (right panel) calculated within the QPA (black triangles). Each point corresponds to a real-time simulation at the given laser frequency (see Sec. III). Comparison is made with results obtained by direct evaluation of the in Ref. Moss et al., 1990 from empirical (grey solid line) and semi-ab-initio tight-binding (brown dashed line).
Figure 9: [Color online] Dispersion in the magnitude of the anisotropy (top panel) and in the relative phase of (bottom panel) calculated within the QPA (black triangles). Each point corresponds to a real-time simulation at the given laser frequency (see Sec. III). Comparison is made with experimental results from Ref. Moss et al., 1989 (blue stars) and results obtained by direct evaluation of the in Ref. Moss et al., 1990 from empirical (grey solid line) and semi-ab-initio tight-binding (brown dashed line).

V Conclusions

We presented an ab-initio real-time approach to calculate nonlinear optical properties of extended systems. The key strengths of the proposed approach are first, the correct treatment of the coupling between electrons and external field and second the possibility of easily include effects beyond the IPA.

Regarding the treatment of the electron-field coupling, following the work of Souza et al.Souza et al. (2004), we started from the Berry-phase formulation for the dynamical polarization—a definition consistent with the PBC—to derive a covariant numerical expression for the dipole operator in the EOMs.

Note that we worked in the length-gauge even if the velocity gauge may appear a more natural choice. In fact, as opposed to the position operator the velocity operator is consistent with the PBC. However, in the velocity gauge even if the position operator disappears from the Hamiltonian, it reappears in the phase factor for the wavefunction, Lamb et al. (1987) so that the problem of re-defining the position operator remains. Furthermore, the velocity gauge is plagued by unphysical numerical divergences for the response to low frequencies. Aversa and Sipe (1995)

Regarding effects beyond the independent-particle approximation, they are included by simply adding the corresponding operator to the single-particle Hamiltonian. This is an easy task when compared with deriving the corresponding expressions for the nonlinear optical susceptibility. Cabellos et al. (2009); Virk and Sipe (2009) As an example, in the present work we have included quasiparticle corrections to the band-gap by adding to the Hamiltonian a scissor operator and crystal local-field effects by adding the time-evolution of the Hartree potential. In principle, one can add as well excitonic effects by adding the time-evolution of the screened exchange self-energy as in the scheme proposed in Ref. Attaccalite et al., 2011; or perform a real-time TD-DFT calculations by adding the time-evolution of the exchange-correlation potential. Being the focus of this work the validation of the proposed approach for calculating nonlinear properties, we leave the inclusion of these correlation effects for future work.

We have proved the validity of our approach for the SHG in bulk SiC and AlAs by showing an excellent agreement between our results, obtained from real-time simulations, and results in the literature obtained from direct evaluation of the second order susceptibility in frequency-domain. For CdTe we have computed the SHG and shown that local field effects are important to reproduce experimental measurements. Finally, our approach is not limited to the SHG and we have computed the phase and anisotropy of the THG in bulk Si and obtained results consistent with existing experimental measurements.

Vi Acknowledgements

The authors acknowledge Angel Rubio, Ivo Souza, Raffaele Resta, and Andrea Marini for useful discussion, and suggestions. M.G. acknowledges the Portuguese Foundation for Science and Technology for funding (PTDC/FIS/103587/2008) and support through the Ciência 2008 program. Computing time has been provided by the national GENGI-IDRIS supercomputing centers at Orsay under contract i2012096655 and by the HECToR Supercomputing Facility through EPSRC grant EP/K013459/1, allocated to the UKCP Consortium.


  • Aryasetiawan and Gunnarsson (1998) F. Aryasetiawan and O. Gunnarsson, Reports on Progress in Physics 61, 237 (1998).
  • Strinati (1988) G. Strinati, Rivista del Nuovo Cimento 11, 1 (1988).
  • Onida et al. (2002) G. Onida, L. Reining,  and A. Rubio, Rev. Mod. Phys. 74, 601 (2002).
  • Sipe and Ghahramani (1993) J. E. Sipe and E. Ghahramani, Phys. Rev. B 48, 11705 (1993).
  • Cabellos et al. (2009) J. L. Cabellos, B. S. Mendoza, M. A. Escobar, F. Nastos,  and J. E. Sipe, Phys. Rev. B 80, 155205 (2009).
  • Luppi et al. (2010) E. Luppi, H. Hübener,  and V. Véniard, Phys. Rev. B 82, 235201 (2010).
  • Cazzanelli et al. (2012) M. Cazzanelli, F. Bianco, E. Borga, G. Pucker, M. Ghulinyan, E. Degoli, E. Luppi, V. Véniard, S. Ossicini, D. Modotto, S. Wabnitz, R. Pierobon,  and L. Pavesi, Nat Mater 11, 148 (2012).
  • (8) S. Botti, F. Sottile, N. Vast, V. Olevano, L. Reining, H.-C. Weissker, A. Rubio, G. Onida, R. Del Sole,  and R. W. Godby, Phys. Rev. B 69, 155112.
  • Leitsmann et al. (2005) R. Leitsmann, W. G. Schmidt, P. H. Hahn,  and F. Bechstedt, Physical Review B 71, 195209 (2005).
  • Chang et al. (2001) E. K. Chang, E. L. Shirley,  and Z. H. Levine, Physical Review B 65, 035205 (2001).
  • Virk and Sipe (2009) K. S. Virk and J. E. Sipe, Phys. Rev. B 80, 165318 (2009).
  • Yabana and Bertsch (1996) K. Yabana and G. F. Bertsch, Phys. Rev. B 54, 4484 (1996).
  • Takimoto et al. (2007) Y. Takimoto, F. D. Vila,  and J. J. Rehr, The Journal of Chemical Physics 127, 154114 (2007).
  • Castro et al. (2004) A. Castro, M. A. L. Marques,  and A. Rubio, The Journal of Chemical Physics 121, 3425 (2004).
  • Meng and Kaxiras (2008) S. Meng and E. Kaxiras, The Journal of Chemical Physics 129, 054110 (2008).
  • Boyd (2008) R. W. Boyd, Nonlinear Optics (Academic Press, 2008).
  • Attaccalite et al. (2011) C. Attaccalite, M. Grüning,  and A. Marini, Phys. Rev. B 84, 245110 (2011).
  • Aversa and Sipe (1995) C. Aversa and J. E. Sipe, Phys. Rev. B 52, 14636 (1995).
  • Resta (1994) R. Resta, Rev. Mod. Phys. 66, 899 (1994).
  • Springborg and Kirtman (2008) M. Springborg and B. Kirtman, Phys. Rev. B 77, 045102 (2008).
  • Virk and Sipe (2007) K. S. Virk and J. E. Sipe, Phys. Rev. B 76, 035213 (2007).
  • Souza et al. (2004) I. Souza, J. Íñiguez,  and D. Vanderbilt, Phys. Rev. B 69, 085106 (2004).
  • Petersen and Pedersen (2012) K. B. Petersen and M. S. Pedersen, “The matrix cookbook,”  (2012), version 20121115.
  • King-Smith and Vanderbilt (1993) R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993).
  • Nunes and Gonze (2001) R. W. Nunes and X. Gonze, Phys. Rev. B 63, 155107 (2001).
  • Andrew and Tan (1998) A. L. Andrew and R. C. Tan, SIAM Journal on Matrix Analysis and Applications 20, 78 (1998).
  • Strinati et al. (1982) G. Strinati, H. J. Mattausch,  and W. Hanke, Phys. Rev. B 25, 2867 (1982).
  • Kohn and Sham (1965) W. Kohn and L. J. Sham, Phys. Rev. 140, A1133 (1965).
  • Hohenberg and Kohn (1964) P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964).
  • Adler (1962) S. L. Adler, Phys. Rev. 126, 413 (1962).
  • Kammerlander et al. (2012) D. Kammerlander, S. Botti, M. A. L. Marques, A. Marini,  and C. Attaccalite, Phys. Rev. B 86, 125203 (2012).
  • Stella et al. (2011) L. Stella, C. Attaccalite, S. Sorella,  and A. Rubio, Physical Review B 84, 245117 (2011).
  • Massidda et al. (1995) S. Massidda, R. Resta, M. Posternak,  and A. Baldereschi, Physical Review B 52, R16977 (1995).
  • Chen and Lee (2011) K.-T. Chen and P. A. Lee, Phys. Rev. B 84, 205137 (2011).
  • Marini et al. (2009) A. Marini, C. Hogan, M. Gruning,  and D. Varsano, Comp. Phys. Comm. 180, 1392 (2009).
  • Gonze et al. (2002) X. Gonze et al., Computational Materials Science 25, 478 (2002).
  • Koonin and C. (2008) S. E. Koonin and M. D. C., eds., Computational Physics: Fortran Version (Perseus Books, 2008).
  • Hübener et al. (2010) H. Hübener, E. Luppi,  and V. Véniard, Phys. Status Solidi B 427, 1984 (2010).
  • Ghahramani et al. (1991) E. Ghahramani, D. J. Moss,  and J. E. Sipe, Phys. Rev. B 43, 9700 (1991).
  • Shoji et al. (1997) I. Shoji, T. Kondo, A. Kitamoto, M. Shirane,  and R. Ito, J. Opt. Soc. Am. B 14, 2268 (1997).
  • Jang et al. (2013) J. I. Jang, S. Park, D. J. Clark, F. O. Saouma, D. Lombardo, C. M. Harrison,  and B. Shim, J. Opt. Soc. Am. B 30, 2292 (2013).
  • Moss et al. (1990) D. J. Moss, E. Ghahramani, J. E. Sipe,  and H. M. van Driel, Physical Review B 41, 1542 (1990).
  • Moss et al. (1989) D. Moss, H. M. van Driel,  and J. E. Sipe, Optics letters 14, 57 (1989).
  • Lamb et al. (1987) W. E. Lamb, R. R. Schlicher,  and M. O. Scully, Phys. Rev. A 36, 2763 (1987).
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description