Ab initio nonlinear optics in solids: linear electro-optic effect and electric-field induced second-harmonic generation

Second-harmonic generation (SHG), linear electro-optic effect (LEO) and electric-field induced second-harmonic generation (EFISH) are nonlinear optical processes with important applications in optoelectronics and photovoltaics. SHG and LEO are second-order nonlinear optical processes described by second-order susceptibility. Instead, EFISH is a third-order nonlinear optical process described by third-order susceptibility. LEO and EFISH are only observed in the presence of a static electric field. These nonlinear processes are very sensitive to the symmetry of the systems. In particular, LEO is usually observed through a change in the dielectric properties of the material while EFISH can be used to generate a “second harmonic” response in centrosymmetric material. In this work, we present a first-principle formalism to calculate second- and third-order susceptibility for LEO and EFISH. LEO is studied for GaAs semiconductor and compared with the dielectric properties of this material. We also present how it is possible for LEO to include the ionic contribution to the second-order macroscopic susceptibility. Concerning EFISH we present for the first time the theory we developed in the framework of TDDFT to calculate this nonlinear optical process. Our approach permits to obtain an expression for EFISH which does not contain the mathematical divergences in the frequency-dependent second-order susceptibility that caused until now many difficulties for numerical calculations.


Introduction
A deep understanding of the nonlinear optical properties of solids [1,2] provides an opportunity to search for Among all nonlinear phenomena existing in nature, a major role is played by second-harmonic generation (SHG), linear electro-optic effect (LEO) and electricfield induced second-harmonic generation (EFISH).
SHG is a process in which two photons of frequency ω are absorbed by the material and a photon of frequency 2ω is emitted. SHG is described by the second-order macroscopic susceptibility χ (2) (−2ω; ω, ω). LEO is a second-order nonlinear optical process as well, observed in the presence of a static electric field. LEO can be seen as a process in which two photons are absorbed, one at frequency ω and one at vanishing frequency and therefore it is described by the second-order macroscopic susceptibility χ (2) (−ω; ω, 0). For this reason, the electro-optic effect is usually observed by a change in the dielectric properties of the system, proportional to the static electric field.
SHG and LEO are only observed for systems which has no inversion symmetry, because χ (2) is equal to zero for centrosymmetric materials. As a consequence, they are highly sensitive to the symmetry of the system. SHG has been largely used as a probe to investigate the properties of bulk, surfaces, interfaces and complex systems [3][4][5][6][7][8].
LEO has many applications in the development of optoelectronic devices. In fact, the nonlinear photoinduced effects in the materials are related to the fast or ultrafast response of charge carriers within the materials and it has been shown that the electro-refractive effect is a promising route to realize efficient high speed optical modulators [9]. Moreover, a giant electro-optic effect has been observed in Ge/SiGe coupled quantum wells, therefore enhancing the performance of optical modulators [10]. Fast LEO can be disentangled by SHG in Si waveguides [11,12]. On the other hand, applying an asymmetric strain in crystals can induce noncentrosymmetry and it has been reported that strain in a Si photonic crystal waveguide leads to strong nonlinearities and enables electro-optic effects, showing that data processing could be potentially performed by allsilicon components [13,14]. Finally, because of its sensitivity to space symmetry, LEO can also be used as a sensitive, nondestructive and noninvasive, probe for studying many kinds of surfaces and interfaces in semiconductors [15,16].
The presence of a static electric field inside a material also enables second harmonic generation, through EFISH, a third order process, described by the thirdorder macroscopic susceptibility χ (3) (−2ω; ω, ω, 0). Contrary to SHG and LEO, EFISH can be used to generate a "second harmonic" response in centrosymmetric material. Moreover, despite the fact that EFISH is a third-order process, expected to be smaller than SHG, in materials such as surfaces and interfaces, EFISH can be competitive to SHG. In fact, SHG and EFISH processes can be present simultaneously and must be distinguished.
EFISH, as well as SHG and LEO, has many applications in the development of optoelectronic devices. Recently, it was also reported that ultrafast recombination together with carrier diffusion can be monitored by EFISH generated by space charge accumulation in the material. [4,17,18].
The theoretical description of SHG, LEO and EFISH requires the knowledge of the first-, second-and thirdorder susceptibility which is a formidable task for an ab initio theory [19][20][21]. Some of the authors developed an ab initio formalism in the framework of the time-dependent density functional theory (TDDFT) to calculate χ (2) in order to describe SHG. Within this approach, SHG is described with different levels of approximation for the electronic correlations : independent particle approximation (IPA), random-phase approximation (RPA) and excitons [21,22]. Moreover, the quasiparticle effects have been included through a scissor operator. Recently, this formalism was extended to describe LEO in IPA and including quasiparticle and excitons [23].
Concerning EFISH, most of the calculations have been performed using phenomenological models [24,25] in some cases integrated with group theory formalism in order to investigate the relation of EFISH with the symmetry of the crystal [26]. A frequencydependent formulation of the third-order susceptibility was presented [27]. However, this equation was proven to be difficult for computational calculations because of numerical divergences and it was not used to calculate frequency-dependent spectra. Another derivation of the EFISH frequency-dependent susceptibility from perturbation theory in independent particle approximation (IPA) was proposed [28]. This formalism was applied only in the energy range close to critical point energy gaps [28].
In this work, we present a first-principle formalism to calculate second-and third-order susceptibility for LEO and EFISH. In the case of LEO we also show how to include the ionic contribution to the second-order frequency-dependent susceptibility. Concerning EFISH we show how it is possible to remove the mathematical divergences in the frequency-dependent susceptibility that caused until now many difficulties for numerical calculations.
The paper is organised as follows: in Sect. 2 we present the theory for LEO and EFISH and in Sect. 3 we present results for LEO and EFISH and also SHG in the case of GaAs, Si and SiC semiconductors. Conclusions are in Sect. 4. The mathematical conventions used for the susceptibilities are given in Appendix 1 and details of the calculation for the ionic contribution to the LEO susceptibility are presented in Appendix 2.

Electronic contribution
The macroscopic polarisation at frequency ω of a material irradiated by a time-dependent electric-field E(ω) together with a static electric-field E E E is P(ω) = P (1) (ω) + P (2) (ω). The i-th component (x , y, z ) of the macroscopic first-order polarisation is related to the j-th component (x , y, z ) of the electric field as where the first-order susceptibility χ where the second-order susceptibility χ (2) ijk is a 27component tensor. Atomic units are used throughout unless otherwise stated.
The number of independent and non-zero components of the susceptibilities χ (1) ij and χ (2) ijk are entirely determined by the symmetry of the material studied. The mathematical conventions used to define the susceptibilities are in Appendix 1.
Considering the electric displacement D(ω) = E(ω)+ 4πP(ω) and replacing P(ω) by the first-and the secondorder macroscopic polarisation (Eqs. (1) and (2)), it is possible to write D(ω) =˜ (ω)E(ω) in terms of an effective dielectric matrix˜ (ω) defined as where ij (ω) is the linear dielectric matrix and depends on the second-order susceptibility in the presence of a static electric field E E E polarised along the k direction.
The quantity E k give rise to additional terms than those already included in the linear dielectric matrix ij (ω). This is easily observed considering systems with zinc-blende symmetry. In this case, the dielectric tensor matrix ij (ω) is diagonal and the presence of a static electric field E E E polarised along z direction gives rise to the off diagonal term Ez xy = 8πχ where xx = yy = zz . The effective dielectric matrix˜ (ω) is related to the LEO coefficients r ijk (ω) as The calculation of the LEO coefficients r ijk (ω) implies therefore the calculation of χ (2) ijk (−ω; ω, 0). As mentioned earlier, some of the authors developed an abinitio formalism in the framework of the TDDFT to calculate χ (2) ijk (−ω; ω, 0) starting from IPA and then including quasiparticle (scissor approximation) and excitonic effects. [23] Within this formalism [23] the χ is expressed in terms of the second-order response function which is where G is a vector of the reciprocal lattice, f n,k are Fermi occupation numbers (1 for occupied states and 0 for unoccupied states) labelled with the number of bands (n) and k-points (k) in the first Brillouin zone, and V is the volume of the cell. Moreover, in the above expression, the factor 2 accounts for the spin and η is a small vanishing quantity. The wave-vectors q 1 and q 2 are along the polarization of the incoming electric fields while q = q 1 + q 2 is in the polarization direction of the outgoing electric field. In Eq. (8) it also appears the quantityρ which is defined as the matrix elementρ nn k+q (q + G) = φ n,k+q |e i(q+G)r |φ n ,k in terms of Bloch wave functions φ n,k . In the case we considered the operator e −i(q+G)r in the calculation we indicated the matrix element asρ nn k+q (−(q + G)). In the low frequency range we consider the optical limit (q → 0) and we used the k.p perturbation theory to evaluate the matrix elements.
To go beyond IPA in order to include excitonic effects, we solved the TDDFT Dyson equation for χ (2) which is [21] 1 − χ (1) where we omitted the explicit dependence on the q and G-vectors for readability. [21] Here, f vxc is the sum of the bare-coulomb potential v and of the exchange-correlation kernel f xc . Besides f xc , another exchange-correlation kernel g xc appears in the secondorder TDDFT Dyson equation. In the framework of TDDFT, the kernels are unknown quantities and have to be approximated. Finally, χ (1) 0 and χ (2) 0 are the linear and second-order response functions in IPA while χ (1) is the linear response function calculated via the following linear-order Dyson equation [21]

Ionic contribution
The calculation of the LEO susceptibility and therefore of the LEO coefficients requires also the ionic contribution in addition to the electronic one.
The ionic contribution to the susceptibility is related to the ionic displacements induced by the static electric field, denoted as R(E E E) and depends on the variation of the dielectric tensor induced by these displacements. The sum of the electronic and ionic terms is referred to as the clamped value and the ratio of these two contributions is the Faust-Henry coefficient C F H = r i /r e . Taking into account the electronic and the ionic parts allows for a direct comparison of the theoretical calculation with the experimental measurements. While some experimental values of C F H have been published for typical semiconductors, such as GaAs or GaN [29][30][31], very few ab initio calculations are available. In a pioneering work, Veithen et al. [32] calculated ab initio both the electronic and the ionic static contributions to the susceptibility. Their work is based on the computation of the total derivative of the static dielectric tensor with respect to the static electric field.
However, the numerical evaluation of the frequencydependent Faust-Henry coefficient is still missing. We have extended this formulation to the dynamical case by considering the dielectric tensor as a function of the frequency ω, the amplitude of the static electric field E and the ionic displacements R(E E E). As the ionic displacements depend on the static electric field E E E, the total derivative of the dielectric tensor with respect to E contains two terms: the first one is obtained by considering the atoms at their equilibrium position R 0 while the second one depends on the ionic displacements where τ nα = R nα − R 0,nα is the displacement of atom n in the direction α. The first term in Eq. 11, obtained by considering the atoms at their equilibrium position R 0 , is the aforementioned electronic contribution while the second term is the ionic contribution and involves the derivative of the dielectric tensor with respect to the atomic displacements τ nα and the first-order electric field induced atomic displacement Following [32], τ nα is expanded in the basis of the zone-center phonon-mode eigendisplacements and we get where the summation runs over the phonon modes m and p mk is the polarity mode. For a detailed derivation see Appendix 2.
We finally obtain that the clamped LEO susceptibility with electronic and ionic terms is

EFISH
The macroscopic polarization of a material irradiated by a time-dependent electric-field E(ω) together with a static electric-field E E E at the frequency 2ω is P(2ω) = P (2) (2ω) + P (3) (2ω). The i-th component (x , y, x ) of the macroscopic second-order polarization is which contains the second-order susceptibility tensor χ (2) ijk and contributes to SHG. The EFISH is instead related to the macroscopic third-order polarization as : where the third-order susceptibility χ (3) ijkl is a 81component tensor. As for SHG, the number of independent and non-zero components are entirely determined by the symmetry of the material studied. In the case of EFISH, the indices j and k are interchangeable χ When the material is irradiated by a time-dependent electric-field E(ω) together with a static electric-field E E E, we obtain an effective susceptibility defined as For simplicity, we consider only systems with cubic symmetry and where the static electric field is applied in the z -direction, i.e. E z . Using the Voigt notation, it is possible to write the effective susceptibility which is a 3 × 3 × 3 tensor in a more convenient way as a 3 × 9 tensor which reads asχ Here the 9 columns correspond to components xx , yy, zz , yz , zy, zx , xz , xy, yx , respectively. Therefore, the effective susceptibility tensor contains the component χ (2) xyz which is also the only second-order component different from zero for zinc-blende symmetry together with the following three independent components induced by the static field In practice, to calculateχ (2) , it is necessary to know χ (2) and χ (3) . We have extended the methodology developed for the calculation of the χ (2) in the framework of TDDFT to the calculation of the χ (3) in the IPA and for the specific case of the EFISH process [33]. In order to calculate the χ (3) for EFISH we used the general expression for the third order response function valid for any frequency of the incoming fields and written in terms of the Bloch wave functions φ n,k , which reads as q1, q2, q3, ω1, ω2, ω3) = 2 V k n,n ,n ,n ṽ nn k (q)ṽ n n k (q1)ṽ n n k (q2)ṽ n nk (q3) ω ω1 ω2 ω3(E nn ,k + ω1 + ω2 + ω3 + 3iη) ω3)), (20) with ω = ω 1 + ω 2 + ω 3 , η being a vanishingly small quantity. In Eq. (20) the limit q → 0 is already understood [21] and f nm,k = f n,k − f m,k is the difference between the occupation numbers f n,k and f m,k , E nm,k = E n,k −E m,k is the difference between the band energies E n,k and E m,k andṽ nn k (q) = φ n,k |qv|φ n ,k is the matrix element of the velocity operator.
In the case of EFISH, the frequencies in Eq. (20) are ω 1 = ω 2 = ω and ω 3 = 0, which leads to an apparent non-physical divergence. The simplest way to remove the divergence is to evaluate the third-order response function for EFISH in the following way as lim δω→0 + = lim δω→0 − . The calculation of the limit δω → 0 is done analytically using a first order expansion of the function χ (3) 0 in terms of δω. The physical interpretation of the EFISH process appears then clearly in Fig. 1. After some tedious calculation, the final expression for the third-order response function is obtained and is given in [33]. This result permits to calculate the third-order susceptibility in IPA as , q 1 , q 2 , q 3 , ω, ω, 0). (22)

Results
We calculated SHG, LEO and EFISH for Si, SiC and GaAs bulk semiconductors. We first calculated the electronic structure of the materials in their ground-state within Density-Functional Theory (DFT) in the localdensity approximation (LDA), using norm-conserving and plane-wave basis set with the ABINIT code [34][35][36]. In the case of Ga, we explicitly included the 3d semicore states as valence states, giving the valence configuration 3d 10 4s 2 4p 1 , while only valence states are considered for the other materials. We used the experimental lattice constants 5.41Å for Si, 4.36Å for SiC and 5.63Å for GaAs and the cutoff energy was 20 Ha for Si, 40 Ha for SiC and 50 Ha for GaAs.
Then, we calculated the electronic contribution of the optical susceptibilities from the 2light code. In this code the third-and the second-order nonlinear optical response, in the framework of TDDFT, is implemented (see Eqs. (8) and (20)) using plane-wave basis set. [21,33] We used 27,000 shifted k points for Si, SiC and GaAs. The number of unoccupied bands was 16 for Si and SiC and 36 for GaAs, 965 G vectors for the wavefunctions for Si and SiC and 5005 for GaAs. Crystal local-field effects are not taken into account. In the case of LEO we used SO = 0.85 eV for GaAs.
For the ionic contribution, the evaluation of the frequency of the phonon mode ω m and the polarity mode p mk are obtained using the ABINIT code [34][35][36], while the derivative of ij with respect to the atomic displacement is evaluated through the calculation of ij (R, ω), with the linear TDDFT code, DP [37], for fixed ionic positions R, obtained by moving the ions along the phonon eigenmodes.

LEO
GaAs has a zinc-blende symmetry and therefore the dielectric tensor (see Eq. (4)) is diagonal, with xx = yy = zz and with χ (2) xyz the only nonvanishing secondorder component. As explained in Sect. 2 by applying an external static electric-field along the z-Cartesian axis, the dielectric tensor acquires an off-diagonal contribution which is Ez = 8πχ In Fig. 2 we show the imaginary part of the linear dielectric tensor and the Ez for an electric field of 4·10 5 V cm −1 , which is chosen to be weak enough in order not to destroy the material (at the limit of the electrical breakdown). The two components are displayed on different scales since the field-induced one is much smaller and would be indistinguishable otherwise. Note that for off-diagonal components the imaginary part does not need to be positive.
In Fig. 3 several contributions for the LEO second order susceptibility are presented. The electronic (green curve) and ionic (red curve) contributions are displayed as a function of the frequency of the electric field. For clarity, only the absolute value is presented, but one should remember that they are both complex quantities. It is important to stress that in the low energy part of the spectrum, where the imaginary part is small, the electronic and ionic contributions have opposite sign and therefore partially cancel each other. This was already known from the experimental value of the Faust-Henry coefficient C F H = −0.51, measured at ω = 1 eV [29]. The absolute value of the clamped susceptibility (electronic + ionic) is compared to the experimental values (blue points) [38], showing a good agreement. The frequency-dependent Faust-Henry coefficient is displayed in the inset. One can conclude from these results that the frequency dependence for all these quantities (ionic contribution and Faust-Henry coefficients) is far from being negligible. One also notes the good agreement between the theoretical and experimental value for C FH at ω = 1 eV.

EFISH
Two components for χ (3) (−2ω; ω, ω, 0) corresponding to the EFISH process for Si and SiC are shown in Figs. 4 and 5. We note that the intensity of the susceptibility is much higher for Si than for SiC, note the change of scale when comparing Figs. 4 and 5. The common feature of the two compounds is the fact that one finds a factor close to 2 between χ zzzz and χ zxxz . The effective second-order susceptibility is χ (2) zzz (−2ω; ω, ω) = 3χ (3) zzzz (−2ω; ω, ω, 0)E z . To compare the intensity of SHG and EFISH, we show on the same plot, Fig. 6 the intensity of the induced second-order componentχ (2) zzz for both Si and 3C-SiC, which displays a factor 10 difference between the two, clearly showing that the EFISH response in silicon is much larger than in SiC. The strength of the static field is 5 · 10 5 V cm −1 and we see that we can actually compare the SHG and EFISH components on the same scale, unlike what was done for LEO, showing that the field-induced component, while smaller than the χ (2) , is far from Fig. 6 Comparison for 3C-SiC between the two χ (2) components xyz coming from SHG (black curve) and zzz coming from EFISH (red curve) with a static field of Ez = 5 · 10 5 V cm −1 . The zzz -component for Si is also shown for the same static field amplitude being negligible in SiC. Note that there is no χ (2) xyz contribution for Si due to centrosymmetry.

Conclusion
In this work we presented the ab initio formalism we developed for the calculation of the second-order macroscopic susceptibility for LEO and of the third-order macroscopic susceptibility for the EFISH. In the case of LEO we have included in our theoretical approach the quasiparticle effect at the level of the scissors correction and the electronic and ionic contributions to the optical response. We have applied our method to the calculation of LEO coefficients for the semiconductor GaAs and we have compared our results with the experimental data presented in literature, finding a good agreement. We can conclude from this comparison that it is important to account for the frequency-dependent ionic contribution. Concerning EFISH we presented our formalism where the mathematical divergences in the third-order susceptibility have been removed. We have presented numerical values for the EFISH susceptibility for Si and SiC, showing in particular the high potential of Si as a non-linear medium.
The ab initio calculation of LEO and EFISH shows that an high accuracy can be reached both for the electronic and the ionic parts of the nonlinear coefficients. Moreover, the analysis and comparison of the SHG, the LEO and the EFISH together with the dielectric properties for different systems open the way for an accurate investigation of complex materials with technological interest. Data availibility statement The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Appendix 1
Writing a time-dependent electric field as the time-dependent second-order polarization induced by E(t) is SHG E 2 0 e −2iωt + e 2iωt + 2χ (2) OR E 2 0 (24) where the second-order susceptibility χ (2) relates the polarization to the electric field. The first term corresponds to the second harmonic generation (SHG) and the second term to the optical rectification (OR).One can define a frequency-dependent SHG an OR polarization as P (2) SHG (2ω) = χ (2) SHG E 2 0 and P (2) OR = 2χ (2) OR E 2 0 . In the general case, where the electric field contains two different frequencies, E(t) = E 0,1 e −iω1t + e iω1t + E 0,2 e −iω2t + e iω2t (25) the time-dependent second-order polarization becomes where the last two terms correspond respectively to the sum frequency generation (SFG) and the difference frequency generation (DFG). For the linear electro-optic effect (LEO), the total field is E(t) = E 0 e −iω1t + e iω1t + E dc (27) where E dc = lim ω→0 E 2 (t) = 2E 0,2 . The polarization corresponding to the linear electro-optic effect is then P LEO (t) = 2χ One can generalize these conventions to the third order case and we get for the EFISH process P (3) EFISH (t) = 3χ 3 EFISH E dc E 2 0 e −2iωt + e 2iωt (29)