# Nonlinear optical properties of glassy TeO2: ab initio modeling of $\chi(3)$ and Hyper Raman spectra

Evgenii M Roginskii, Raghvender Raghvender, Olivier Noguera, Philippe Thomas, Olivier Masson, Assil Bouzid

► To cite this version:

Evgenii M Roginskii, Raghvender Raghvender, Olivier Noguera, Philippe Thomas, Olivier Masson, et al.. Nonlinear optical properties of glassy TeO2: ab initio modeling of $\chi(3)$ and Hyper Raman spectra. Physica B: Condensed Matter, 2025, 699, pp.416768. 10.1016/j.physb.2024.416768. hal-04963214

HAL Id: hal-04963214

https://cnrs.hal.science/hal-04963214v1

Submitted on 27 Feb 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Nonlinear optical properties of glassy TeO₂: *ab initio*
modeling of χ⁽³⁾ and Hyper Raman spectra

Evgenii M. Roginskii*^a$, Raghvender Raghvender^b$, Olivier Noguera^b$, Philippe
Thomas^b$, Olivier Masson^b$, Assil Bouzid*^b

<sup>a</sup>Ioffe Institute, Polyteknicheskaya 26, 194021, St. Petersburg, Russia

<sup>b</sup>Institut de Recherche sur les Céramiques, UMR 7315 CNRS - Université de Limoges,
Centre Européen de la Céramique, 12 rue Atlantis, 87068, Limoges Cedex, France

Abstract

The nonlinear optical properties of amorphous tellurium oxide is studied in a
framework of computational technique based on the *ab initio* quantum-chemical
calculations. The calculated values of refractive index and third-order nonlinear
susceptibility are in a very good agreement with experimental data. The origin
of the high nonlinear properties of TeO₂ is attributed, in general, to the hy-
perpolarizability of Te atoms electron lone pairs and to the hyperpolarizability
of electronic density localized on bonds of the asymmetric Te-O-Te bridges in
part. The computational technique is extended on Hyper Raman spectroscopy
including scattering on longitudinal modes. The spectra was calculated for both
crystalline and glassy TeO₂ with a very good reproduction of experimental data.
An extra bands inactive in classical Raman spectroscopy is studied using this
approach. The most intense high-frequency band in the spectra of paratellurite
is attributed to anti-symmetric stretching of long-short bonds in TeO₄ disphe-
noid with A₂ symmetry. The influence of the structural unit configuration on
the band position in Hyper Raman spectra is established. The best agreement
between calculated phonon spectrum and experimental data is obtained using
hybrid functional approach, which explains strong correlation on valence elec-
tron states of Te atoms. The strong dependence of LO modes intensity in Hyper
Raman spectra on third order nonlinear susceptibility is obtained and explained.
The reported technique can be easily extended on a big family of binary and
ternary compounds.

Keywords: Raman, nonlinear optics, DFT

1. Introduction

Tellurium dioxide-based glasses feature linear and nonlinear optical prop-
erties (NLOP) that make them very promising materials for applications in

Email addresses: e.roinskii@mail.ioffe.ru (Evgenii M. Roginskii*),
assil.bouzid@cnrs.fr (Assil Bouzid*)
telecommunication technologies. In particular, the third-order non linear susceptibilities ($\chi^3$) of TeO$_2$ is the highest among the known oxide glasses with values 50 to 100 times higher than that of fused silica glass [1]. Numerical investigations of the NLOP in TeO$_2$ based materials are limited to the study of crystalline phases or to small isolated clusters. Nevertheless, the cluster approach suffers from the bond breaking issue for chain like structures which leads to significant influence of surface on the simulated optical properties [2]. The crystalline phase calculations of nonlinear properties based on Berry phase approach [3] hardly expands on supercell approach that is necessary when simulating disordered systems, due to the high computational cost in terms of central processor (CPU) requirement and/or memory requirement of large size plane wave basis sets as discussed in previous research [4]. The third order nonlinear susceptibility can be simulated also by calculations of fourth order total energy derivative within Coupled Perturbed Kohn-Sham approach [5] which requires as well a high computational resources even for medium size crystalline systems and in practice implemented only in case of linear combination of atomic orbitals basis set technique. While instructive, these studies remain limited and do not provide an overall understanding of the origins of the NLOP of glassy materials. From a computational point of view, computing NLOP on glassy models faces two main issues: first, reliable glassy models should be produced, second, a proper computational technique suited for periodic and disordered systems is still missing. The present paper tackles both issues. The technique reported here takes advantage of glassy TeO$_2$ models which were recently produced by resorting to *ab initio* molecular dynamics (MD) [6], and whose Raman spectra where calculated reproducing experimental data in excellent agreement allowing thereby an interpretation of the origins of the Boson peak in this material [7]. The glassy models were used as a model objects for the calculation of NLO properties and the obtained results are compared to their experimental counterparts. The calculations of the NLOP on glassy oxides revealed size effects, which are associated with the size of simulation cell. Additionally, formalism on calculations of the Hyper-Raman (HR) spectra is explored on both crystalline and amorphous TeO$_2$. The low frequency band in HR spectra is observed in glassy structure and explained. The methods reported in this paper can be adopted to other binary and ternary glassy systems including systems with modifiers which can be used to achieve a better analysis of the experimental data or to predict the optical properties of amorphous materials in a way to guide their experimental design.

## 2. Computational details

### 2.1. DFT setup

*Ab initio* calculations of optical and dynamical properties were carried out within density functional theory (DFT) as implemented in CP2K software package [8] using the generalized gradient approximation with Perdew-Burke-Ernzerhof (PBE) functional [9] and pseudopotential approximation with core-valence interactions taken into account according to Goedecker-Teter-Hutter scheme [10].
The Gaussian Plane Waves method [11] was used to form wavefunctions basis set with cutoff energy $E_c = 500$ Ha for plane wave part and using triple-zeta valence basis for local combination molecular orbitals part. Geometry relaxation was performed until the forces on atoms are below $10^{-5}$ Ha/Bohr. Unless explicitly stated, the computational setup is identical to that used for computing Raman spectra in Ref. [7].

# 3. Results and Discussion

## 3.1. Nonlinear optical properties

Nonlinear optics is a branch of physics that describes the nonlinear dependence of the polarizability $P$ with respect to an applied electric field $E$. According to the perturbation theory [12] the total energy of the crystalline system is decomposed as a series of derivatives with respect to the electric field strength $E$ such as:

$$ E(\mathcal{E}) = E(0) - \sum_t \mu_t \mathcal{E}_t - \frac{1}{2!} \sum_{tu} \alpha_{tu} \mathcal{E}_t \mathcal{E}_u - \frac{1}{3!} \sum_{tuv} \beta_{tuv} \mathcal{E}_t \mathcal{E}_u \mathcal{E}_v - \frac{1}{4!} \sum_{tuvw} \gamma_{tuvw} \mathcal{E}_t \mathcal{E}_u \mathcal{E}_v \mathcal{E}_w, \quad (1) $$

where, $\mu_t = -\frac{\partial E}{\partial \mathcal{E}_t}$, $\alpha_{tu} = -\frac{\partial^2 E}{\partial \mathcal{E}_t \partial \mathcal{E}_u}$, $\beta_{tuv} = -\frac{\partial^3 E}{\partial \mathcal{E}_t \partial \mathcal{E}_u \partial \mathcal{E}_v}$, $\gamma_{tuvw} = -\frac{\partial^4 E}{\partial \mathcal{E}_t \partial \mathcal{E}_u \partial \mathcal{E}_v \partial \mathcal{E}_w}$, are the dipole moment, the linear polarizability, and the first and second order hyperpolarizabilities, respectively. The derivative of Eq. (1) with respect to $\mathcal{E}$ gives the dipole moment of the system as follows:

$$ \mu_i = \mu_i^0 - \sum_j \alpha_{ij} \mathcal{E}_j - \frac{1}{2} \sum_{jk} \beta_{ijk} \mathcal{E}_j \mathcal{E}_k - \frac{1}{6} \sum_{jkl} \gamma_{ijkl} \mathcal{E}_j \mathcal{E}_k \mathcal{E}_l, \quad (2) $$

where, $\mu^0$ is the dipole moment of the system at zero electric field. Differentiating Eq. (2) with respect to the electric field yields:

$$ \alpha_{ij} = \alpha_{ij}^0 - \sum_k \beta_{ijk} \mathcal{E}_k - \frac{1}{2} \sum_{kl} \gamma_{ijkl} \mathcal{E}_k \mathcal{E}_l, \quad (3) $$

where, $\alpha_{ij}^0$ is the spontaneous polarisibility, which is linked to the first-order nonlinear susceptibility via the following relation: $\chi^{(1)} = \frac{4\pi}{\Omega_0} \bar{\alpha}$. In this relation, $\bar{\alpha}$ is the isotropic average of the polarisibility at zero electric field, and $\Omega_0$ is the volume of the unitcell. The hyperpolarizability $\beta_{ijk}$ is related to the second-order nonlinear susceptibility through $\chi_{ijk}^{(2)} = \frac{2\pi}{\Omega_0} \beta_{ijk}$ and is equal to zero for centrosymmetric systems. Finally, the second order hyperpolarizability $\gamma_{ijkl}$ is linked to the third order-order nonlinear susceptibility $\chi^{(3)}$ as follows: $\chi_{ijkl}^{(3)} = \frac{2\pi}{3\Omega_0} \gamma_{ijkl}$. In practice, $\chi^{(1)}$ and $\chi^{(3)}$ can be easily obtained by fitting the polarizability dependence on the electric field using a second order polynomial expansion. To this end, one have to calculate for various applied electric field
intensities and polarizations their corresponding polarizability $\alpha_{ij}$ values by resorting to linear response theory within density functional formalism.

In this work, an electric field [13] with intensities ranging from 0 to $E_{max} = \frac{E_g}{e \times max(a,b,c)}$ ($E_g$ is the electronic bandgap, a, b, and c are the lattice parameters) with a step of 0.00008 a.u. was applied to the studied systems. For each $\mathcal{E}$ value, three directions of electric field are considered: (1,0,0), (0,1,0) and (0,0,1). Hence, the obtained data provides a full description of the dependence of the polarizability as a function of the electric field. Consequently, fitting the polarizability components $\alpha_{xx}$, $\alpha_{yy}$ and $\alpha_{zz}$ with a second order polynomial expansion according to Eq. (3), give us a direct access to the values of hyperpolarisibilities.

The calculation of the NLO properties were performed on glassy TeO$_2$. The atomistic models produced in Ref. [6] where used where first-principles molecular dynamics has been used to produce accurate atomic scale models of the TeO$_2$ glass yielding a close agreement with experimental results. In particular, two amorphous models with 270 and 480 atoms [6] are considered. For each model 10 representative snapshots extracted from the T = 300 K trajectory of TeO$_2$ obtained from Ref. [6] were treated in order to get the optical properties along MD trajectory.

In order to study the dependence of the calculated $\chi^{(3)}$ on the model size, a 2 × 2 × 2 supercell was constructed by duplicating the 270 atoms model. With the aim of ensuring the consistency of calculations, all the snapshots of the three 270, 480 and 2160 atoms models were fully optimized by allowing the atomic positions and the cell dimensions to relax at fixed angles (due to lack of the symmetry in glass) before computing the polarizability tensor. It is noteworthy, that such technique was successfully applied in calculation of Raman spectra of amorphous systems [7]. It is noteworthy that although the 2160 atoms model was not produced by a full equilibration of the supercell at finite temperature, it is expected to give a reasonable estimation of the polarizability tensor as this later have shown a little dependence on the statistical fluctuations of the glassy structure. In addition to TeO$_2$, models with 270 and 456 atoms amorphous SiO$_2$ previously obtained in Ref. [7] were also considered here and a similar approach of building a 2160 atoms amorphous model was applied.

The results of the NLOP calculations for amorphous TeO$_2$ and SiO$_2$ are presented in Table 1 and compared to experimental data [14, 15, 16]. In the case of amorphous system the mean values of $\chi^{(1)}$, diagonal $\chi_d^{(3)}$ and non-diagonal $\chi_d^{(3)}$ components are obtained as arithmetical mean, while the mean value of nonlinear third order susceptibility is obtained as follows:

$$ \langle \chi^{(3)} \rangle = \frac{1}{5} \left( \chi_{xxxx}^{(3)} + \chi_{yyyy}^{(3)} + \chi_{zzzz}^{(3)} + 2\chi_{xxyy}^{(3)} + 2\chi_{xxzz}^{(3)} + 2\chi_{yyzz}^{(3)} \right) \quad (4) $$

Despite a slight overestimation, the calculations yield values of refractive index $n_0$ for TeO$_2$ and SiO$_2$ in good agreement with experiments indicating a good estimation of $\chi^{(1)}$. Coming to $\chi^{(3)}$, while showing the same experimental order of magnitude, the absolute values calculated on the glassy TeO$_2$ and SiO$_2$
<table><caption>Table 1: Linear refractive index, $n_0$, diagonal ($\chi_d^{(3)}$), non diagonal ($\chi_{nd}^{(3)}$) and average $\langle\chi^{(3)}\rangle$ third order nonlinear susceptibility $\chi^{(3)}$ in $10^{-22} \text{ m}^2/\text{V}^2$ as a function of number of atoms N in simulated cell calculated for glassy TeO$_2$ and glassy SiO$_2$ and their ratio. Experimental reference values are taken from Refs. [14, 15]. In the case of TeO$_2$, the experimental results correspond to $\alpha$-TeO$_2$.</caption><thead><tr><th rowspan="2">Prop. N</th><th colspan="4">SiO<sub>2</sub></th><th colspan="4">TeO<sub>2</sub></th><th colspan="4">Ratio TeO<sub>2</sub>/SiO<sub>2</sub></th></tr><tr><th>270</th><th>456</th><th>2160</th><th>Expt. [15]</th><th>270</th><th>480</th><th>2160</th><th>Expt. [14]</th><th>270</th><th>480</th><th>2160</th><th>Expt. [14]</th></tr></thead><tbody><tr><td>n<sub>0</sub></td><td>1.49</td><td>1.5</td><td>1.51</td><td>1.458</td><td>2.25</td><td>2.34</td><td>2.34</td><td>2.22</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>&chi;<sub>d</sub><sup>(3)</sup></td><td>1.02</td><td>1.22</td><td>1.69</td><td>—</td><td>38.01</td><td>56.05</td><td>84.58</td><td>95.1</td><td>37.26</td><td>45.94</td><td>50.05</td><td>45.4</td></tr><tr><td>&chi;<sub>nd</sub><sup>(3)</sup></td><td>0.41</td><td>0.48</td><td>0.6</td><td>—</td><td>17.31</td><td>24.26</td><td>33.29</td><td>42.0</td><td>16.97</td><td>19.88</td><td>19.7</td><td>20.04</td></tr><tr><td>&chi;<sub>d</sub><sup>(3)</sup>/&chi;<sub>nd</sub><sup>(3)</sup></td><td>2.49</td><td>2.54</td><td>2.82</td><td>—</td><td>2.20</td><td>2.31</td><td>2.54</td><td>2.26</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>&langle;&chi;<sup>(3)</sup>&rangle;</td><td>1.10</td><td>1.31</td><td>1.73</td><td>2.095</td><td>41.49</td><td>62.67</td><td>90.69</td><td>102.8</td><td>37.72</td><td>47.84</td><td>52.42</td><td>49.1</td></tr></tbody></table>

models slightly underestimate their experimental counterpart. This behavior might be ascribed to the underestimation of the bandgap when GGA exchange and correlation functionals are used [4]. Furthermore, it is noteworthy that experimental data obtained for $\alpha$-TeO$_2$ are only obtained for components in the (a,b) plane, therefore, due to experimentally undefined components of nonlinear susceptibility tensor which include z-component, the experimental mean value of $\chi^{(3)}$ in the crystal is slightly underestimated. It is noteworthy, that regardless of the studied system the absolute values of the NLOP show a strong size effect where they increase with increasing size of the simulation cell.

Experimentally, the nonlinear susceptibility of a given material is usually measured with respect to another reference material. A comparative approach is also used in a computational research, where the same calculation parameters and the size of the cell are the same for both materials [17]. In the case of paratellurite single crystal, the measurements were performed with respect to a fused silica reference glass plate [14]. Hence, by focusing on the ratio of the $\chi^{(3)}$ obtained from TeO$_2$ and SiO$_2$ at similar model sizes one grants a legitimacy to the computed non linear susceptibility directly comparable with their experimental counterpart. Table 1 shows the calculated components of $\chi^{(3)}$ on TeO$_2$ versus amorphous SiO$_2$ in comparison to experimental results from Ref. [14].

It is noteworthy, that other computational technique based on finite field perturbation method applied to paratellurite [17] derives significantly overestimated values of nonlinear susceptibility with values equal to $\langle \chi^{(3)} \rangle = 320 \cdot 10^{-22} \text{ m}^2/\text{V}^2$. Contrarily, the molecular approach based on calculations of finite size clusters [2] underestimate the values of third order nonlinear susceptibility ($\langle \chi^{(3)} \rangle \sim 40 \cdot 10^{-22} \text{ m}^2/\text{V}^2$). Interestingly, the results of the calculations reported in this paper yields a third-order nonlinear susceptibility in much better agreement with experimental data than what is reported in the literature. In particular, the $\chi^{(3)}$ ratio shows a minimal size effect for glassy models larger than 480 atoms and converges towards that experimentally measured on $\alpha$-TeO$_2$. Therefore, one can infer that the NLO response of the
material only weakly depends on the structural disorder. This result can be explained by the origins of the high NLO properties of TeO₂ that were associated to the hyperpolarisibility of Te atom electron lone pair and to the asymmetric Te-O-Te bridges. Both features are maintained in the glassy phase [6, 18, 19] and could explain the maintenance of the strong NLO properties of glassy TeO₂.

## 3.2. Hyper-Raman spectra calculations

### 3.2.1. Theoretical aspects

In addition to NLO properties, other nonlinear processes involving phonon-phonon scattering, such as Raman scattering and Infrared absorption, can provide detailed information on the structure of the studied material thanks to the phonons sensitivity to the local structure deviations from the ground state. In practice, when the symmetry of a crystal is high, the phonon spectra can be decomposed into irreducible representations (IRreps) to devise selection rules for the analysis of Raman and Infrared spectra. For example, one can consider the structure of α-TeO₂ which belongs to the tetragonal symmetry with space group P4₁2₁₂ (No. 92) in which Te atoms occupy 4a Wyckoff positions (WP), local symmetry 2 and oxygen atoms occupy 8b WP. The optical phonon spectrum of α-TeO₂ at Γ point can be decomposed as follows:

$$ \Gamma_{\text{optic}} = 4A_1 \oplus 4A_2 \oplus 5B_1 \oplus 4B_2 \oplus 8E \quad (5) $$

According to table 2, the A₂ and E modes are active in Infrared spectrum as they show a linear dependence on the dipole moment which is responsible for the light absorption process. Another way of exploiting the lattice vibrations is the investigation of the Raman spectra. Raman scattering is a nonlinear process based on light scattering on phonon modes being characterised by a 2nd rank Raman tensor [20]:

$$ \mathcal{A}_{ij}^{m} = \sqrt{\Omega_0} \sum_{n,\beta} \frac{\partial \chi_{ij}^{(1)}}{\partial \tau_{n\beta}} u_m(n, \beta), \quad (6) $$

where, $u_m(n, \beta)$ is the displacement of the nᵗʰ atom in the mᵗʰ vibrational mode along cartesian direction $\beta$. Consequently, the selection rules of the Raman scattering process are determined by a quadratic function transformation of the IRreps as shown in Table 2. From this table, one finds that A₂ modes which are active in the Infrared absorption spectrum are “silent” in Raman spectra. Here, “silent modes” refer to vibrational modes that are not observed in Raman scattering under given experimental conditions. This phenomenon arises due to the symmetry properties of the vibrational modes and the selection rules governing Raman scattering.

Looking at other crystalline systems, such as rutile TiO₂, one can identify vibrational modes, such as A₂g modes, that are “silent” in both infrared absorption and Raman spectra [21].

Nevertheless, despite being silent, these modes can provide valuable information about the structure and the symmetry of the studied system. For instance,
<table><caption>Table 2: Character table for point group 422</caption><thead><tr><th>IRreps</th><th>E</th><th>2C<sub>4</sub>(z)</th><th>C<sub>2</sub>(z)</th><th>2C'<sub>2</sub></th><th>2C''<sub>2</sub></th><th>Functions</th></tr></thead><tbody><tr><td>A<sub>1</sub></td><td>+1</td><td>+1</td><td>+1</td><td>+1</td><td>+1</td><td>x<sup>2</sup> + y<sup>2</sup>, z<sup>2</sup></td></tr><tr><td>A<sub>2</sub></td><td>+1</td><td>+1</td><td>+1</td><td>-1</td><td>-1</td><td>z, Rz, z<sup>3</sup>, z(x<sup>2</sup> + y<sup>2</sup>)</td></tr><tr><td>B<sub>1</sub></td><td>+1</td><td>-1</td><td>+1</td><td>+1</td><td>-1</td><td>x<sup>2</sup> - y<sup>2</sup>, xyz</td></tr><tr><td>B<sub>2</sub></td><td>+1</td><td>-1</td><td>+1</td><td>-1</td><td>+1</td><td>xy, z(x<sup>2</sup> - y<sup>2</sup>)</td></tr><tr><td>E</td><td>+2</td><td>0</td><td>-2</td><td>0</td><td>0</td><td>(x, y), (Rx, Ry), (xz, yz), (xz<sup>2</sup>, yz<sup>2</sup>), (xy<sup>2</sup>, x<sup>2</sup>y), (x<sup>3</sup>, y<sup>3</sup>)</td></tr></tbody></table>

the absence of Raman scattering from certain vibrational modes can be indicative of the presence of specific symmetries within the molecule. Therefore, analyzing silent modes can aid in the determination of crystalline structure and can assist one to gain insights into the composition, phase, and properties of materials.

Complementary to Raman scattering, Hyper Raman (HR) scattering, which is a higher-order nonlinear optical process that involves the simultaneous absorption of two photons by a molecule, resulting in one scattered photon with a frequency equal to the sum or difference of the absorbed photon frequencies, is an interesting scheme to investigate silent modes, which are inactive both in Raman and IR spectra. HR process is characterized by the third-rank Hyper Raman tensor featuring symmetry based selection rules that are different from those of Raman scattering [22]:

$$B_{ijk}^{m} = \sqrt{\Omega_0} \sum_{n,\beta} \frac{\partial \chi_{ijk}^{(2)}}{\partial \tau_{n\beta}} u_m(n\beta), \quad (7)$$

For longitudinal optical (LO) phonon modes, the second order susceptibility derivative includes an additional term which accounts for phonon induced long range electric field [20, 23]:

$$\frac{d\chi_{ijk}^{(2)}}{d\tau_{n\lambda}} = \frac{\partial\chi_{ijk}^{(2)}}{\partial\tau_{n\lambda}}|_{\epsilon=0} - \frac{12\pi}{\Omega_0} \frac{\sum_l Z_{n\lambda l}^* q_l}{\sum_{l,l'} q_l \epsilon_{ll'} q_{l'}} \sum_l \chi_{ijkl}^{(3)} q_l, \quad (8)$$

where, $\epsilon_{ll'}$ is the macroscopic high frequency dielectric tensor and $Z_{n\lambda l}^*$ represents the Born effective charge of n$^{th}$ atom defined as the derivative of the polarization with respect to atomic displacement: $Z_{ij}^* = \Omega_0 \frac{\partial P_i}{\partial \tau_j}$. The values of $Z_{ij}^*$ tensor can be obtained using finite difference method by calculating dipole moments for each atomic displacement. In the case of amorphous materials, as there is no preferential crystal orientation, Hyper Raman spectra should be simulated in assumption of powder sample, i.e. unoriented one, or as an average spectrum integrated over all possible polarizations. This can be done by considering rotational invariants as follows [24]:
$$ \begin{aligned} \langle \beta^2 \rangle = & \frac{19}{105} \sum_i \beta_{iii}^2 + \frac{13}{105} \sum_{i \neq j} \beta_{iii} \beta_{ijj} + \\ & \frac{44}{105} \sum_j \beta_{ijj}^2 + \frac{13}{105} \sum_{ijk, \text{cyclic}} \beta_{iij} \beta_{jkk} + \frac{5}{7} \beta_{123}^2 \end{aligned} \quad (9) $$

In practice, the calculations of Hyper Raman spectra are performed in two steps: First, phonon modes are obtained from the dynamical matrix. Second, the susceptibility polarization are derived with respect to atomic displacements. The Hyper Raman activities is then obtained by applying Eq. (7).

### 3.2.2. The case of paratellurite α-TeO₂

In order to check the validity of this approach, the calculations of the HR spectra for paratellurite α-TeO₂ were performed by resorting to PBE and PBE0 exchange and correlation functionals.

The calculated values of the optical phonon modes are provided in Table 3 and compared to their experimental counterpart. By looking at the PBE results, it appears that while low-frequency modes are relatively well described, the high-frequency range exhibits a significant blue shift with respect to experimental data up to 100 cm⁻¹ leading to an overall root mean square error (RMSE) of 67 cm⁻¹ compared to experiments. By accounting for a fraction of the exact Hartree-Fock exchange (PBE0 with α=1/2 or α=1/3), it is notable that the RMSE is drastically reduced down to 19.4 cm⁻¹ for α=1/3 calculations. Interestingly, an overall improvement in both low- and high-frequency ranges is observed. This behavior is ascribed to the improved description of the Te-O bond lengths in hybrid exchange and correlation functionals.

According to the above analysis (see Table 3), the experimental band positions in HR spectra (which are related to phonon frequencies) are expected to be well reproduced when hybrid exchange and correlation functional is used with the best agreement found in the case of PBE0 (α=1/3). Nevertheless, the value of LO-TO splitting, which demonstrate the accuracy of dielectric properties calculations, is nearly the same when comparing the PBE and PBE0 results. Therefore, in order to ensure a balance between accuracy and computational cost, the HR spectra are calculated by relying on the band positions (i.e. eigen values and eigen vectors $u_m(n, β)$ used in Eq.7) obtained from calculations using the PBE0(α=1/3) approximation and derivatives of dielectric susceptibility ($\frac{\partial\chi}{\partial\tau}$) obtained from PBE calculations. This scheme is referred to hereafter as PBE0freq-PBEdielectric. Following this scheme, the calculated HR spectra of α-TeO₂ is compared to experimental data in Figure 1 for a set of the most significant polarizations.

According to the selection rules, the HR spectra in XZZ and XYY polarizations should reveal modes with E(TO+LO) symmetry. From Fig.1, it is quite clear that the calculated spectra reproduce most of the experimental peak positions up to a frequency shift while showing discrepancies on relative intensities. The low frequency band in XYY polarization is mostly related to E(TO) mode,
<table><caption>Table 3: Frequencies of optical phonon modes of α-TeO₂, LO modes are given in parentheses</caption><thead><tr><th>IRreps</th><th>PBE</th><th>PBE0 (α=1/2)</th><th>PBE0 (α=1/3)</th><th>Experimental</th></tr></thead><tbody><tr><td>B1</td><td>78.84</td><td>48.36</td><td>72.84</td><td>62</td></tr><tr><td>A2</td><td>99.77 (126.60)</td><td>99.78 (131.61)</td><td>105.81 (133.31)</td><td>82</td></tr><tr><td>E</td><td>143.19 (144.65)</td><td>130.36 (133.43)</td><td>149.07 (151.44)</td><td>122 (122)</td></tr><tr><td>B1</td><td>132.47</td><td>138.98</td><td>148.76</td><td>125</td></tr><tr><td>A1</td><td>134.98</td><td>146.80</td><td>155.33</td><td>148</td></tr><tr><td>B2</td><td>131.67</td><td>150.20</td><td>155.69</td><td>150</td></tr><tr><td>B1</td><td>156.36</td><td>172.65</td><td>182.95</td><td>177</td></tr><tr><td>E</td><td>172.73 (185.49)</td><td>182.71 (202.95)</td><td>193.84 (209.20)</td><td>173 (197)</td></tr><tr><td>E</td><td>191.41 (230.90)</td><td>212.28 (250.47)</td><td>220.44 (259.39)</td><td>207 (239)</td></tr><tr><td>A1</td><td>195.72</td><td>221.62</td><td>241.17</td><td>233</td></tr><tr><td>B1</td><td>238.38</td><td>229.67</td><td>242.10</td><td>210</td></tr><tr><td>A2</td><td>255.26 (259.06)</td><td>273.30 (273.52)</td><td>278.71 (278.81)</td><td>262</td></tr><tr><td>B2</td><td>243.82</td><td>275.39</td><td>291.12</td><td>287</td></tr><tr><td>E</td><td>268.06 (301.48)</td><td>293.20 (326.41)</td><td>303.77 (342.68)</td><td>296 (337)</td></tr><tr><td>A2</td><td>289.77 (347.71)</td><td>312.23 (377.87)</td><td>326.58 (393.76)</td><td>321</td></tr><tr><td>E</td><td>308.65 (313.37)</td><td>328.24 (335.49)</td><td>343.62 (355.08)</td><td>336 (?)</td></tr><tr><td>E</td><td>320.46 (360.80)</td><td>338.97 (382.94)</td><td>358.48 (394.82)</td><td>392 (418)</td></tr><tr><td>A1</td><td>372.21</td><td>394.13</td><td>412.69</td><td>393</td></tr><tr><td>B2</td><td>387.57</td><td>413.05</td><td>435.24</td><td>415</td></tr><tr><td>B1</td><td>451.99</td><td>529.97</td><td>555.76</td><td>590</td></tr><tr><td>A2</td><td>457.62 (618.41)</td><td>536.66 (692.67)</td><td>562.19 (711.56)</td><td>592</td></tr><tr><td>E</td><td>526.98 (583.66)</td><td>601.31 (656.54)</td><td>626.24 (679.50)</td><td>645 (720)</td></tr><tr><td>A1</td><td>531.13</td><td>607.34</td><td>631.34</td><td>648</td></tr><tr><td>E</td><td>640.16 (731.38)</td><td>726.25 (823.42)</td><td>751.82 (844.69)</td><td>766 (812)</td></tr><tr><td>B2</td><td>652.18</td><td>742.18</td><td>765.18</td><td>784</td></tr><tr><td>RMSE</td><td>67.2</td><td>27.0</td><td>19.4</td><td></td></tr></tbody></table>

while other high frequency bands are related to LO modes and depend on the value of the third order susceptibility as it can be seen in Eq. 8. Therefore, difference in the calculated $\chi^{(3)}$ values with respect to experimental data leads to differences in TO/LO modes intensity ratio. At the same time the classical Raman scattering intensity of LO modes depends on values of second order susceptibility [20], which is zero in α-TeO₂ due to tetragonal symmetry [26], and therefore the intensity of TO modes in Raman spectra obtain the same amplitude as LO ones. The latter explains the relative high intensity of the high frequency bands in hyper Raman spectrum compared to classical Raman one.

The bands in HR spectra in ZZZ polarization corresponds to modes with symmetry $A_2$ and they are quite well reproduced in the calculated spectrum compared to experimental data, except bands that are marked with asteriks in the experimental spectrum (see Figure 1). These particular bands are forbidden in ZZZ polarization E modes and their experimental manifestation might be due to a slight depolarization effect. Overall, these results show that the scheme of hyper Raman calculation is successful in the case of paratellurite crystal, and can be easily adopted to amorphous system as described in the following section.
<figure><img src="image_2.png" /><figcaption>Figure 1: Hyper Raman spectra of α-TeO₂ in XZZ (a) XYY (b) and ZZZ (c) polarizations computed using the PBE0*freq-PBE*dielectric scheme. Cyan curves show the calculated spectra and the red lines show experimental data from Ref. [25]. The bands marked with asterisks in ZZZ-spectrum (c) denotes forbidden E(LO) modes.</figcaption></figure>

### 3.2.3. The case of amorphous TeO₂

The calculations of Hyper Raman spectra for glassy materials is much more computational demanding due to a large number of atoms in the simulation cell [27]. First, the calculations were performed using PBE exchange and correlation functional and the HR spectra was derived as an average over 5 glassy TeO₂ snapshots [6]. The obtained results, shown in Figure 2, reveal band positions shifted towards to low frequencies, similar to the case of paratellurite crystal discussed above. To improve band positions, the vibrational modes eigenvalues of TeO₂ glass were calculated within the hybrid PBE0(α=1/3) functional while the polarizabilities derivatives are calculated with PBE functional (PBE0*freq-PBE*dielectric scheme) in order to ensure an effective computational cost. Accordingly, the constructed HR spectra is plotted in Figure 2, that obviously reproduces quite well the spectrum obtained within the PBE approximation while it feature a better description of the bands positions at higher frequencies.

The most significant feature in the HR spectra is observed in the low frequency range and corresponds to the Boson peak. The origins of this peak in glassy TeO₂ were previously discussed in the framework of classical Raman scattering [7] where it has been shown that this band is related to vibrations of the structural units joined in nanoclusters as a whole. Therefore, the calculations with PBE0*freq-PBE*dielectric scheme reveal that low frequency HR spectra cor-
<figure><img src="image_3.png" /><figcaption>Figure 2: Calculated Raman and Hyper Raman spectra of amorphous TeO₂ calculated in PBE, and PBE0*freq-PBE*dielectric and approximations.</figcaption></figure>

relates very well with classical Raman scattering in amorphous solids. Coming to the middle frequency range, the HR spectra (Figure 2) reveals a slow decay wing of Boson peak which includes weak bands at c.a. 350 and 450 cm⁻¹ related to anti-symmetric and symmetric stretching-bending modes of TeO₄ disphenoid structural units which corresponds to “silent” A₂ modes in crystalline α-TeO₂ (Figure 2). These results are inline with the analysis of the Hyper Raman spectra of TeO₂-TaO₅/₂-ZnO glassy system in Ref. [28]. Finally, the high frequency range of the spectrum reveals a number of bands which fairly correlates with the spectrum of crystalline α-TeO₂. Namely, the band at 600 cm⁻¹ can be attributed to anti-symmetric stretching of long-short bonds in TeO₄ disphenoid which corresponds to A₂ vibrational mode in crystalline TeO₂. The bands in range 700 – 800 cm⁻¹ are related to a combination of stretching and scissor vibrations in distorted disphenoids and breathing modes in TeO₃ units. The high frequency mode bands are much more pronounced in hyper Raman spectra compared to bands in classical Raman spectra, which can be explained by high third order susceptibility which determine the intensity of the most high frequency LO bands as discussed in the case of paratellurite HR spectra.

4. Conclusions

The nonlinear properties of crystalline and amorphous tellurium oxide were calculated using a combination of linear response theory and finite difference method. The results obtained by reported computational scheme are in a good agreement with experimental data. Moreover, the obtained values of third order nonlinear susceptibility of crystalline TeO₂ significantly improved compared to a previously reported calculations results.

The influence of the cell size on value of nonlinear susceptibility was revealed with convergence for models with 2160 atom in the unitcell. Such size effect can be minimized when calculations are based on reference system like silica. The method reported in this work yields a third-order nonlinear susceptibility of TeO₂ ≈ 50 times higher than that of fused silica, in a very good agreement with experiments. The previously reported technique to calculate Raman spectra [7]
was extended to take into account the second order nonlinear susceptibility thereby enabling the calculation of Hyper Raman spectra of both crystalline and glassy materials. The calculated HR spectra of α-TeO₂ reproduce all the peculiarities of experimental one including modes with A₂ symmetry (forbidden in classical Raman), while the spectra of glass demonstrate the same shape of classical Raman spectrum in the low-frequency range. Therefore the low frequency band known as Boson peak in classical Raman spectra is observed in the Hyper Raman spectroscopy as well. Particular attention in analysis of the spectra is given to LO modes, which in contrast to classical Raman, strongly depends on nonlinear properties of TeO₂ glass in HR spectroscopy. The latter give an extension of the HR spectroscopy as a method to study third order nonlinear properties of the solid state. Finally, the calculated phonon spectrum within hybrid functional method is found as the most accurate technique reproducing experimental frequencies in the whole range, therefore strong correlation on valence electron states of Te atoms takes place. The results give rise to understanding the nature of high nonlinear properties in both crystalline and amorphous tellurium oxides which related mostly to lone pair electron states localized on tellurium atoms.

**Declaration of competing interest**

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Data availability**

Data will be made available on request.

**5. Author Information**

**5.1. ORCID**

Evgenii M. Roginskii: 0000-0001-7511-2910
Raghvender Raghvender: 0000-0001-7511-2910
Olivier Noguera: 0000-0003-1021-1752
Philippe Thomas: 0000-0002-8431-8272
Assil Bouzid: 0000-0002-9363-7240

**5.2. CRediT authorship contribution statement**

**Evgenii M. Roginskii:** Conceptualization, Methodology. **Raghvender Raghvender:** Data curation, Visualization, Software. **Olivier Noguera:** Formal analysis, Validation. **Olivier Masson:** Project administration, Resources, Methodology. **Philippe Thomas:** Supervision, Validation. **Assil Bouzid:** Conceptualization, Writing - original draft, Methodology.
### 5.3. Notes

Authors declare no competing financial interest.

## Acknowledgments

This work was partially supported by the ANR the project AMSES (ANR-20-CE08-0021). Computational resources were provided, in part, by the super-computer facility at the Ioffe Institute and resources provided by the GENCI (project No. AX0910832), and by région Nouvelle Aquitaine via the CANaMIAS project AAPR2021-2020-11779110.

## References

[1] B. Jeansannetas, S. Blanchandin, P. Thomas, P. Marchet, J. Champarnaud-Mesjard, T. Merle-Mejean, B. Frit, V. Nazabal, E. Fargin, G. Le Flem, et al., Glass structure and optical nonlinearities in thallium (i) tellurium (iv) oxide glasses, J. Solid State Chem. 146 (2) (1999) 329-335. doi: 10.1006/jssc.1999.8355.

[2] O. Noguera, M. Smirnov, E. Roginskii, P. Thomas, Microscopic insight into structural and optical properties of glassy TeO₂ derived from finite-cluster ab initio modelling, Physica B: Condensed Matter 652 (2023) 414616. doi: 10.1016/j.physb.2022.414616.

[3] R. Resta, Macroscopic polarization in crystalline dielectrics: the geometric phase approach, Reviews of Modern Physics 66 (3) (1994) 899-915. doi: 10.1103/revmodphys.66.899.

[4] E. M. Roginskii, V. G. Kuznetsov, M. B. Smirnov, O. Noguera, J.-R. Duclère, M. Colas, O. Masson, P. Thomas, Comparative analysis of the electronic structure and nonlinear optical susceptibility of α-TeO₂ and β-TeO₃ crystals, J. Phys. Chem. C 121 (22) (2017) 12365-12374. doi:10.1021/acs.jpcc.7b01819.

[5] M. Ferrero, M. Rérat, B. Kirtman, R. Dovesi, Calculation of first and second static hyperpolarizabilities of one- to three-dimensional periodic compounds. implementation in the CRYSTAL code., The Journal of Chemical Physics 129 (24) (2008) 244110. doi:10.1063/1.3043366.

[6] R. Raghvender, A. Bouzid, S. Cadars, D. Hamani, P. Thomas, O. Masson, Structure of amorphous teo₂ revisited: A hybrid functional ab initio molecular dynamics study, Phys. Rev. B 106 (17) (2022) 174201. doi:10.1103/PhysRevB.106.174201.

[7] E. M. Roginskii, R. Raghvender, O. Noguera, P. Thomas, O. Masson, A. Bouzid, Ab initio study of the Raman spectra of amorphous oxides: Insights into the Boson peak nature in glassy TeO₂, physica status solidi (RRL) Rapid Research Letters 17 (4) (Feb. 2023). doi:10.1002/pssr.202200505.
[8] T. D. Kühne, M. Iannuzzi, M. Del Ben, V. V. Rybkin, P. Seewald, F. Stein, T. Laino, R. Z. Khaliullin, O. Schütt, F. Schiffmann, D. Golze, J. Wilhelm, S. Chulkov, M. H. Bani-Hashemian, V. Weber, U. Bortnik, M. Taillefumier, A. S. Jakobovits, A. Lazzaro, H. Pabst, T. Müller, R. Schade, M. Guidon, S. Andermatt, N. Holmberg, G. K. Schenter, A. Hehn, A. Bussy, F. Belleflamme, G. Tabacchi, A. GlöSS, M. Lass, I. Bethune, C. J. Mundy, C. Plessl, M. Watkins, J. VandeVondele, M. Krack, J. Hutter, CP2K: An electronic structure and molecular dynamics software package - Quick-step: Efficient and accurate electronic structure calculations, The Journal of Chemical Physics 152 (19) (May 2020). doi:10.1063/5.0007045.

[9] J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77 (18) (1996) 3865–3868, erratum Phys. Rev. Lett. 78, 1396 (1997). doi:10.1103/physrevlett.77.3865.

[10] S. Goedecker, M. Teter, J. Hutter, Separable dual-space gaussian pseudopotentials, Phys. Rev. B 54 (3) (1996) 1703. doi:10.1103/PhysRevB.54.1703.

[11] J. VandeVondele, J. Hutter, Gaussian basis sets for accurate calculations on molecular systems in gas and condensed phases, J. Chem. Phys. 127 (11) (2007) 114105. doi:10.1063/1.2770708.

[12] X. Gonze, Adiabatic density-functional perturbation theory, Physical Review A 52 (2) (1995) 1096–1114. doi:10.1103/PhysRevA.52.1096.

[13] P. Umari, A. Pasquarello, Ab initio molecular dynamics in a finite homogeneous electric field, Phys. Rev. Lett. 89 (15) (2002) 157602. doi:10.1103/PhysRevLett.89.157602.

[14] J.-R. Duclère, T. Hayakawa, E. Roginskii, M. Smirnov, A. Mirgorodsky, V. Couderc, O. Masson, M. Colas, O. Noguera, V. Rodriguez, et al., Third order nonlinear optical properties of a paratellurite single crystal, J. Appl. Phys. 123 (18) (2018) 183105. doi:10.1063/1.5020646.

[15] I. H. Malitson, Interspecimen comparison of the refractive index of fused silica, Journal of the Optical Society of America 55 (10) (1965) 1205. doi:10.1364/JOSA.55.001205.

[16] N. Uchida, Optical properties of single-crystal paratellurite (TeO$_{2}$), Physical Review B 4 (10) (1971) 3736. doi:10.1103/physrevb.4.3736.

[17] M. B. Yahia, E. Orhan, A. Beltrán, O. Masson, T. Merle-Méjean, A. Mirgorodski, P. Thomas, Theoretical third-order hyperpolarizability of paratellurite from the finite field perturbation method, J. Phys. Chem. B 112 (35) (2008) 10777–10781. doi:10.1021/jp805050s.

[18] A. Gulenko, O. Masson, A. Berghout, D. Hamani, P. Thomas, Atomic simulations of teo$_{2}$-based glasses: interatomic potentials and molecular dynamics, Phys. Chem. Chem. Phys. 16 (27) (2014) 14150–14160. doi:10.1039/C4CP01273A.
[19] L. Torzuoli, A. Bouzid, P. Thomas, O. Masson, An enhanced core-shell in-teratomic potential for te-o based oxides, Mater. Res. Express 7 (1) (2020) 015202. doi:10.1088/2053-1591/ab6128.

[20] M. Veithen, X. Gonze, P. Ghosez, Nonlinear optical susceptibilities, Raman efficiencies, and electro-optic tensors from first-principles density functional perturbation theory, Physical Review B 71 (2005) 125107. doi:10.1103/PhysRevB.71.125107.

[21] U. Balachandran, N. Eror, Raman spectra of titanium dioxide, Journal of Solid State Chemistry 42 (3) (1982) 276-282. doi:10.1016/0022-4596(82)90006-8.

[22] V. Denisov, B. Mavrin, V. Podobedov, Hyper-Raman scattering by vibrational excitations in crystals, glasses and liquids, Physics Reports 151 (1) (1987) 1-92. doi:10.1016/0370-1573(87)90053-6.

[23] W. D. Johnston, Nonlinear optical coefficients and the Raman scattering efficiency of LO and TO phonons in acentric insulating crystals, Physical Review B 1 (8) (1970) 3494-3503. doi:10.1103/PhysRevB.1.3494.

[24] S. K. Kurtz, T. T. Perry, A powder technique for the evaluation of nonlinear optical materials, J. Appl. Phys. 39 (8) (1968) 3798-3813. doi:10.1063/1.1656857.

[25] V. Rodriguez, M. Couzi, F. Adamietz, M. Dussauze, G. Guery, T. Cardinal, P. Veber, K. Richardson, P. Thomas, Hyper-raman and raman scattering in paratellurite TeO₂, J. Raman Spectrosc. 44 (5) (2013) 739-745. doi: 10.1002/jrs.4251.

[26] R. W. Boyd, Nonlinear Optics, Academic Press Inc, 2007.

[27] P. Umari, A. Pasquarello, Hyper-raman spectrum of vitreous silica from first principles, Phys. Rev. Lett. 98 (17) (2007) 176402. doi:10.1103/PhysRevLett.98.176402.

[28] V. Rodriguez, G. Guery, M. Dussauze, F. Adamietz, T. Cardinal, K. Richardson, Raman gain in tellurite glass: How combination of IR, Raman, Hyper-Raman and Hyper-Rayleigh brings new understandings, J. Phys. Chem. C 120 (40) (2016) 23144-23151. doi:10.1021/acs.jpcc.6b07627.