Nonlinear optical properties of TeO$_2$ crystalline phases from first principles

We have computed second and third nonlinear optical susceptibilities of two crystalline bulk tellurium oxide polymorphs: $\alpha$-TeO$_{2}$ (the most stable crystalline bulk phase) and $\gamma$-TeO$_{2}$ (the crystalline phase that ressembles the more to the glass phase. Third order nonlinear susceptibilities of the crystalline phases are two orders of magnitude larger than $\alpha$-SiO$_{2}$ cristoballite, thus extending the experimental observations on glasses to the case of crystalline compounds. While the electronic lone pairs of Te contribute to those large values, a full explanation of the anisotropy of the third order susceptibility tensor requires a detailed analysis of the structure, in particular the presence of helical chains, that seems to be linked to cooperative non-local polarizabilty effects. Our results demonstrate that first-principles simulations are a powerful predictive tool to estimate nonlinear optical susceptibilitites of materials.


I. INTRODUCTION
Tellurium oxide glasses arouse lots of interest in the field of nonlinear optics (NLO) since their unusual nonlinear optical indices have been noticed. The third order optical susceptibilities χ (3) exhibited by pure TeO 2 glasses, of the order of 14 × 10 −13 esu, 1 are indeed among the highest observed for oxide glasses (50 times larger than in pure silica glasses) and they thus are of great interest in both fundamental science and technological applications as optical modulators and frequency converters. [1][2][3] The origin of these high values is not fully established yet. Using a combination of experimental techniques (interferometric measurements) and ab initio calculations (within the restricted Hartree-Fock scheme), some authors 4 suggested that the highly polarizable 5s 2 electronic lone pair of Te (IV) could be responsible for the high NLO indices. Other theoretical works, 5,6 based on the use of hybrid functionals within the density functional theory (DFT), reinforced this idea by demonstrating the importance of the Te (IV) lone pair on the hyperpolarizabilities of isolated TeO 4 and TeO 3 structural units. This conclusion was extrapolated, by extension, to the case of TeO 2 based glasses. However, another series of theoretical studies, 7-10 also carried out in the framework of DFT with hybrid functionals on (XO 2 ) n (X = Si or Te) polymer clusters of different shapes (chains, rings and cage geometries) and sizes (monitored by the number n of XO 2 units) suggested another origin for the unusually high values of the NLO susceptibilities, highlighting the relevance of the structural features themselves, in particular how the structural blocks (i. e. TeO 2 units) are linked together. Only one type of such molecules, the linear chains, seem to be capable of realistically reproducing the high hypersusceptibility values for the tellurium oxides. This was attributed to an exceptionally strong nonlocality of the electronic polarization in these chains, much more important for the Te than for the Si oxides.
Nevertheless, previous works were based on hypothetical fragments that were supposed to be likely found in the glass. The question about the high nonlinear susceptibility in the solid phases was not addressed. Clearly, further studies are needed to achieve a deeper understanding in the origin of these properties and notably to gauge the relative importance of the lone pair versus the structural features. A new way to tackle this problem is to treat the case of TeO 2 -based crystalline compounds. This would prevent the recourse to hypothetical structural fragments to feed the first-principles calculations. Besides, it would allow to study how the anisotropic nature in the crystalline phases translates into the variations of dielectric susceptibilities with crystalline directions.
Unfortunately, the situation for crystals is different than for molecules and two main problems arise in the first-principles calculation of the hypersusceptibilities. The optical susceptibilities are derivatives of the bulk macroscopic polarization with respect to electric field. If the polarization can easily be expressed in terms of the charge distribution for molecules (finite systems), it can not be obtained that way for crystals (infinite systems treated periodically). The polarization in a periodic system would indeed depend on the choice of the unit cell. 11 Solutions arised in the early 1990s and are often referred to as the "modern theory of polarization". 12 The basic idea is to consider the change in polarization 13 of a crystal as it undergoes some slow change, e.g. a slow displacement of one sublattice relative to the others, and relate it to the current that flows during this adiabatic evolution of the system. 14 The second problem lies in the nature of the applied electric field that is macroscopic. The scalar potential of a macroscopic homogeneous electric field is non-periodic (so the Bloch theorem does not apply), and unbounded from below (so the energy of the system can always be lowered transferring electrons to regions sufficiently far away, hampering the aplicability of traditional variational methods). 15 The first approach to circumvent this problem in first-principles simulations, due to Kunc and Resta,16 was to consider "sawtooth" potentials in a supercell. We have previously tested this scheme, as implemented in the Crystal06 program, 17 on the computation of hypersusceptibilities of TeO 2 crystalline oxides. 18 Although this method gives satisfactory results, it requires defining a supercell for keeping the periodicity along the applied field direction. The dimension of the studied system is thus very quickly limiting in terms of expansiveness in time and computational requirements. A more recent variational alternative, firmly rooted on the modern theory of polarization, was due to Souza,Íñiguez and Vanderbilt. 15,19 It is based on the minimization of a electric enthalpy functional with respect to a set of polarized Bloch functions, thus including the effect of the electric field directly inside the unit cell. This approach, recently implemented in the Siesta method 20,21 is the one used in the present work.
In this paper we compute the second and third order optical susceptibility tensors in two bulk TeO 2 polymorphs: the α-TeO 2 phase (known as paratellurite) that is the most stable one, and the γ-TeO 2 phase, whose structure ressembles the more to the glass. The estimations of the nonlinear susceptibility data provided in the present work intend to fill the gap in the reported values of these quantities. Unfortunately, there is a cruel lack of experimental nonlinear susceptibility data for crystalline phases, mainly due to the difficulty of growing sufficiently large single crystals. To our knowledge, only the χ (2) susceptibility tensor elements for the α-TeO 2 phase has been measured, while there are no experimental values of the χ (3) tensor elements for any crystalline phase. In addition, third order susceptibility tensor of α-SiO 2cristobalite, which is structurally similar to α-TeO 2 , is computed for comparison.
The rest of the paper is organized as follows. After presenting the computational details in Sec. II, and the structure characteristics of the different polymorphs in Sec. III, we describe the methodology used to compute the nonlinear susceptibilities in Sec. IV. Second-order susceptibility values are then calculated in Sec. V A, and compared to experimental results in order to test the validity and the limitations of the method. Finally, the method is used as a predictive tool through the calculation of the third order susceptibilities in Sec. V B, and clues are given for exploring relevant features responsible for large variations of the dielectric susceptibilities with crystalline directions.

II. COMPUTATIONAL DETAILS
We have carried out density functional first-principles simulations based on a numerical atomic orbital method as implemented in the Siesta code. 20 All the calculations have been carried out within the Generalized Gradient Approximation (GGA), using the functional parametrized by Perdew, Burke and Ernzerhof (PBE) 22 to simulate the electronic exchange and correlation. Core electrons were replaced by ab initio norm conserving pseudopotentials, generated using the Troullier-Martins scheme, 23 in the Kleinman-Bylander fully nonlocal separable representation. 24 The 5s and 5p electrons of Te, 2s and 2p electrons of O, and 3s and 3p electrons of Si were considered as valence electrons and explicitly included in the simulations. In order to avoid the spiky oscillations close to the nucleus that oftenly appear in GGA-generated pseudopotentials, we have included small partial core corrections 25 for all the atoms. Te pseudopotential was generated scalar relativistically. The reference configuration, cutoff radii for each angular momentum shell, and the matching radius between the full core charge density and the partial core charge density for the non-linear-core-corrections (NLCC) for the pseudopotentials used in the present work can be found in Table I. The one-electron Kohn Sham eigenvalues were expanded in a basis of strictly localized 26 numerical atomic orbitals. 20,27 The size of the basis set chosen was doubleζ plus polarization for the valence states of all the atoms. All the parameters that define the shape and the range of the basis functions were obtained by a variational optimization of the energy in the α-cristobalite polymorph of SiO 2 , and of the enthalpy (with a pressure P = 0.2 GPa) in the α-phase of TeO 2 , following the recipes given in Refs. 28  The electronic density, Hartree, and exchange correlation potentials, as well as the corresponding matrix elements between the basis orbitals, were calculated in a uniform real space grid. An equivalent plane wave cutoff of 400 Ry was used to represent the charge density. During the geometry optimizations, we used a 6 × 6 × 6 Monkhorst-Pack mesh 32 for all the Brillouin zone integrations. The macroscopic polarization and its derivatives with respect to an external electric field, depend highly on the number of k-points used. To quantify this dependence, we have refined the Monkhorst-Pack meshes and followed the evolution of the field induced polarization with increasing number of k-points. Further details will be given in Sec. IV B.
For the structural characterization in the absence of an external electric field, the atoms were allowed to relax until the maximum component of the force on any atom was smaller than 0.01 eV/Å, and the maximum component of the stress tensor was smaller than 0.0001 eV/Å 3 . α-SiO 2 cristobalite crystallizes in the same space group (P 4 1 2 1 2, D 4 4 , no. 92) and with the same independent atomic positions as paratellurite, α-TeO 2 . Two atoms are independent by symmetry: one Si atom at position (u,u,0), and one O atom at position (x,y,z). The unit cell contains four formula units (twelve atoms). In the case of α-cristobalite, the SiO 4 entities are almost regular tetrahedra (see Fig. 1) and the Si-O distances are all close to 1.60Å (see Table-II). As shown in Fig. 1, the tetrahedra are organized as to form an helical chain along the z-direction. While in the directions x and y, the structure is made by zig-zag chains. Theoretical lattice parameters and Wickoff positions are reported in Table-II, together with some experimental values for comparison. Although the data summarized in Table-II include results obtained with different implementations of the density functional theory (differences in the electrons included explicitly in the calcula-tion, with and without pseudopotentials, different basis sets, different ways of sampling the Brillouin zone), a general trend that can be observed is that both PBE-GGA and B3LYP hybrid functional yield an overestimation of the lattice constant of up to 3 % with respect to the experimental values. The Wyckoff positions obtained with the different DFT-methods are in very good agreement between them [the maximum difference being in the O(y) position], and are perfectly comparable with the experimental results. Indeed, this good performance in a traditionally complicated system like SiO 2 (very sensible to many approximations, in particular to basis sets 28 ) validates the Siesta basis sets and pseudopotentials used in the present work.
2. α-TeO2 phase α-TeO 2 crystallizes in a tetragonal unit cell with the space group symmetry P 4 1 2 1 2 (D 4 4 , no. 92) as discussed before. 30 An schematic view of the unit cell is depicted in Fig. 2. The tellurium atom occupies the center of triangle bipyramids whose basis is formed by two oxygen atoms and by the tellurium electronic lone pair, and whose apexes are also oxygen atoms. Therefore, the Te atoms is coordinated with four O atoms. This TeO 4 bypiramidal unit, building block of the tellurium oxides discussed in the present work, is referred to as a disphenod. Two different Te-O bond can be distinguished within the disphenod, being the equatorial O atoms closer to Te than the apical O atoms (experimental distances of 1.87Å and 2.12Å, respectively). As in α-SiO 2 cristobalite, the polyhedra are connected by vertices to form a three dimensional network, ressembling to an helical chain along z direction and zig-zag chains along x and y directions (see Fig. 2

.)
Unit cell lattice parameters and internal coordinates are reported in Table-III. The structural parameters obtained with Siesta are in very good agreement with those obtained with a plane wave code with ultrasoft pseudopotentials and an energy cutoff of 30 Ry, showing the good performance of the basis set used in the present work. As usual, the standard overstimation of the experimental equilibrium volume by the generalized gradient approximation is found. The calculated independent bondlengths at the theoretical equilibrium lattice parameters also overstimates the experimental numbers. Nevertheless, the difference between the short equatorial and the long axial Te-O bond lengths is preserved within Siesta.
3. γ-TeO2 phase γ-TeO 2 crystallizes in an orthorhombic unit cell with the space group P 2 1 2 1 2 1 (D 4 2 , no. 19). 39 A schematic view of the unit cell is represented in Fig. 3. This phase is metastable at normal conditions, and has been recently identified by x-ray powder diffraction of recrystallized amorphous TeO 2 doped with oxides. The unit cell contains four formula units (twelve atoms), with three atoms independent by symmetry: one Te atom located at (u,v,w), and two oxygen atoms labeled as O I and O II . As in the α-phase, the structure of the γ-phase can be considered as a 3D network of corner-sharing TeO 4 disphenods. However, in the γ-phase the disphenods are strongly deformed, so the length of the four Te-O bonds are rather different (experimentally the bond lengths range from 1.86Å to 2.20Å), with a much larger spread (0.34Å) than in the α-phase (0.24Å). Indeed, if we assume that the longest Te-O distance (marked with a dashed line in  In Fig. 3, despite the fact that the bridge between Te atoms along the z-oriented chains is always through a O II atom, the length of the Te-O II bond is different, ranging between 1.94Å for one of the bonds of the chain to 2.02Å for the second bond. The theoretical lattice parameters and independent positions of the atoms are reported in Table-IV. The good agreement between Siesta and plane waves results confirms and highlights the transferability of our basis set, that was optimized for the α-TeO 2 structure. Again, the PBE-GGA functional overstimates the experimental equilibrium volume, although the deviation in this case (9 %) is slightly larger than usual. This overstimation translates also in a slight overestimation of the spread of the bond lengths in the case of Siesta (0.40Å).
Regarding the first-principles simulations on the structure of the TeO 2 phases, we can summarize that the disphenoidal configuration of the TeO 4 entities are respected both in the α and γ phases. The lone pair sterical effect is thus conserved in our geometry optimization. The good comparison between our structural parameters and the ones obtained with plane waves 38 support the use of the numerical atomic orbital method implemented in Siesta in the present study.

A. Methodology
When an intense light (a powerful laser) goes through an insulator, this medium responds nonlinearly. Its polarization P can be expressed as a Taylor expansion of the electric field E.
where i, j, k and l refer to cartesian directions, and P s i is the zero-field (spontaneous) polarization. The coefficients of the previous expansion, derivatives of the polarization with respect to the electric field of increasing order, are related with the electrical susceptibilities of the material. In particular, the coefficient of the linear term is directly proportional to the linear dielectric susceptibility (a second-order rank tensor), where ε 0 is the permittivity of free space. Unless otherwise is specified, we will use the SI system of units throughout the paper. In the same way, the coefficients of the quadratic and cubic terms are used to define the second-order and third-order nonlinear susceptibilities, Replacing Eqs. (2)-(4) into Eq. (1), then the expansion of the polarization can be written as In general, the polarization depends both on the valence electrons and the ions. In the present work, we deal only with the electronic contribution. The main reason behind this approximation is that the second (SHG) and third (THG) harmonic generation experiments leading to the second and third order susceptibilities are performed at frequencies high-enough to get rid of ionic relaxations but low enough to avoid electronic excitations. 40 This constraint implies that both the frequency of E and its second and third harmonics are lower than the fundamental adsorption gap. Indeed, SHG and THG experiments are pump-probe settings which are sensitive only to the electronic contributions.
For a completely anisotropic crystal, Eq. (5) implies respectively 9, 27 and 81 elements for χ (1) , χ (2) and χ (3) tensors. These numbers can be reduced significantly by considering the specific crystal class. 41 The nonvanishing elements for the symmetry groups of α-TeO 2 and γ-TeO 2 can be found in Table-V and Table-VI, respectively. Thus, combining the symmetry-allowed components of the susceptibility tensors with a judicious choice of the direction of the applied electric field, we can restrict the expansion in Eq. (5) so that only a given component of the susceptibility is present. Then, its value can be obtained by fitting the polarization versus electric field dependence as obtained from first-principles computations of the response of the bulk materials to an homogeneous electric field. 15,19 To better fix this procedure, let us develop the way the χ (2) yxz of γ-TeO 2 is computed. Taking into account the nonvanishing optical susceptibilities of its crystal class, P 2 1 2 1 2 1 , the expansion of the Eq. (5) in this case is reduced to TABLE V. Symmetry-allowed components of the linear dielectric susceptibility, χ (1) , and the second-order, χ (2) , and third-order, χ (3) , nonlinear optical susceptibilities tensors for the P 41212 space group, in which both the α-TeO2 and the α-SiO2 cristoballite crystallize.
χ (1) χ (2) χ (3) xx xyz = −yxz xxxx = yyyy yy xzy = −yzx zzzz zz zxy = −zyx yyzz = xxzz yzzy = xzzx xxyy = yyxx yzyz = xzxz xyxy = yxyx zzyy = zzxx zyyz = zxxz zyzy = zxzx xyyx = yxxy To isolate the χ yxz component, we can apply a field with only x and z components, E = (E, 0, E), so the expansion in Eq. (6) reduces to Assuming the Kleinman symmetry conditions, 42 that states that the nonlinear optical susceptibility tensor, as defined in Eq. (3), is symmetric under a permutation of the i, j, k indices so χ yzx , then Eq. (7) transforms into Up to now we have applied only basic definitions and symmetry properties. The ingredient from firstprinciples comes from the computation of the fieldinduced electronic polarization (the ionic cores are clamped at the equilibrium zero-field positions) as a function of the external electric field. Here, we have used the method of Refs. 15 and 19. In these milestone works, Souza,Íñiguez and Vanderbilt showed how to compute stationary states of a periodic system at a finite, constant electric field E on a uniform discrete k-point mesh, providing that the magnitude of the field does not exceed a critical value E c (N ), that decreases as the number of k-points N increases. The algorithm, described in Sec. V of Ref. 19, is based on the diagonalization of a field-dependent Hermitian operator,T k (E) in the notation of Ref. 19, for every k-point in the mesh. The fielddependent operator depends explicitly on a given k, but implicitly couples neighboring k-points. Therefore, to know the occupied Bloch-like eigenstates, the diagonalization has to be iterated until the procedure converges at all the k-points and the occupied subspace stabilizes. Once self-consistency has been achieved, the overlap matrix between the cell periodic parts of the Bloch functions at neighboring k-points can be used to compute the component of the polarization parallel to a vector of the reciprocal lattice, as usual in the discretized formulation of the polarization of a solid as a Berry-phase. 12 Once a set of polarization versus electric field data have been obtained from first principles, one can make a choice between different methods to extract a value for the nonlinear susceptibility. The first one is a direct fit of the raw data (see the symbols of Fig. 5) to expressions like Eq. (8). The second one is to compute the derivatives of the macroscopic polarization with respect the electric field by finite differences, using the Richardson's extrapolation to estimate the limit h → 0 from calculations with two different step sizes: where D (1) is given by the centered finite difference expresion D (2) is given by and D (3) is given by Here, in the results quoted below, we have used a field step of h = 0.04 V/Å. Comparison between the results obtained with the two methods will be presented in the next Sec. IV B.

B. Convergence studies
It is well known that the total energy ground state calculations for insulators converge expeonentially fast with respect to k-point sampling. However, while a given kpoint sample might be perfectly adequate for some properties, it might constitute too severe approximation for others. 43 In other words, the convergence with respect the k-point sampling is property-dependent. That is the case of the computation of the polarization and its derivatives, when a discretized Berry phase polarization expression is used. 44 In the formalism developed in Ref. 19, and summarized in Sec. IV A, both the calculation of the stationary states of the periodic system at a finite constant electric field, and the Berry-phase polarization is made on a uniform discrete k-point mesh in the first-Brillouin zone, and the convergence for finite field simulations is considerably slower than in total energy or charge density calculations. 44 The situation does not improve significantly when a perturbation approach is used instead of the finite field method. Previous firstprinciples simulations on the nonlinear optical susceptibilities in the framework of the density functional perturbation theory, 40 show that the convergence of the results is quite poor with the number of spetial k points, at least when the discretization of the formula for the Berry phase of the occupied bands is performed after the perturbation expansion (although the situation improves dramatically when the perturbation expansion is performed after the discretization). We can wonder how rapidly our results converge with respect to the k-point sampling.
In Fig. 6 we represent the χ The results are shown for the two different methods used to obtain the derivatives of the macroscopic polarization with respect the electric field: the direct fit to the polarization versus field curve and the Richardson extrapolation to compute finite differences derivatives. The results provided by both methods are in good agreement, with differences smaller than 2 % for reasonable sizes of the Monkhorst-Pack mesh. For the rest of the paper all the reported results have been obtained with the Richardson extrapolation method.
In any case, as we can see from Fig. 6, the convergence with respect the number of k-points included in the simulations is rather slow. A way of predicting converged values 44 for a given magnitude p would be to explapolate to N → ∞ through a least square fit against an analytical formula of the form p = p ∞ + a/N b (where p ∞ , a, and b are adjustable parameters). However, this would require several calculations at different N . Here, we have proceed in a different way, computing only the values of the NLO susceptibilities for N = 20. This would lead to results with errors of usually around 2% for the secondorder and up to 15 % for the third-order susceptibilities. However, a similar trend is expected for estimations of the susceptibilities in different phases with approximately the same unit cell volume and number of atoms per unit cell. Therefore, for the consequence of the main goal of this work, that is, the comparison of the nonlinear optical susceptibilities between phases to ascertain the origin of their large values, the previous approach is justified.

A. Second order NLO susceptibilities
The values for the symmetry allowed components of the second order susceptibilities of α-TeO 2 and γ-TeO 2 are gathered in Table-VII, where they are presented in terms of the d tensor, often used in the literature of nonlinear optics, and defined as where the first index, i, refers to a cartesian direction (1 for x, 2 for y, and 3 for z), while the second index l contracts the two other cartesian indices j and k in Voigt notation (see Table VIII). According to (i) the conditions imposed by the space groups (see Tables V and VI), and (ii) the Kleinman symmetry, which means that χ (2) ijk is symmetric under a permutation of i, j, and k, all the elements of the NLO susceptibility tensor of α-SiO 2 and α-TeO 2 should vanish, and all the independent elements of γ-TeO 2 should be equal. To illustrate this, let us take, for instance, the d 14 element for α-TeO 2 . The crystal symmetry imposes that

Material Element Theory Experiment
while following the Kleinman symmetry The only way to satisfy Eqs. (14) and (15) simultaneously is d 14 = d 25 = 0. A similar reasoning for γ-TeO 2 , where in principle all the crystal symmetry allowed elements of the second order nonlinear susceptibility might be independent, shows that the Kleinman rule forces them to be equal. A permutation of indices impose that Both results are clearly visible in Table VII. Unfortunately, the comparison with the experiment is not straightforward. The only experimental results published up to now on crystals concerns the α-SiO 2 critoballite, 47 and α-TeO 2 phase, 48 resulting in a weak SHG response (indicating the presence of a second order susceptibility) with the SHG efficiency of α-TeO 2 five times larger than in SiO 2 . The metastable nature of the γ-TeO 2 phase makes impossible to grow single crystals big enough to be optically characterized.
It can be surprising at first sight to obtain χ (2) values for the α-TeO 2 where, as explained before, the combination of the space group symmetry and the Kleinman's rules 42 should inhibit the presence of a SHG signal. Different mechanisms have been called to explain this apparent inconsistency. First, Kleinman's relations are expected to breakdown if the second harmonic frequency is close to an absorption band of the material, but this was not the case for α-TeO 2 under the measurement conditions. 45 Second, although the Kleinman's rule is always satisfied in non-dispersive media, it can be broken in dispersive materials, as might be the case for the experimental α-TeO 2 samples. As the present calculations were done in the static finite field limit, no dispersive effect can be taken into account so that the Kleinman's conditions are verified.

B. Third order NLO susceptibilities
A glance to Tables V and VI reveals that there are 21 elements of the third order nonlinear susceptibility allowed by spatial symmetry considerations. A firstprinciples estimation of all of the independent elements would be a rather cumbersome task, that is out of the scope of this work. Instead, we have computed the xxxx, yyyy, and zzzz elements of the third order susceptibilities, as they are informative enough for the structure/properties relationship considerations. For those particular components both the polarization and the electric fields involved are directed along the same cartesian direction [see Eq. (4)]. Results for α-SiO 2 , α-TeO 2 and γ-TeO 2 are summarized in Table-IX. The results obtained with Siesta using the variational implementation of the finite field in periodic system approach are in very good agreement with those obtained with Crystal06, where the field was introduced as a sawtooth potential. However, the last approach requires the use of a supercell to adapt the periodicity of the potential and, therefore, the increase of the computational cost of the simulation.
Unfortunately, the comparison with the experiment is not so straightforward, since there are no measurements yet concerning the crystalline phases of TeO 2 or SiO 2 . Therefore, the values reported in Table-IX can be considered as purely predictive. The only χ (3) available data for tellurium oxides has been measured on the corresponding glass. 1 However, Raman spectroscopy measurements 8 have shown that the structure of γ-TeO 2 is close to the structure of this glass, so it is tempting to think that the order of magnitude of the third-order susceptibility should be the same in both the glass and the crystalline γ-TeO 2 phase.
As can be drawn from Table-IX, the third order susceptibilities of TeO 2 and SiO 2 crystals obtained from our first-principles simulations are of the same order of magnitude as the experimental values for the relevant glasses. The very large χ (3) predicted for the two tellurium oxide polymorphs, on the range of 10 −12 esu, are nearly two order of magnitudes larger than for crystallized silica, with an average χ (3) (TeO 2 ) / χ (3) (SiO 2 ) ratio close to 50, thus extending the experimental observations on glasses to the case of crystalline compounds.
We can wonder now about the origin of the large values for χ (3) in TeO 2 . Despite the fact that α-SiO 2 and α-TeO 2 crystallizes in the same space group, their third order susceptibilities are different by two orders of magnitude. This fact points out that structural arrangement by itself can not be responsible for the remarkable NLO properties of TeO 2 .
A more detailed analysis of the χ (3) tensor elements reveals that the highest value are for the z direction in the case of α-SiO 2 and α-TeO 2 , and for the y direction for the γ-TeO 2 . As stated in Sec. III, those are precisely the crystalline directions where the chains display an helical shape. Moreover, for the two structures that  are structurally similar (α-SiO 2 and α-TeO 2 ) nearly the same ratio χ xxxx of around 1.8 is obtained. Since no lone pair is present in SiO 2 compounds, Te lone pair effect (orientation of the LP with respect to the electric field for example) cannot be responsible for these large variations.
This suggests that high χ (3) values are structurally induced and that helical chains are much more favourable than the zig-zag chains structures shown along the x direction for these materials. Indeed, Mirgorodsky et al. 9 and Soulis and coworkers 10 have shown how there is a strong link between the structure of polymer TeO 2 molecules and its nonlinear optical properties, with the chain-like species developing a drastic increase in their third-order hyperpolarizability with increasing chain length. The Te-O-Te bridges play a dominant role in the polarization properties of the long TeO 2 chains. This was attributed to an exceptionally strong nonlocality of the electronic polarization, that is, assuming that the electric field applied at a given point would induce a dipole moment not only at the very point but in the vicinity of this point (extending up to eight-ten neighbors from the point of perturbation). Now, let us turn our atention to the zig-zag chains of the two TeO 2 compounds (along the x and y direction for the α-TeO 2 phase and along the x and z direction for the γ-TeO 2 phase.) They are all very similar in shape and made of asymmetric single Te-O-Te bridges with the bond length sequence −S − L − S − L− (where S and L stand for short and long bond lengths, respectively). We can define an asymmetry index as AI = 100 × (L − S)/L, whose value amounts to 11 for the chains parallel to x or y in α-TeO 2 , 15 parallel to x in γ-TeO 2 , and 4 parallel to z in γ-TeO 2 . It is interesting to note that the AI values evolve as the inverse of the χ (3) values, suggesting that, for a given chain, the more symmetrical the bridge the higher the third order nonlinear optical susceptibility value.

VI. CONCLUSION
The second and third order nonlinear optical susceptibilities of two bulk crystalline phases of TeO 2 , and α-SiO 2 cristoballite have been computed using the variational approach to compute the response of a periodic solid to an external electric field.
The third order nonlinear susceptibilities are in good agreement with previous more expensive theoretical predictions, were the electric field is introduced by means of a sawtooth potential, and with the experimental results for related glass phases. Indeed we were could reproduce the large values for the χ (3) tensor elements of the tellurium oxides, 50 times larger than in pure silica glasses.
Our results provide some clues to explain the origin for the high hypersusceptibilities and the large variations with respect to the crystalline directions. In particular these properties could be attributed to the presence of helical chains in the structure.
A next step to be taken in order to through some light into the origin of the large values of the nonlinear optical susceptibilities would require the calculation of the center of the localized Wannier functions, and their variation with the external fields.
Our results demonstrate how first-principles calculations are an efficient tool to estimate nonlinear optical susceptibilities of crystalline solids. These might contribute to fill the gap usually left by experimental measurements due to the difficulty in growing single crystals big enough to be optically characterized.