You are an expert at summarizing long articles. Proceed to summarize the following text: 
nowadays , density functional theory  @xcite ( dft ) is considered as a standard method in condensed matter physics , to study electronic , structural and macroscopic properties of solids from first principles . combined with adiabatic perturbation theory , it allows _ a priori _ the computation of derivatives of the energy and related thermodynamic potentials up to any order . at the second order , this approach has been applied to compute linear response functions such as phonon frequencies or born effective charges with an accuracy of a few percents . the third - order derivatives are related to non - linear properties such as phonon lifetimes , raman tensors or non - linear optical susceptibilities . the linear - response formalism has been implemented in various first - principles codes and is routinely applied to an increasing number of systems [ see for example ref . [ ] and references therein ] . by contrast , the non - linear response formalism has been mostly restricted to quantum chemistry problems . although the hyperpolarizabilities of a huge number of molecules have been computed , taking into account both electronic and vibrational ( ionic ) contributions  @xcite , applications in condensed matter physics have focused on rather simple cases  @xcite .    here , we present a methodology for the computation of nonlinear response functions and related physical quantities of periodic solids from density functional perturbation theory ( dfpt ) . we focus on perturbations characterized by a zero wavevector and involving either three electric fields , or two electric fields and one atomic displacement . following nunes and gonze  @xcite , our approach makes use of the @xmath1 theorem applied to an electric - field - dependent energy functional  @xcite . we report the local density approximation ( lda ) expressions , as implemented within the abinit package  @xcite . our paper is organized as follows . in sec . [ sec_form ] , we describe the theoretical background related to the @xmath1 theorem and the electric field perturbation . in sec . [ sec_nlo ] , we describe the computation of the non - linear optical susceptibilities , the non - resonant raman scattering efficiencies of both transverse ( to ) and longitudinal ( lo ) zone - center optical phonons and the linear electrooptic ( eo ) tensor . in sec . [ sec_res ] , we illustrate the validity of the formalism by applying our methodology to some semiconductors and ferroelectric oxides and we briefly discuss the effect of a scissors correction on the eo coefficients and non - linear optical susceptibilities . some of the tensors we consider in this work depend on static electric fields : they include contributions of both the electrons and the ions . other quantities imply only the response of the valence electrons : they are defined for frequencies of the electric field high enough to get rid of the ionic contributions but sufficiently low to avoid electronic excitations . for clarity , we adopt the following convention . static fields will be labeled by greek indices ( @xmath2 , @xmath3 , ... ) while we will refer to optical fields with roman symbols ( @xmath4 , @xmath5 , ... ) . to simplify the notations , we will also drop labels such as `` @xmath6 '' for quantities that do not involve the response of the ions . using this convention , we can write @xmath7 and @xmath8 , respectively , for the optical and static dielectric tensor , respectively , and @xmath9 for the linear eo ( pockels ) tensor that involves two optical and one static electric field . in this section , we present the general framework of the computation of third order energy derivatives based on the @xmath0 theorem  @xcite . using the notations of refs . [ ] , we consider three hermitian perturbations labeled @xmath10 , @xmath11 and @xmath12 . the mixed third - order derivatives @xmath13 can be computed from the ground - state and first - order wavefunctions    @xmath14    with @xmath15 - \sum_{\alpha,\beta } \lambda_{\beta \alpha}^{\lambda_2 } \langle \psi_{\alpha}^{\lambda_1 } | \psi_{\beta}^{\lambda_3 } \rangle \nonumber \\ % & & + % \frac{1}{6 } \int d\textbf{r } d\textbf{r } ' d\textbf{r } '' \frac{\delta^3 e_{hxc}[n^{(0 ) } ] } { \delta n(\textbf{r } ) \delta n(\textbf{r } ' ) \delta n(\textbf{r } '' ) } n^{\lambda_1}(\textbf{r } ) n^{\lambda_2}(\textbf{r } ' ) n^{\lambda_3}(\textbf{r } '' ) \nonumber \\ % & & + % \frac{1}{2 } \int d\textbf{r } d\textbf{r } ' \left . \frac{d}{d \lambda_{2 } } \frac{\delta^2 e_{hxc}[n^{(0 ) } ] } { \delta n(\textbf{r } ) \delta n(\textbf{r } ' ) } \right |_{\lambda = 0 } n^{\lambda_1}(\textbf{r } ) n^{\lambda_3}(\textbf{r } ' ) \nonumber \\ % & & + % \frac{1}{2 } \int d\textbf{r } \left . \frac{d^2}{d \lambda_{1 } d \lambda_{3 } } \frac{\delta e_{hxc}[n^{(0 ) } ] } { \delta n(\textbf{r } ) } \right |_{\lambda = 0 } n^{\lambda_2}(\textbf{r } ) + \left . \frac{1}{6 } \frac{d^3 e_{hxc}[n^{(0 ) } ] } { d \lambda_{1 } d \lambda_{2 } d \lambda_{3 } } \right |_{\lambda = 0}.\end{aligned}\ ] ]    @xmath16 is the kinetic energy and @xmath17 ( resp . @xmath18 ) is the sum of the hartree and exchange - correlation energy ( resp . potential ) . the first - order potential @xmath19 can be computed as a second - order functional derivative of @xmath17  @xcite @xmath20}{\delta n(\bm{r } ) \delta n(\bm{r } ' ) } n^{\lambda_2}(\bm{r } ' ) \ ; d \bm{r } ' + \left . \frac{d}{d \lambda_{2 } } \frac{\delta e_{hxc}[n^{(0)}]}{\delta n(\bm{r } ) } \right |_{\lambda = 0}.\ ] ] within the parallel gauge , the first - order lagrange multipliers are given by @xmath21 as a consequence of the `` 2n + 1 '' theorem , the evaluation of eq . ( [ eq_3dte_exp ] ) requires no higher order derivative of the wavefunctions than the first one . these first - order wavefunctions are nowadays available in several first - principles codes . they can be computed from linear response by minimizing a stationary expression of the second - order energy  @xcite or equivalenty by solving the corresponding sternheimer equation  @xcite . it follows that the computation of third - order energies does not require additional quantities than the calculation of second - order energy derivatives . ( [ eq_3dte_exp ] ) is the general expression of the third - order energy derivatives . in case at least one of the perturbations does not affect the explicit form of the kinetic energy or the hartree and exchange - correlation energy , it can be simplified : some of the terms may be zero . this is the case for the electric field perturbations treated in this work as well as for phonon type perturbations . further simplifications can be made in case pseudopotentials without non - linear exchange - correlation core - correction are used . as mentioned in the introduction , special care is required in case one of the perturbations @xmath22 is a macroscopic electric field @xmath23 . in fact , as discussed in the litterature , for infinite periodic solids , usually treated with born - von karman boundary conditions , the scalar potential @xmath24 breaks the periodicity of the crystal lattice . moreover , it is unbound from below  : it is always possible to lower the energy by transferring electrons from the valence states to the conduction states in a distant region ( zener breakdown ) . however , for sufficiently small fields , the tunneling current through the band gap can be neglected and the system is well described by a set of electric field dependent wannier functions @xmath25 . as shown by nunes and vanderbilt  @xcite , these wannier functions minimize the energy functional @xmath26 = e_0 \left [ w_n \right ] - \omega \bm{\mathcal{e } } \cdot \bm{\mathcal{p}}\ ] ] where @xmath27 is the kohn - sham energy under zero field and @xmath28 the macroscopic polarization that can be computed from the wannier function centers . it is important to note that these wannier functions do not correspond to the true ground - state of the system but rather to a long lived metastable state .    in practical applications , it is not mandatory to evaluate the functional eq . ( [ eq_energy ] ) in a wannier basis . it can equivalently be expressed using bloch functions @xmath29 related to @xmath30 by a unitary transform . in this case , the polarization can be computed as a berry phase of the occupied bands  @xcite @xmath31 where bz is the brillouin zone , @xmath32 the absolute value of the electronic charge and the factor of 2 accounts for spin degeneracy . the bloch functions are chosen to satisfy the periodic gauge condition @xmath33 in order to use eq . ( [ eq_ksvb ] ) in practical calculations , the integration over the bz , as well as the differentiation with respect to * k * , have to be performed on a discrete mesh of @xmath34 * k*-points . in case of the ground - state polarization , the standard approach is to build strings of * k*-points parallel to a vector of the reciprocal space @xmath35 . the polarization can then be computed as a string - averaged berry phase . unfortunately , the adaptation of this method to the computation of the energy derivatives is plagued with several difficulties , like the following . the general form of the non - linear optical susceptibility tensor of a compound is imposed by its symmetry . for example , in zinc - blende semiconductors , this tensor , expressed in _ cartesian coordinates _ reduces to @xmath36 where @xmath37 is the levi - civita tensor . it follows that the _ reduced coordinates _ formulation of @xmath38 satisfy the relation @xmath39 where at least one of the three indices @xmath4 , @xmath5 and @xmath40 are different from the two others . when we tried to use strings of * k*-points to compute @xmath38 , eq . ( [ eq_chizb ] ) was not satisfied . however , we were able to avoid these problems , by using the finite difference formula of marzari and vanderbilt  @xcite on a regular grid of special * k*-points ( instead of strings ) @xmath41\ ] ] where * b * is a vector connecting a * k*-point to one of its nearest neighbours and @xmath42 is a weight factor . the sum in eq . ( [ eq_mv ] ) includes as many shells of first neighbours as necessary to satisfy the condition @xmath43 where @xmath44 are the reduced coordinates of * b * and @xmath45 is the metric tensor associated to the real space crystal lattice .    in case of the ground - state polarization , we can not apply the discretization eq . ( [ eq_mv ] ) directly to eq . ( [ eq_ksvb ] ) . as shown by marzari and vanderbilt , the discretization of eq . ( [ eq_ksvb ] ) does not transform correctly under the gauge transformation @xmath46 since eq . ( [ eq_gauge ] ) is equivalent to a shift of the origin by one lattice vector @xmath47 , @xmath28 must change accordingly by one polarization quantum . in order to obtain a discrete expression that matches this requirement , we must combine eq . ( [ eq_mv ] ) with the king - smith and vanderbilt formula  @xcite @xmath48\ ] ] where @xmath49 is the overlap matrix between bloch functions at * k * and @xmath50 : @xmath51    as discussed by nunes and gonze  @xcite , when we apply the perturbation expansion of the preceeding section to the energy functional eq . ( [ eq_energy ] ) , we can adopt two equivalent approaches . the first possibility is the use of eq . ( [ eq_ksvb ] ) for the polarization and a discretization after having performed the perturbation expansion . the second possibility is to apply the @xmath0 theorem directly to eq . ( [ eq_poldisc ] ) in which case no additional discretization is needed . using the notations of nunes and gonze , we will refer to the first case as the dape ( discretization after perturbation expansion ) formulation and to the second one as the pead ( perturbation expansion after discretization ) formulation of the third - order energy . in the following sections , we will discuss both expressions . in addition , in sec . [ sec_res_sc ] , we will compare their convergence with respect to the * k*-point sampling on a realistic example . the perturbation expansion of the first term ( @xmath27 ) of eq . ( [ eq_energy ] ) can easily be performed as it is described in the sec . [ sec_der ] . at the opposite , the expansion of the second term ( @xmath52 ) is more tricky since it explicitely depends on the polarization . in the two sections that follow , we will focus on the @xmath52 term of eq . ( [ eq_energy ] ) . it will be referred to as @xmath53 .      according to the formalism of the preceeding section , the @xmath54 term acts as an additional external potential that has to be added to the ionic one . the @xmath54 perturbation is linear in the electric field and does not depend explicitely on other variables such as the atomic positions . it just enters the terms of eq . ( [ eq_3dte_exp ] ) that involve the first derivative of @xmath55 with respect to @xmath23 . in other words , the only terms in eq . ( [ eq_perm ] ) that involve the expansion of @xmath28 are of the form @xmath56 where @xmath57 and @xmath58 represent an arbitrary perturbation such as an electric field or an atomic displacement . the dape expression of the third - order derivative of e@xmath59 writes as follows @xmath60 are the projection of the first - order wavefunctions on the conduction bands . the complete expression of various third - order energy derivatives , taking into account the expansion of both e@xmath61 and e@xmath59 , are reported in sec . [ sec_nlo ] . ( [ eq_dape ] ) was derived first by dal corso and mauri  @xcite in a slightly different context : they performed the perturbation expansion of the energy functional eq . ( [ eq_energy ] ) using a wannier basis . the resulting expression of the third - order energy was expressed in terms of bloch functions by applying a unitary transform to the wannier orbitals . using the finite difference expression of marzari and vanderbilt eq . ( [ eq_mv ] ) , eq . ( [ eq_dape ] ) becomes @xmath62 where @xmath63 is a basis vector of the reciprocal lattice . applying directly the @xmath0 theorem to eq . ( [ eq_poldisc ] ) we obtain the alternative pead formulation of the third - order energy : @xmath64\end{aligned}\ ] ] where @xmath65 is the inverse of the overlap matrix @xmath49 and @xmath66 its first - order perturbation expansion @xmath67 in the previous section we have discussed the general expressions of third energy derivatives . we now particularize them to the computation of selected non - linear properties .      in an insulator the polarization can be expressed as a taylor expansion of the macroscopic electric field @xmath68 where @xmath69 is the zero - field ( spontaneous ) polarization , @xmath70 the linear dielectric susceptibility ( second rank tensor ) and @xmath38 the second - order non - linear optical susceptibility ( third rank tensor ) . in the litterature on non - linear optics , one often finds another definition of the non - linear optical susceptibility : instead of @xmath38 , it is more convenient to rely on the @xmath71-tensor defined as @xmath72    in general , the polarization depends on valence electrons as well as ions . in the present section , we deal only with the electronic contribution : we will consider the ionic cores as clamped to their equilibrum positions . this constraint will be relaxed in the following sections where we allow for ionic displacements . experimentally , the electronic contribution to the linear and non - linear susceptibilities corresponds to measurements for electric fields at frequencies high enough to get rid of the ionic relaxation but low enough to avoid electronic excitations . in case of the second - order susceptibilities , this constraint implies that both the frequency of @xmath23 , and its second harmonic , are lower than the fundamental absorption gap . the general expression of the electronic non - linear optical susceptibility depends on the frequencies of the optical electric fields [ see for example ref . [ ] ] . in the present context of the @xmath0 theorem applied within the lda to ( static ) dft , we neglect the dispersion of @xmath38 . as a consequence , @xmath38 satisfies kleinman s  @xcite symmetry condition which means that it is symmetric under a permutation of @xmath4 , @xmath5 and @xmath40 . in order to be able to investigate its frequency dependence , one would need either to apply the formalism of time - dependent dft  @xcite or to use expressions that involve sums over excited states  @xcite . following the work of dal corso and co - workers  @xcite we can relate the non - linear optical susceptibilities to a third - order derivative of the energy with respect to an electric field @xmath73 where @xmath74 is defined as the sum over the permutations of the three perturbations @xmath75 ( [ eq_perm ] ) . using the pead formulation of sec . [ sec_elfd ] we can express these terms as follows    @xmath76 \nonumber \\ % & & % + \frac{2}{n_k } \sum_{\bm{k } } \sum_{n , m}^{occ } \bigglb [ \delta_{m , n } \langle u_{n \textbf{k}}^{\mathcal{e}_{i } }   | v_{hxc}^{\mathcal{e}_{j } }    - \langle u_{m \textbf{k}}^{(0 ) } | v_{hxc}^{\mathcal{e}_{j } }    \langle u_{n \textbf{k}}^{\mathcal{e}_{i } } | u_{m \textbf{k}}^{\mathcal{e}_{l } }   \rangle \bigglb ] \nonumber \\ % & & + % \frac{1}{6 } \int d\textbf{r } d\textbf{r } ' d\textbf{r } '' \frac{\delta^3 e_{xc}[n^0 ] } { \delta n(\textbf{r } ) \delta n(\textbf{r } ' ) \delta n(\textbf{r } '' ) } n^{\mathcal{e}_i}(\textbf{r } ) n^{\mathcal{e}_j}(\textbf{r } ' ) n^{\mathcal{e}_l}(\textbf{r}'').\end{aligned}\ ] ]      we now consider the computation of raman scattering efficiencies of zone - center optical phonons . in the limit @xmath77 , the matrix of interatomic force constants @xmath78 can be expressed as the sum of an analytical part and a non - analytical term  @xcite @xmath79 the analytical part corresponds to the second - order derivative of the energy with respect to an atomic displacement at @xmath80 under the condition of vanishing macroscopic electric field . the second term is due to the long - range electrostatic interactions in polar crystals . it is at the origin of the so - called lo - to splitting and can be computed from the knowledge of the born effective charges @xmath81 and the electronic dielectric tensor  @xcite @xmath7 . the phonon frequencies @xmath82 and eigendisplacements @xmath83 are solution of the following generalized eigenvalue problem @xmath84 where @xmath85 is the mass of atom @xmath86 . as a convention , we choose the eigendisplacements to be normalized as @xmath87    in what follows we consider non - resonant raman scattering where an incoming photon of frequency @xmath88 and polarization @xmath89 is scattered to an outgoing photon of frequency @xmath90 and polarization @xmath91 by creating a phonon of frequency @xmath82 ( stokes process ) . the scattering efficiency  @xcite ( cgs units ) corresponds to @xmath92 where @xmath93 is the speed of light in vacuum and @xmath94 the boson factor @xmath95 the raman susceptibility @xmath96 is defined as @xmath97 where @xmath98 is the electronic linear dielectric susceptibility tensor . @xmath99 is the angle of collection in which the outgoing photon is scattered . due to snell s law , @xmath99 is modified at the interface between the sample and the surrounding medium . experimentally , the scattering efficiencies are measured with respect to the solid angle of the medium while eq . ( [ eq_efficiency ] ) refers to the solid angle inside the sample . in order to relate theory and experiment , one has to take into account the different refractive indices of the sample and medium . for example , in the isotropic case , eq . ( [ eq_efficiency ] ) has to be multiplied  @xcite by @xmath100 where @xmath101 and @xmath102 are respectively the refractive indices of the sample and the medium .    for pure transverse optical phonons , @xmath103 can be computed as a mixed third - order derivative of the energy with respect to an electric field , twice , and to an atomic displacement under the condition of zero electric field @xmath104 in case of longitudinal optical phonons with wavevector @xmath77 in a polar crystal , eq . ( [ eq_ramansus ] ) must take into account the effect of the macroscopic electric field generated by the lattice polar vibration . this field enters the computation of the raman susceptibilities at two levels . on one hand , it gives rise to the non - analytical part of the matrix of interatomic force constants eq . ( [ eq_ifc ] ) that modifies the frequencies and eigenvectors with respect to pure transverse phonons . on the other hand , the electric field induces an additional change in the dielectric susceptibility tensor related to the non - linear optical coefficients @xmath105 . for longitudinal optical phonons , eq . ( [ eq_dchidtau ] ) has to be modified as follows  @xcite : @xmath106    the mixed third - order derivatives ( [ eq_dchidtau ] ) can be computed from various techniques including finite differences of the dielectric tensor  @xcite or the second derivative of the electronic density matrix  @xcite . here , we follow an approach similar to deinzer and strauch  @xcite based on the @xmath1 theorem . the third - order energy can be computed as the sum over the 6 permutations eq . ( [ eq_perm ] ) of @xmath107 , @xmath108 and @xmath109 . according to the discussion of sec . [ sec_elfd ] , we have to distinguish between the terms that involve the discretization of the polarization such as @xmath110 or @xmath111 and those that can be computed from a straightforward application of the @xmath1 theorem such as @xmath112 . the former ones show an electric field as second perturbation . they can be computed from an expression analogous to eq . ( [ eq_e3_elfd ] )    @xmath113 \nonumber \\ % & & % + \frac{2}{n_k } \sum_{\bm{k } } \sum_{n , m}^{occ } \bigglb [ \delta_{m , n } \langle u_{n \textbf{k}}^{\tau_{\kappa \lambda } }   | v_{hxc}^{\mathcal{e}_{i } }    - \langle u_{m \textbf{k}}^{(0 ) } | v_{hxc}^{\mathcal{e}_{i } }    \langle u_{n \textbf{k}}^{\tau_{\kappa \lambda } } | u_{m \textbf{k}}^{\mathcal{e}_{j } }   \rangle \bigglb ] \nonumber \\ % & & + % \frac{1}{6 } \int d\textbf{r } d\textbf{r } ' d\textbf{r } '' \frac{\delta^3 e_{xc}[n^0 ] } { \delta n(\textbf{r } ) \delta n(\textbf{r } ' ) \delta n(\textbf{r } '' ) } n^{\tau_{\kappa \lambda}}(\textbf{r } ) n^{\mathcal{e}_{i}}(\textbf{r } ' ) n^{\mathcal{e}_{j}}(\textbf{r}'').\end{aligned}\ ] ]    we obtain a similar expression for @xmath111 . the remaining terms do not require any differentiation with respect to * k*. they can be computed from the first - order change of the ionic ( pseudo ) potential with respect to an atomic displacement @xmath114 @xmath115 \nonumber \\ % & & % + \frac{1}{2 } \int d\textbf{r } d\textbf{r } ' \left . \frac{d}{d \tau_{\kappa \lambda } } \frac{\delta^2 e_{hxc } } { \delta n(\textbf{r } ) \delta n(\textbf{r } ' ) } \right |_{n^{(0 ) } } n^{\mathcal{e}_{i}}(\textbf{r } ) n^{\mathcal{e}_{j}}(\textbf{r } ' ) \nonumber \\ % & & % \frac{1}{6 } \int d\textbf{r } d\textbf{r } ' d\textbf{r } '' \frac{\delta^3 e_{xc}[n^0 ] } { \delta n(\textbf{r } ) \delta n(\textbf{r } ' ) \delta n(\textbf{r } '' ) } n^{\tau_{\kappa \lambda}}(\textbf{r } ) n^{\mathcal{e}_{i}}(\textbf{r } ' ) n^{\mathcal{e}_{j}}(\textbf{r}'').\end{aligned}\ ] ]    in pseudopotential calculations , the computation of the first - order ionic potential @xmath114 requires the derivative of local and non - local ( usually separable ) operators . these operations can be performed easily without any additional workload as described in ref . [ ] .    in spite of the similarities between eqs . ( [ eq_e3_tee ] ) and ( [ eq_e3_ete ] ) and the expression proposed by deinzer and strauch we can quote few differences . first , our expression of the third - order energy makes use of the pead fomulation for the expansion of the polarization . moreover , eq . ( [ eq_e3_ete ] ) is more general since it allows the use of pseudopotentials with non - linear core correction through the derivative of the second - order exchange - correlation energy with respect to @xmath107 ( third term ) .      as in the cases of the born effective charges and of the dynamical matrix  @xcite , the coefficients @xmath116 must vanish when they are summed over all atoms in the unit cell . @xmath117 physically , this sum rule guarantees the fact that the macroscopic dielectric susceptibility remains invariant under a rigid translation of the crystal . in practical calculations , it is not always satisfied although the violation is generally less severe than in case of @xmath78 or @xmath118 . even in calculations that present a low degree of convergence , the deviations from this law can be quite weak . they can be corrected using similar techniques as in case of the born effective charges  @xcite . for example , we can define the mean excess of @xmath116 per atom @xmath119 and redistribute it equally between the atoms @xmath120      the optical properties of a compound usually depend on external parameters such as the temperature , electric fields or mechanical constraints ( stress , strain ) . in the present section we consider the variations of the refractive index induced by a static or low - frequency electric field @xmath121 . at the linear order , these variations are described by the linear eo coefficients ( pockels effect ) @xmath122 where @xmath123 is the inverse of the electronic dielectric tensor and @xmath9 the eo tensor .    within the born and oppenheimer approximation , the eo tensor can be expressed as the sum of three contributions : a bare electronic part @xmath124 , an ionic contribution @xmath125 and a piezoelectric contribution @xmath126 . the electronic part is due to an interaction of @xmath121 with the valence electrons when considering the ions artificially as clamped at their equilibrium positions . it can be computed from the non - linear optical coefficients . as can be be seen from eq . ( [ eq_nlo_pol ] ) , @xmath127 defines the second - order change of the induced polarization with respect to @xmath121 . taking the derivative of eq . ( [ eq_nlo_pol ] ) , we also see that @xmath127 defines the first - order change of the linear dielectric susceptibility , which is equal to @xmath128 since the eo tensor depends on @xmath129 rather than @xmath130 , we have to transform @xmath130 by the inverse of the zero field electronic dielectric tensor  @xcite @xmath131 using eq . ( [ eq_transfo ] ) we obtain the following expression for the electronic eo tensor @xmath132 eq . ( [ eq_eo_el ] ) takes a simpler form when expressed in the principal axes of the crystal under investigation  @xcite @xmath133 where the @xmath134 coefficients are the principal refractive indices . the origin of the ionic contribution to the eo tensor is the relaxation of the atomic positions due to the applied electric field @xmath121 and the variations of @xmath7 induced by these displacements . it can be computed from the born effective charges @xmath135 and the @xmath136 coefficients introduced in sec . [ sec_raman ] . as shown in appendix [ sec_app ] [ see also refs . ] , the ionic eo tensor can be computed as a sum over the transverse optic phonon modes at @xmath80 @xmath137 where @xmath96 is the raman susceptibility of mode @xmath138 [ eq . ( [ eq_ramansus ] ) ] and @xmath139 the mode polarity @xmath140 that is directly linked to the modes oscillator strength @xmath141 for simplicity , we have expressed eq . ( [ eq_eo_ion ] ) in the principal axes while a more general expression can be derived from eq . ( [ eq_transfo ] ) . finally , the piezoelectric contribution is due to a relaxation of the unit cell shape due to the converse piezoelectric effect  @xcite . as it is discussed in appendix [ sec_app ] , it can be computed from the elasto - optic coefficients @xmath142 and the piezoelectric strain coefficients @xmath143 @xmath144    in the discussion of the eo effect , we have to specify whether we are dealing with strain - free ( clamped ) or stress - free ( unclamped ) mechanical boundary conditions . the clamped eo tensor @xmath145 takes into account the electronic and ionic contributions but neglects any modification of the unit cell shape due to the converse piezoelectric effect  @xcite @xmath146 experimentally , it can be measured for frequencies of @xmath121 high enough to eliminate the relaxations of the crystal lattice but low enough to avoid excitations of optical phonon modes ( usually above @xmath147 100 mhz ) . to compute the unclamped eo tensor @xmath148 , we have to add the piezoelectric contribution to @xmath145 @xmath149 experimentally , @xmath148 can be measured for frequencies of @xmath121 below the ( geometry dependent ) mechanical body resonances of the sample  @xcite ( usually below @xmath147 1 mhz ) . our calculations have been performed within the local density approximation ( lda ) to the density functional theory  @xcite ( dft ) . we used the abinit  @xcite package , a planewave , pseudopotential dft code  @xcite in which we have implemented the formalism presented above . for reasons that will become obvious below , we chose the pead formulation eq . ( [ eq_pead ] ) to perform the differentiation with respect to * k*. for the exchange - correlation energy @xmath150 we relied on the parametrization of perdew and wang  @xcite as well as the parametrization of goedecker , teter and hutter  @xcite . these expressions have the advantage to avoid any discontinuities in the functional derivative of @xmath150 .    in case of the semiconductors si , alas and alp , we used at @xmath151 grid of special k - points , a plane - wave kinetic energy cutoff of 10 hartree and troullier - martins  @xcite norm - conserving pseudopotentials . these calculations have been performed at the theoretical lattice constant . to perform the finite difference calculations of the raman polarizabilities , changes of the electronic dielectric tensor were computed for atoms displaced by @xmath152 of the unit cell parameter along the cartesian directions .    in case of rhombohedral batio@xmath153 , we used a @xmath154 grid of special k - points , a plane - wave kinetic energy cutoff of 45 hartree and extended norm - conserving pseudopotentials  @xcite . since the ferroelectric instability is quite sensitive to the volume of the unit cell and tends to disappear due to the volume underestimation of the lda  @xcite , we chose to work at the experimental lattice constants . at the opposite to the lattice parameters , the atomic positions have been relaxed : the residual forces on the atoms were smaller than @xmath155 hartree / bohr . it was shown by gonze , ghosez and godby  @xcite that an accurate functional for the exchange - correlation energy in extended systems should depend on both the density and the polarization . the lda used here neglects this polarization dependence and may consequently indroduce significant relative errors when studying the response of a solid to an electric field . in case of the second - order derivatives , the lda usually yields an overestimate of the dielectric tensor ( as large as 20 % in batio@xmath156 )  @xcite . at the opposite , no clear trends have been reported yet concerning non - linear optical properties such as @xmath127 . @xcite    in lda calculations , it is common practice to apply a _ scissors correction _ @xcite to compensate the lack of polarization dependence of the exchange - correlation functional . in case of non - linear optical properties , such a correction can be applied at different levels . on the one hand , we can compute the non - linear optical susceptibilities ( eq . ( [ eq_e3_elfd ] ) ) using a scissors operator for the first - order wavefunctions  @xcite . on the other hand , in the computation of the eo coefficients , we can use a scissors corrected refractive index in eqs ( [ eq_eo_el1 ] ) and ( [ eq_eo_ion ] ) . the influence of these corrections will be discussed below .      in order to illustrate the computation of third - order energy derivatives , we performed a series of calculations on various cubic semiconductors . in these compounds , the non - linear optical susceptibility tensor @xmath157 and the raman susceptibility tensor @xmath158 only have one independent element d@xmath159 and @xmath160 . also , instead of @xmath160 it is customary to report the raman polarizability  @xcite defined as @xmath161 where @xmath162 is the reduced mass of the two atoms in the unit cell .    the formalism of sec . [ sec_form ] involves an integration over the bz and a differentiation with respect to * k*. in practical calculations , these operations must be performed on a discrete mesh of special * k*-points . as we explained in sec . [ sec_form ] , the discretization can either be performed before ( pead ) or after ( dape ) the perturbation expansion of the energy functional eq . ( [ eq_energy ] ) . up to know , the applications of the present formalism to real materials  @xcite made use of the dape formula of the third - order energy . the only application of the pead formula has been reported by nunes and gonze  @xcite on a one - dimensional model system . these authors observed that the pead formula converges better with respect to the * k*-point sampling than the dape formula . in order to compare the performance of these two approaches on a realistic case , we applied both of them to compute the non - linear optical susceptibility d@xmath159 of alas . we performed a series of calculations on a @xmath163 grid of special * k*-points . as can be seen on figure [ fig_convergence ] the pead formula converges much faster than the dape formula . therefore , the pead formulation has been applied to obtain the results presented below . it is the one that is actually available in the abinit code .     non - linear optical susceptibility d@xmath159 ( pm / v ) of alas for various grids of @xmath163 special * k*-points.,width=226 ]    in table [ tab_scnl ] , we report the non - linear optical susceptibilities of the cubic semiconductors alas and alp . our results are in close agreement with the values obtained by dal corso and co - workers  @xcite who applied the @xmath0 theorem within the dape formalism , the results of levine and co - workers  @xcite who used a `` sum over excited states '' technique and the values obtained by souza and co - workers  @xcite who followed a finite electric field approach . the values in the lower part of table [ tab_scnl ] have been obtained using a scissors correction . our methodology provides a correction similar to what is reported by levine and allan  @xcite . the scissors correction decreases the value of the non - linear optical susceptibilities in agreement with the discussion of ref . . to the authors knowledge , no experimental data are available for alas and alp . for other cubic semiconductors , it is however not clear that the use of a scissors correction improves the agreement with the experiment  @xcite and will even have a negative effect when the lda underestimates the experimental value . in addition , it is not straightforward to isolate the error of the lda on the non - linear response functions from other sources of errors . other factors have a similar strong influend on @xmath127 as the scissors correction . for example , the values of the non - linear optical susceptibilities strongly depend on the pseudopotential  @xcite or on the error on the unit cell volume  @xcite that is usually underestimated in lda calculations . .[tab_scnl ] non - linear optical susceptibilities @xmath164 ( pm / v ) of some semiconductors . the values in the lower part of the table have been obtained using a scissors ( sci ) correction . [ cols="<,<,^,^",options="header " , ] in this paper , we presented the general framework for the computation of third - order energy derivatives within dft . our formalism makes use of the @xmath165 theorem and the modern theory of polarization . focusing on derivatives that are characterized by a zero wavevector and that involve either three electric fields or two electric fields and one atomic displacement , we described the computation of non - linear optical susceptibilities , of raman scattering efficiencies of to and lo phonons and of the eo tensor . the computation of the berry phase polarization involves a derivative of the wavefunctions with respect to their wavevector . in practice , this differentiation is computed on a grid of special * k*-points . the perturbation expansion can either be performed before ( dape ) or after ( pead ) the discretization , leading to two mathematically distinct expression of the third - order energies . we used both of them to compute the non - linear optical susceptibility of alas and we have shown that the pead formulation converges faster with respect to the * k*-point sampling .    we have computed the non - linear optical susceptibilities and raman polarizabilities of some cubic semiconductors as well as the eo tensor in the rhombohedral phase of batio@xmath153 .    finally , we have studied the effect of a scissors correction on the eo coefficients and the non - linear optical susceptibilities . at the opposite to the dielectric tensor , we did not find a systematic improvement of the results by using this correction . we can figure out several applications of the methodology presented in this work . combined with the calculation of phonon frequencies and infrared intensities , the computation of raman efficiencies can be a useful complementary tool for the interpretation of experimental spectra . furthermore , the computation of the eo tensor from first - principles can guide the tuning of the eo properties and help designing new efficient eo materials . this could reveal particularly helpful since accurate optical measurements require high quality single crystals not always directly accessible . the authors are grateful to m. d. fontana , p. bourson , b. kirtman and b. champagne for helpful discussions . m. v. and x. g. acknowledge financial support from the fnrs belgium . this work was supported by the volkswagen - stiftung within the project nano - sized ferroelectric hybrids `` ( i/77 737 ) , the region wallonne ( nomade , project 115012 ) , the communaut francaise de belgique - actions de recherche concertes , the pai / uiap phase 5 ' ' quantum size effects in nanostructured materials `` , fnrs - belgium through grants 9.4539.00 and 2.4562.03 , the european union through the research and training network ' ' exciting " ( hprm - ct-2002 - 00317 ) . as it is discussed in sec . [ sec_eo ] , the optical properties of a compound are modified by an electric field @xmath121 or a mechanical constraint ( a stress @xmath166 or a homogeneous strain @xmath167 ) . at the linear order , the variations of @xmath168 can be described using either the variables @xmath169 or @xmath170  @xcite      where @xmath145 and @xmath148 are resp . the clamped ( strain free ) and unclamped ( stress free ) eo coefficients , @xmath142 are the elasto - optic ( strain - optic ) coefficients and @xmath172 are the piezo - optical ( stress - optical ) coefficients . in order to relate eqs . ( [ eq_reta ] ) and ( [ eq_rsigma ] ) , we can express the strain as beeing induced by the stress or by the electric field ( converse piezoelectric effect ) @xmath173 where @xmath174 are the elastic compliances and @xmath143 the piezoelectric strain coefficients . if we assume for example that the unit cell is free to relax within the electric field ( stress - free mechanical boundary conditions ) we can either use eq . ( [ eq_rsigma ] ) ( in which case the second term of the right - hand - side is zero ) or eq . ( [ eq_reta ] ) to compute @xmath175 . in the latter case , the strain induced by the electric field can be obtained from the second term of the right - hand - side of eq . ( [ eq_strain ] ) @xmath176 using this identity , we obtain the following relation between the unclamped and the clamped eo coefficients @xmath177      in order to derive the expressions of the clamped and unclamped eo tensor of sec . [ sec_eo ] , we use a taylor expansion of the electric enthalpy  @xcite @xmath178 . similar developments have already been applied to determine the lattice contribution of the static dielectric tensor and of the piezoelectric tensor  @xcite . they are based on an expansion of @xmath178 up to the second order in the atomic coordinates @xmath179 , the homogeneous strain @xmath167 and the macroscopic electric field @xmath121 . in this section , we extend these developments to the third order . the electric enthalpy of a solid in an electric field is obtained by the minimization @xmath180 we denote @xmath181 , @xmath182 the atomic positions and the strain that minimize @xmath178 at constant @xmath23 and @xmath183 , @xmath184 (= 0 ) their values at @xmath185 . for small fields , we can expand the function @xmath186 in powers of @xmath23 around @xmath185 : @xmath187 where @xmath99 is the volume of the primitive unit cell in real space and @xmath188 , @xmath189 and @xmath190 are the macroscopic polarization , electronic dielectric tensor and non - linear optical coefficients at zero macroscopic electric field and for a given configuration ( @xmath191 , @xmath192 ) . at non - zero field , these quantities are defined as partial derivatives of @xmath178 with respect to @xmath23 . for example , the electric field dependent electronic dielectric tensor can be computed from the expression @xmath193    let @xmath194 be the displacement of atom @xmath86 along direction @xmath2 and @xmath195 ( resp . @xmath196 ) the first - order modification of the atomic position ( resp . strain ) induced by a perturbation @xmath197 @xmath198 in the discussion that follows , we will study the effect of an electric field perturbation and a strain perturbation on the electric enthalpy @xmath178 in order to obtain the formulas to compute the elasto - optic coefficients as well as the clamped and the unclamped eo tensors . the elasto - optic tensor can be computed from the _ total _ derivative of the dielectric tensor with respect to @xmath167 at zero electric field @xmath200 the derivative in the first term of the right hand side is computed considering the ionic cores as artificially clamped at their equilibrum positions . the remaining terms represent the ionic contribution to the elasto - optic tensor . they involve derivatives of the linear dielectric susceptibility @xmath98 with respect to the atomic positions that have to be multiplied by the first - order strain induced atomic displacements @xmath201 [ eq . ( [ eq_tau ] ) ] . to compute these quantities we use the fact that @xmath178 is minimum at the equilibrum for an imposed strain @xmath192 . this condition implies @xmath202 since we are interested in first - order atomic displacements we can write @xmath203 solving the extremum equation ( [ eq_tau_eq ] ) to the linear order in @xmath192 , we obtain @xmath204 the second derivatives on the left side of eq . ( [ eq_tau_disp ] ) define the matrix of interatomic force constants at zero macroscopic electric field which enables the computation of the transverse phonon frequencies @xmath82 and eigendisplacements @xmath83 . by decomposing @xmath201 in the basis of the zone - center phonon - mode eigendisplacements @xmath205 and using eqs . ( [ eq_eigen ] ) , ( [ eq_norm ] ) we derive the following expression for the first - order strain induced atomic displacements @xmath206 where @xmath207 if we introduce eqs . ( [ eq_ion_decomp ] ) and ( [ eq_tau_disp1 ] ) into eq . ( [ eq_depsdeta ] ) and use the definition of the raman susceptibility eq . ( [ eq_ramansus ] ) and the transformation eq . ( [ eq_transfo ] ) , we finally obtain the formula to compute the elastooptic tensor @xmath208 to simplify , we write eq . ( [ eq_elastoopt ] ) in the principal axes of the crystal under investigation . a more general expression can be obtained from eq . ( [ eq_transfo ] ) . ( [ eq_elastoopt ] ) is different from the approach used previously by detraux and gonze to study the elasto - optic tensor in @xmath2-quartz  @xcite . the authors of ref . used finite differences with respect to strains to compute the the total derivative of @xmath7 . in their approach , the atoms where relaxed to their equilibrum positions in the strained configurations . in case of eq . ( [ eq_elastoopt ] ) , the first term of the right - hand - side is computed at clamped atomic positions while the effect of the strain - induced atomic relaxations is taken into account by the second term . the clamped eo tensor can be computed from the _ total _ derivative of the electric field dependent dielectric tensor eq . ( [ eq_ion_diel ] ) with respect to @xmath210 @xmath211 the derivative in the first term is computed considering the ionic cores as artificially clamped at their equilibrum positions . this term represents the bare electronic contribution to the eo tensor that can be computed from the non - linear optical coefficients @xmath212 related to a third - order partial derivative of @xmath178 @xmath213 the remaining terms in eq . ( [ eq_clamp_chidc ] ) represent the ionic contribution to the eo tensor . they involve derivatives of the linear dielectric susceptibility @xmath98 with respect to the atomic positions that have to be multiplied by the first - order electric field induced atomic displacements @xmath214 [ eq . ( [ eq_tau ] ) ] . to obtain these quantities , we proceed the same way as in case of the elasto - optic tensor . using the equilibrum condition @xmath215 and expanding @xmath216 to the first - order in the electric field , we obtain @xmath217 this expression is similar to eq . ( [ eq_tau_disp ] ) . the second - order derivatives of @xmath178 on the left side are the interatomic force constatnts and the derivative of the zero field polarization with respect to @xmath216 on the right side is the born effective charge tensor @xmath218 of atom @xmath86 . decomposing @xmath214 in the basis of the zone - center phonon - mode eigendisplacements [ eq . ( [ eq_ion_decomp ] ) ] and using the orthononormality constraint eq . ( [ eq_norm ] ) we derive the following expression for the first - order electric field induced atomic displacements @xmath219 if we introduce eqs . ( [ eq_elec_ent ] ) and ( [ eq_clamp_disp1 ] ) into eq . ( [ eq_clamp_chidc ] ) we finally obtain the formula to compute the total derivative of the dielectric tensor @xmath220 using the definition of the raman susceptibility [ eq . ( [ eq_ramansus ] ) ] , the mode polarity [ eq . ( [ eq_eo_pol ] ) ] and the transformation [ eq . ( [ eq_transfo ] ) ] we obtain the expression of the clamped eo tensor @xmath221 as in case of the elasto - optic tensor [ eq . ( [ eq_elastoopt ] ) ] , we have written eq . ( [ eq_clamp_rijgamma ] ) in the principal axes of the crystal under investigation .      in order to compute the unclamped eo tensor , we have to take into account both the electric field induced atomic displacments @xmath214 and the electric field induced strain @xmath223 when computing the total derivative of @xmath7 @xmath224 the electronic contribution [ first term of eq . ( [ eq_unclamp_chidc ] ) ] is the same as for the clamped eo tensor . it can be computed from the non - linear optical coefficients [ eq . ( [ eq_elec_ent ] ) ] . to compute @xmath214 and @xmath223 , we can use an equilibrum condition similar to eq . ( [ eq_clamp_eq ] ) where we require that the first - order derivatives of @xmath178 with respect to @xmath216 and @xmath167 vanish . expanding @xmath216 and @xmath167 to the first - order in the electric field , we obtain the system of coupled equations [ see also ref . ]      because of the coupling between @xmath214 and @xmath223 , defined by the mixed second - order derivatives @xmath226 the second term of the right hand side of eq . ( [ eq_unclamp_chidc ] ) is different from that of eq . ( [ eq_clamp_chidc ] ) . that means that the sum of the first and second term of eq . ( [ eq_unclamp_chidc ] ) is not identical to the clamped eo coefficients @xmath145 . moreover , the third term of eq . ( [ eq_unclamp_chidc ] ) is different from the piezoelectric contribution of sec . [ sec_app_mac ] .    in order to obtain the decomposition of @xmath148 into electronic , ionic and piezoelectric contributions defined previously , we can solve eq . ( [ eq_taueta1 ] ) for @xmath214 . in the basis of the zone - center phonon mode eigendisplacements we can write @xmath227 if we insert this relation into eq . ( [ eq_unclamp_chidc ] ) and use the transformation eq . ( [ eq_transfo ] ) we obtain the following expression of the unclamped eo tensor in the principal axes @xmath228 \eta_{\mu \nu}^{\mathcal{e}_{\gamma}}\end{aligned}\ ] ] the sum of the first and second term of the right - hand side of eq . ( [ eq_unclamp_1 ] ) is equal to the clamped eo coefficient @xmath145 . the product of the conversion factor times the bracket in the third term of eq . ( [ eq_unclamp_1 ] ) is equal to the elasto - optic coefficient @xmath142 [ eq . ( [ eq_elastoopt ] ) ] . finally , by definition of the converse piezoelectric effect , @xmath223 is equal to the piezoelectric strain coefficient @xmath143 . we obtain thus the following expression of the unclamped eo coefficients that is equal to the one derived in sec . [ sec_app_mac ] from pure macroscopic arguments @xmath177 it is worth noticing that the so - called piezoelectric contribution does not only take into account the change of the linear optical susceptibility with strain ( third term of the right hand side of eq . ( [ eq_unclamp_chidc ] ) ) but also includes the modification of the ionic contribution , with respect to the clamped case , that is associated to the modification of the ionic relaxation induced by the strain . x. gonze , j .- m . beuken , r. caracas , f. detraux , m. fuchs , g .- m . rignanese , l. sindic , m. verstraete , g. zerah , f. jollet , m. torrent , a. roy , m. mikami , ph . ghosez , j .- y . raty and d.c . allan , computational materials science * 25 * , 478 ( 2002 ) . [ url www.abinit.org ]
the non - linear response of infinite periodic solids to homogenous electric fields and collective atomic displacements is discussed in the framework of density functional perturbation theory . 
 the approach is based on the @xmath0 theorem applied to an electric - field - dependent energy functional . 
 we report the expressions for the calculation of the non - linear optical susceptibilities , raman scattering efficiencies and electrooptic coefficients . 
 different formulations of third - order energy derivatives are examined and their convergence with respect to the * k*-point sampling is discussed . 
 we apply our method to a few simple cases and compare our results to those obtained with distinct techniques . 
 finally , we discuss the effect of a scissors correction on the eo coefficients and non - linear optical susceptibilities .