Cost-effective method for computational prediction of thermal conductivity in optical materials based on cubic oxides

In this paper we report on a computationally cost-effective method designed to estimate the thermal conductivity of optical materials based on cubic oxide including mixed ones, i.e. solid solutions of different oxides. The proposed methodology take advantage from Density Functional Theory (DFT) calculations to extract essential structural parameters and elastic constants which represent the inputs for revised versions of Slack and Klemens equations relating thermal conductivity to elastic constants. Slack equation is modified by the introduction of a corrective factor that incorporates the Grüneisen parameter γ, while in the revised Klemens equation a distortion parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d$$\end{document}d accounting for the impact of point defects on lattice symmetry is added, which is a critical factor in determining thermal conductivity in optical materials with mixed compositions. The theoretical results were found in good agreement with experimental data, showing the reliability of our proposed methodology.

The thermal conductivity of a solid-state system quantifies the heat flux through the propagation of vibrational energy from one atom to adjacent atoms in the lattice without the transport of matter.In essence, it measures a material ability to conduct heat.Thermal transport properties are important for different categories of materials, including thermo-electrics, opto-electronics, photovoltaic and photo-electrochemical cells, batteries.The accurate knowledge of thermal properties is particularly crucial for the design of high-intensity laser systems, where an ideal gain material should exhibit high thermal conductivity to promote efficient heat removal and prevent thermal stresses 1 .The incorporation of lasing dopants ions, such as Yb 3+1-7 or Tm 3+8,9 , into a host material (crystal or ceramics) leads to the formation of localized lattice defects, which can reduce the overall thermal conductivity of the material, eventually affecting its laser emission performance.The selection of a host material with a high intrinsic thermal conductivity is a critical consideration in the development of mixed laser ceramics, as it can mitigate the thermal degradation caused by the presence of these dopants.The material should possess sufficient thermal stability and conductivity even in the presence of high concentrations of dopant, so that it can maintain a high level of thermal conductivity during the demanding conditions of laser operation.
Of course, thermal properties of materials can be experimentally assessed by direct measurements, but this requires the availability of good quality and large samples and the execution of complex and time-consuming measurements.For this reason, the availability of accurate computational methods for the prediction of thermal properties is attracting attention as a convenient alternative to experimental properties screening, in particular when dealing with new, scarcely available and poorly characterized materials and compounds.
Currently, the prediction of thermal conductivity is based on three categories of methods: Anharmonic Lattice Dynamics (ALD) in combination with phonon transport calculation using the Boltzmann transport equation (BTE) and Fourier's law 10 ; Equilibrium Molecular Dynamics (EMD) using the Green-Kubo formula 11 and the direct evaluation of the heat flux by Non-Equilibrium Molecular Dynamics (NEMD) 12,13 .These current state-of-the-art methods are computationally expensive and time-consuming, limiting their practical applicability especially for complex systems.This is due to the need to solve complex equations and perform multiple

Ab initio calculations of lattice structures
The structural calculations were performed using the plane wave periodic DFT implemented in CASTEP 16 employing the PBEsol exchange-correlation functional 17 with Grimme D3 dispersion correction 18 .The ultrasoft pseudopotentials from the internal QC5 library of CASTEP were used for the Y, Tm, Lu, Al, Sc, and O atoms, with a plane wave cut-off of 410 eV, a Self-Consistent Field convergence threshold equal to 10 −8 eV/atom , a k-point grid with a fine k-point separation of 0.04 Å −1 and a convergence criterion for the maximum force component equal to 0.01 eV/Å.The Limited memory Broyden-Fletcher-Goldfarb-Shanno algorithm (LBFGS) 19 , which scales linearly with respect to the system size, was employed to perform the geometric optimization.This optimizer utilizes a limited number of inverse Hessian updates to construct a new Hessian, a process which requires significantly less computational effort than the traditional BFGS algorithm, which scales quadratically with the system size 20 .
Furthermore, Periodic Boundary Conditions (PBCs) were employed to account for the long-range periodicity inherent to crystalline solid materials.
We started our work from conventional unit cells belonging to the cubic space group Ia3 and containing 16 cell formula units (Z = 16) for the study of all the systems of our interest.The optimization of the cell is performed by searching for the Potential Energy Surface (PES) local minimum (closest to the starting structure) by varying all the parameters of the cell, i.e. unconstrained model, obtaining the optimized sides and angles of the cell, as well as the coordinates of all the atoms.The aforementioned computational methodology enabled the calculation of the cell dimensions, bond lengths, and inter-ionic distances, as well as the determination of local structural distortions of the studied systems.Furthermore, the Vegard's law was employed to estimate unit cell sides 21 according to the following equation: This empirical rule provides a description of the variation of lattice sides ( a i ) in a solid solution (i = C), which is composed of two components (A and B) with relative weights x and (1 − x) , respectively.Although this kind of empirical rules may provide useful insights it is important to note that they are not rigorously predictive and should be utilized with caution.However, in the absence of experimental data, they may serve as a useful point of comparison with results from ab initio calculations.

Ab initio calculations of elastic constants
For the elastic constant calculations we used CASTEP software packages 16 performing calculations at DFT level with PBEsol functional.The linear response of the stress vector to a given strain vector is described as: where σ and ε are the symmetric stress and strain tensors respectively, and C is the symmetric 6 × 6 matrix of elastic constants.In a cubic crystal, only three elements, C 11 , C 12 , and C 44 , are independent.To obtain the matrix C , small deformations, δ (of the order of 1% for cell sides and 2% for cell angles), are applied to the simulation (1) cell along the strain vectors, and the resulting stress vectors are calculated.The maximum force and stress are set to 0.001 eV/Å and 0.001 eV/Å 3 respectively, by employing the LBFGS optimizer.From the C 11 , C 12 , and C 44 elastic constants, several important material properties can be calculated using the Voigt-Reuss-Hill scheme 22,23 .These properties include the Bulk modulus (B), the Voigt Shear modulus (G V ), the Reuss Shear modulus (G R ), the Voigt-Reuss-Hill Shear modulus (G), Young modulus (E), and Poisson ratio (υ), which can be computed using the following Eqs.(3-8): 5][26] Furthermore, the Grüneisen parameter (γ) 27 , the Debye temperature ( θ D ) 24,25 and acoustic Debye temperature ( θ a ) 28 can be determined from well-established Eqs.(9-14): where ρ is the density of the system, n is the number of atoms in the unit cell, n p is the number of atoms in the primitive cell, β is the volume thermal expansion coefficient, V m is the molar volume, C v is the molar heat capacity, h is the Plank constant and k is the Boltzmann constant.All these quantities are expressed in International System units (SI).
Regarding the Grüneisen parameter γ (which takes into account of the average anharmonicity of the bonds), we preferred not to use the original Grüneisen formula (12), because the evaluation of Cv requires the computational evaluation of the whole phonon spectrum, which is a computational intensive and demanding task.Rather we took advantage from the Leontiev's approach, further elaborated by Belomestnykh, and Sanditov in References [29][30][31] that established a relationship between the average value of the Grüneisen parameter γ that characterizes the degree of anharmonicity of interatomic forces and the velocities of sound in an isotropic, spatially unbounded elastic medium, i.e.: (3) Using this formula, γ can be derived from experimental data of sound velocity, which are much more readily available than phonon spectra.The relationship above was verified over a broad range of crystalline solid 29 , and its validity was recently reviewed in 32 .Besides, it was validated against a broad set of experimental data 29,31,32 so we deemed it can be reliably adopted as a proxy for the original Grüneisen formula (Eq.12).

Pure crystalline solids
In a crystalline solid, thermal energy can be conveyed through the motion of both electrons and phonons, i.e. quasiparticles which describe lattice vibrations.In insulators and semiconductors, lacking free charge carriers, thermal conduction is solely mediated by phonons.As it is well known, while a hypothetical, idealized crystal would possess an infinite lattice thermal conductivity, real solids are inevitably imperfect and subject to anharmonic oscillations.These factors have a direct impact on the lifespan and behavior of phonons, ultimately affecting the thermal conductivity of the material.
In literature we can find a multitude of approaches for calculating lattice thermal conductivity.One of the most popular analytical model of lattice thermal conductivity is the Slack equation 33 : where T is the temperature, A = 0.849 • 3 3 √ 4 /(20π 3 ) is an empirical parameter, M av is the mass in kg of the unit formula divided by the number of constituent atoms and V p with δ 3 the volume per atom in m 3 .The Slack equation is a theoretical model that describes the thermal conductivity of a material at various temperatures.It assumes that three-phonon scattering is the dominant mechanism of heat transfer, meaning that the heat is transferred through the vibrations of the material atoms.This equation considers the average mass of the atoms in the material, the speed of sound in the material, and the Debye temperature, which is a measure of the thermal excitation of the material lattice vibrations.Although the Slack formula depends on the empirical parameter A (see Eq. 16), it was chosen as a reference for the calculation of the thermal conductivities of pure compounds because it is one of the most widely used approximations for the estimation of the thermal conductivity of non-metallic compounds.However, in this work we propose a functional form for A depending on γ, derived from the analysis of the thermal properties of a set of compounds based on cubic oxides reported in Table 1.In particular, the experimental values of the thermal conductivity (k exp ) of these compounds have been obtained from literature while γ values for each compounds have been calculated from the parameters B and G Table 1.Data used for the calculation of the corrective term of Slack equation in cubic oxides.The values of γ have been calculated by the bulk (B) and shear modulus (G) 34 .The thermal conductivities are expressed in Wm −1 K −1 .reported in Material Project Database 34 .As it is possible to see by inspecting Fig. 1 there is a correlation between k exp /k Slack and γ for the cubic oxides dataset analysed.This correlation can be modeled as a power function of γ giving rise to a revised Slack equation.
The revised version of Slack equation has been applied to the dataset of the compounds reported in Table 1 which reports the experimental (k exp ) and the Slack (k Slack ) thermal conductivities employed for the estimation of the corrective term A, as a function of γ calculated as reported in Eq. 16.
Thus, the modified Slack equation for cubic oxides becomes: The 0.289003•γ 4.156157 correction considers the influence of phonon scattering processes due to anharmo- nicity, characteristic of the specific material symmetry, which are not fully captured by the basic Slack model.This opens the possibility of using the Grüneisen's parameter (which depend on crystal symmetry and chemical composition) to obtain specific multiplicative function for different symmetry classes of materials.
This approach is reasonable because this parameter describes how the i-th phonon frequencies (ω i ) are affected by the cell volume (V), indeed: Equation 18shows that γ is influenced by the underlying crystal symmetry (i.e. the vibrational modes symmetry) and by the chemical composition.Eq. 17 has been applied to the compounds reported in Table 2 demonstrating a substantial improvement in predictive accuracy with respect to the classic Slack equation.The revised Slack equation exhibits a markedly higher Pearson coefficient (R 2 = 0.944), a mean relative percentage error of 19.539%, a higher correlation coefficient (μ = 0.972) and a considerably lower Mean Squared Error (MSE = 9.987) with respect to the original Slack equation (R 2 = 0.620, μ = 0.788, MSE = 122.862,and mean relative percentage error of 56.348%).These results highlight the increased precision and accuracy achieved through the refined Slack equation.

Doped solids
The presence of dopant ions in optical materials, such as crystals and ceramics, leads to the creation of lattice defects, which can in turn affect the thermal conductivity of the material.For this reason, we developed a method to predict thermal conductivity in this kind of doped materials.The method is based on a revised version of Klemens formula.
As it is well known, the original version of Klemens equation [45][46][47] predicts the thermal conductivity of lattice having some point defects considering only the mass variation on phonon transport (see Eq. 19).The Klemens equation can be expressed as: where V 0 is the volume per atom, Ŵ is the scattering parameter (see Eqs. 21-23) and x is the concentration of dopant atoms, k 0 and k 100 are thermal conductivities of solids with the compositions x = 0 and x = 1 respectively.At each composition, the parameters k m , v a , and V 0 undergo gradual adjustments through linear interpolation, using the properties of pure host A and B. Klemens firstly proposed the estimation of the scattering parameter Γ through the following Eq. 48: where the scattering parameter Ŵ is essentially the average mass variance in the system, (M i − M) 2 , relative to the square of the average mass M 2 .Another term which can consider the point defect system perturbation is the average change in atomic radius (ΔR) and the variation of the harmonic force constant (ΔK) mainly due to a change of the bond strengths with respect to the non-perturbed sites.Indeed, the introduction of dopant ions creates lattice point defects, which perturb the Hamiltonian of the system, inducing a rearrangement of the phonon probability density because of changes in the site mass (ΔM), harmonic force constants (ΔK), ionic radii (ΔR), and crystal symmetry, as depicted in Fig. 2 (left).These perturbations affect the potential energy of the material, enhancing phonon scattering due to a shortened mean free path, and ultimately decreasing the thermal conductivity.ΔM, ΔK, and ΔR are responsible for specific contributions: ΔM changes the material density, ΔK changes the lattice stiffness and ΔR changes the distance between atoms, all together influence the phonon propagation and thermal conductivity.The interplay between these factors creates a complex picture of the material thermal behavior.The relationships between force constants and atomic volumes could be described by a phenomenological fitting parameter.These were initially proposed by Abeles 49 and then modified by Wan et al. 50based on the derivation of Klemens and Callaway: where Ŵ i considers the defects due to the atoms of the chemical specie i.
For a mixed compound containing several defect atoms the total scattering parameter is given by: We insert a term in Ŵ , which considers the symmetry loss due to local distortions, see Fig. 2 (right), as the concentration of dopant cations increases and the cohesion of the material changes.
This corrective term is a power function of the parameter d: where − → R i with i = A or B is the average displacement of the atoms of the partially doped system C, in compari- son to the undoped host lattice A or B respectively; σ is the average bond length, for the compound under consideration.The variables a i are the average side of the cell with i = A , B , C referring respectively to the pure compounds A , B and to the mixed compositions samples, C, with weight x for A and is used as a normalization factor ensuring the consistency of the model and comparability across different materials.
Accordingly, Ŵ i can be written as: All the corrective terms in the case of pure compound assume a null value.Equation 25shows that the term d is integrated into Ŵ i as part of a power function with coefficient a and exponent r .The parameters a and r are extrapolated from the fitting of the difference between the value of Ŵ deriving from experimental values and Ŵ i obtained from original Klemens equation form.Taking advantage from the corrective parameter d , the proposed revised Klemens model considers the impact of point defects on lattice symmetry and the impact of the overall atoms displacement as depicted in Fig. 2.
The modified Klemens equation is thus given by:

Cost effective method for thermal conductivity prediction: overview and applications to mixed composition laser materials
This method provides a cost-effective approach for lattice thermal conductivity predictions in optical materials based on cubic oxides.We start from the pre-relaxed pure compound structure from the Materials Project database 51 that is refined through the structural relaxation in the DFT framework with PBEsol functional and Grimme D3 dispersion correction, providing a high-accurate representation of the system.The elastic constants, which characterize the material mechanical properties, are calculated using the same ab initio approach, providing a comprehensive and internally consistent set of parameters.These elastic constants are then utilized to calculate the Grüneisen parameter, Debye acoustic temperature, and acoustic velocity, which collectively characterize the thermal properties of the material (according to Eqs. 9-15).To achieve a more accurate prediction of the thermal conductivities of cubic oxide optical materials, a modified Slack equation is employed, which incorporates the multiplicative corrective factor of γ, as reported in Eq. 16.Regarding the doped compounds, starting from the pure compound optimized geometry, a doped structure is generated, and DFT-PBEsol relaxation is employed to refine its atomic coordinates.The refined structure is then utilized to calculate the corrective geometric factor, which quantifies the influence of defects on the crystal lattice symmetry and material cohesion.In order to accurately model the effects due to the point defects on thermal conductivity, a modified Klemens equation is applied (see Eq. 26).This equation is utilized to accurately predict the impact of point defects on thermal conductivity by also considering the distortion parameter d, which is a quantitative measure of the deviations from the ideal lattice symmetry due to the presence of defects.More precisely, the multiplicative factor a and the exponent r results from the best fit obtained by the minimization of the MSE between the experimental and calculated thermal conductivity values of both mixed compositions garnets and sesquioxides series of samples.Figure 3 reports in a schematic way the calculation steps of our method.
In the following paragraphs we test our method on two classes of cubic oxides materials with different compositions.In the first case study we estimate thermal conductivity in binary garnet systems while in the second case we analyze a more complex series of ternary mixture of sesquioxides both employed as laser materials.

Calculation of the lattice parameters
We calculated structural parameters for the series YAG, Lu 33 , Lu 50 , Lu 67 and LuAG.
The arrangement of the ions leads to a specific type of chemical bonding and a distinctive distribution of phonon frequencies, which have a significant impact on the thermal transport properties of these materials.As already mentioned, the addition of Lu 3+ dopants in YAG host introduces structural distortions and point defects (i.e.phonon scattering centers), leading to a decrease in the lattice thermal conductivity.The calculations of the unit cell sides are in good agreement with the experimental data 15 as reported in Table S1 (Supplementary Information).The substitution of Y 3+ with Lu 3+ ions in YAG structure disrupts the perfect symmetry of the Y 3+ ions coordination environment, leading to a modified crystal field symmetry.Also in this case, the structural variations are attributed to the difference in radius between Lu 3+ ions (85 pm) and Y 3+ ions (89 pm), leading to a contraction of the unit cell and the formation of point defects.There are two types of dodecahedral Y-O (Lu-O) bonds which are distinguished with subscripts 1 and 2 (see Table 3).Furthermore, the ions Al 3+ occupy both octahedral (Al 1 ) and tetrahedral (Al 2 ) sites.

Determination of elastic constants
We calculated the elastic constants of Lu (E = 283 GPa) and Poisson ratio (υ = 0.25) reported in ref. 52 are in quite good agreement with our calculated data (see Table 5) for Y 3 Al 5 O 12 .As regard Lu 3 Al 5 O 12 single-crystal, the reported data in Table 15 are in agreement with data reported in Ref. 53 (C 11 = 342 GPa, C 12 = 112 GPa, and C 44 = 115 GPa; B = 189 GPa, G = 115 GPa, E = 287 GPa and υ = 0.247).Table 5 reports the calculated structural and phonon parameters of the materials under study.These parameters were used to calculate the thermal conductivities of pure compounds (Y 3 Al 5 O 12 , Lu 3 Al 5 O 12 ) by modified Slack equation.

Thermal conductivities models
Through the revised Slack equation, we obtained thermal conductivities in quite alignment with experimental values as reported in Table 6.
To evaluate the effects of increasing Lu 3+ content in the samples Lu 33 , Lu 50 , Lu 67 , we used for k 0 the thermal conductivity of YAG and for k 100 that of LuAG.For the calculation of the corrective term Ŵ (Eq.24) we set c i as percent fraction of the i-th compound (i = LuAG or YAG), M i as the molar mass of the i-th compound and R i −R R  www.nature.com/scientificreports/ was replaced with R i −R mx R mx , where R i is the average ionic radius in the i-th compound, and R mx is the average ionic radius in mixed ceramic under study.In Table 7, we report the calculated distortion paramaters d YAG and d LuAG (Eq.24), the average bond length ( σ ) and the average cell side length ( a).
Table 8 reports the thermal conductivity values calculated with Klemens Eq. ( 19) and modified Klemens Eq. ( 26) compared with experimental data.In Eq. 26 we set a = −0.00505and r = −0.11348 .We indicate with Mod.k LK (d) the revised Klemens equation with YAG and LuAG experimental thermal conductivities and rev.kLK (d) the revised Klemens equation with YAG and LuAG thermal conductivities calculated by revised Slack equation.
It can be seen from the data of Table 8 and Fig. 5 that Mod.k LK equation provides a more accurate evaluation of the thermal conductivity than the original Klemens model.Regarding the method rev.kLK , its accuracy is apparently lower, but we have to take into account that in this case both the original Klemens model and Mod.k LK are seeded with the experimental data of thermal conductivity of YAG and LuAG, while in rev.kLK relies on thermal conductivity values calculated from the elastic constants.We want to point out that, while the evaluation provided by rev.kLK cannot compete in accuracy with the other two, this evaluation scheme can be used even when the thermal conductivity data of the pure compounds are not know, as long as the elastic constants and lattice parameters are available.
The statistical analysis of the data reported in Table 8 highlights the performance of each model.In particular we obtained a Percentage Relative Error of 6.606% for the original Klemens equation, 1.196% for Mod.k LK and 9.940% for rev.kLK .This data shows a significant accuracy increase of our Mod.kLK .For the rev.kLK it must be taken into account that this model is affected by the errors of both revised Klemens and revised Slack equations.
Table 5. Structural and phonon parameters of the materials under study.The density, ρ is expressed in (kg/ m 3 ), v L , v S and v a are expressed in (m/s), θ D and θ a are expressed in K and γ is dimensionless.

Calculation of the lattice parameters
To assess the accuracy of the PBEsol functional with Grimme D3 dispersion correction, we compared the modeled structural parameters of the pure Y 2 O 3 crystal to the experimental data reported in 54 .This data served as an established benchmark against which the predictive power of the functional could be evaluated.Y 2 O 3 exhibits a body-centered cubic lattice structure (see Fig. 6) with the space group Ia3 and Y 3+ sites with symmetry C 3i (Y 1 ) and C 2 (Y 2 ), and O 2-ions without symmetry (O).The PBEsol functional with Grimme D3 dispersion correction produced a cubic unit cell with a lattice side of 10.61 Å, which is in agreement with experimental data 54 .In particular, the mean bond lengths of the three different Y 1 -O bonds in the C 3i sites were determined to be 2.33 Å, 2.27 Å, and 2.24 Å, while the mean bond length of the Y 2 -O bonds in the C 2 sites was 2.28 Å.These results suggest that the PBEsol functional with Grimme D3 dispersion correction is a reliable tool for accurately modelling the structural properties of Y 2 O 3 crystals and other related materials.
In a preliminary investigation, we performed DFT calculations with PBEsol functional on 5at.%Tm:Y 2 O 3 .The Tm 3+ cations replace Y 3+ at both C 2 (Y 2 ) and C 3i (Y 1 ) sites.We observed very small perturbations in lattice cell parameters (< 0.1%), bonds lengths and inter-cationic distances (<1%).This result was not unexpected as Tm 3+ cations have an ionic radius (85.8 pm) very similar to the ionic radius of Y 3+ (89.3 pm).
To reduce computational time while maintaining structural accuracy, we utilized a structural relaxation technique wherein the unit cell angles were fixed at 90°.This approach assumes that deviations from a perfect  cubic lattice are negligible, thus allowing for an efficient optimization of the atomic positions within the crystal structure.This assumption is supported by the experimental data reported in 55 .
The lattice vector values obtained from the structural relaxation of 5at.%Tm:Y 2 O 3 are shown in Table 9.
In Table 10 the X 1 -O experimental 54 and calculated average values of bond lengths are reported.X 1 is referred to Y 3+ in Y 1 site for experimental data 54 and to Tm 3+ which replaced Y 3+ in Y 1 site for PBEsol calculated data.
Looking at Table 10 a slight contraction of two bond lengths can be observed when Tm 3+ replaced Y 3+ in the Y 1 site.However, the bond lengths variations are smaller than 1% and this is principally due to the similar ionic radius between Y 3+ and Tm 3+ .PBEsol functional have been used also to study the structural variations when an increasing concentration of Sc 3+ is added to the 5 at.%Tm:Y 2 O 3 sample.When Y 3+ is substituted with Sc 3+ a notable contraction of both the cell parameters and of the bond lengths (inter-cationic distances) is observed.This is mainly due to the differences in ion radius between Y 3+ and Sc 3+ as the ionic radius of Sc 3+ (i.e.74.5 pm) is 19% smaller than Y 3+ (i.e.89.3 pm).Due to the complexity of these ternary mixture sesquioxides, in this case, we performed a more accurate structural analysis of DFT output going to deeply inspect local symmetry losses through the Pair Distribution Function (PDF) and X-Ray Diffraction (XRD) simulations (see Fig. S1, Fig. S2 and Fig. S3 in Supplementary Information).The reported data refer to 5at.%Tm:Y 2 O 3 55 , Sc 12 , Sc 25 and Sc 50 .The calculations are carried out with 50% Sc 3+ sample (and not with 49.8% Sc 3+ as in the experiment) only for practical convenience in the cell construction.
The increasing concentration of Sc 3+ results in a contraction of the lattice vectors and, in turn, of the volume (see Table 11).The trends obtained by DFT calculations were compared with Vegard's law 21 .The lattice side of 5at.%Tm:Y 2 O 3 is a 5at.%Tm:Y2O3 =10.60 Å while the lattice side of Sc 2 O 3 is a Sc2O3 = 9.79 Å 56 .By using the Vegard's law, we have: The agreement between the lattice values calculated with Vegard's law, the ones calculated at DFT theory level and experimental values is very good, see Table 12.  9. Sides (a, b, c) angles (α, β, γ) and volume (V) of the 5 at.%Tm:Y 2 O 3 unit cell compared with experimental structural parameters 55 .All sides of the cell are given in Å, the volume in Å 3 and the angles in degrees.We also report the percentage relative deviation (Δ%) of the calculated data from the experimental values.www.nature.com/scientificreports/ The increase of the concentration of Sc 3+ increase the entropy of the system due to the formation of new Sc-O bonds, which add new configurations to the system and increase the number of possible arrangements of the atoms (see Fig. S1); in other words, the number of complexions rises.Accordingly, the partial loss of symmetry of the C 3i and C 2 sites in Sc-doped 5at.%Tm:Y 2 O 3 can be attributed to the atomic-scale distortions caused by the Sc 3+ ions.The substitution of Y 3+ with Sc 3+ introduces local perturbations in the crystal lattice structure, which lead to a decrease in the crystal field symmetry around these sites, resulting in a lower degree of the symmetry coordination with the surrounding oxygen atoms.Three types of bond lengths in the doped crystal are possible, see Table 13.
It is clearly observed a local change of the structure with increasing Sc 3+ concentration through the X-ray atomic PDF; the latter is an X-ray scattering technique which can be used to study the local structure of materials on an atomic scale.It is based on the concept of radial distribution function, which measures the probability of finding an atom at a certain distance from another atom.By analyzing the PDF output, it is possible to extract information on the average interatomic distances, bond lengths, bond angles, coordination numbers, as well as local symmetry of the material.We simulated PDF with EXPO 57 , using as input the optimized structures of 5at.%Tm:Y 2 O 3 , Sc 12 , Sc 25 and Sc 50 because PDF is a powerful tool for probing the structural properties of crystalline and amorphous materials.In our calculation PDF is used to describe the distribution of ion pairs (Y1-O, Y2-O, X-O, Sc-O) within the volume occupied by the system.
We show how the system partially loses the recognition of the C 3i and C 2 sites by increasing the concentration of Sc 3+ leading to an intensification of the formation probability of the Sc-O bond with length in the interval 2.126-2.153Å (see Fig. S1).
Accordingly, the calculated XRD patterns of the samples Y 2 O 3 and Sc 50 show the shift of the most intense peaks and the appearance of new ones (see Fig. S2).
The comparison between the experimental 58 and calculated spectra of the Sc 50 sample (see Fig. S3) confirms the good reliability of the structural data calculated with PBEsol which returns an accurate picture of the local disorder introduced by the Sc 3+ cations as well as the effect of contraction of the cell parameters.

Determination of elastic constants
By using the same computational framework 16 59 are in good agreement with our calculated data (see Table 14) for 5at.%Tm:Y 2 O 3 ; moreover, they are also comparable with data reported in Ref. 60 (B = 153.8,G = 62.4,E = 165.0 and υ = 0.320).Table 15 summarizes the calculated structural and phonon parameters of the materials under study.It is worth to note these parameters will be used to calculate the thermal conductivities of pure compounds (Y 2 O 3 , Sc 2 O 3 ) by modified Slack equation.

Thermal conductivities models
The results of the thermal conductivities of pure Yttria and Scandia calculated by using the modified Slack equation (see Eq. 16) shows a good agreement with the experimental data 36,61   three compounds are involved, contrasting with the simpler binary systems in Lu-doped Yttrium Aluminum Garnet.We utilize a value of a = 0.027855 and r = -0.11352.For the calculation of the corrective term Ŵ (see Eq. 25) c i was set as percent fraction of dopant, M i as the molar mass of the dopant and R i −R R was replaced with R i −R mx R mx , where R i is the average radius of the ions in the host, R mx is the average radius of the ions in the mixed ceramic under study.In Table 17, the distortion parameters d Y 2O3 , d Sc2O3 and d Tm2O3 (calculated according to Eq. 24), the average bond length ( σ ) and the average cell side length ( a ) are listed.The modified Klemens and original Klemens models show trends in good agreement with the experimental data 61 .Anyway both the proposed Mod.Klemens and rev.Klemens method show the best agreement with respect to the experimental trend.In our vision, the introduction of the parameter d in the model allows a more accurate description of these mixed composition compounds because it takes into account the lattice deformation due to the point defects.This is further confirmed by the statistical analysis of the data reported in Table 18.The estimation of the Percentage Relative for the original Klemens equation is 16.793% which is significantly higher with respect to Mod.k LK (0.687%) and rev.kLK (2.562%).We can confirm that our implemented method allows a more accurate prediction of thermal conductivity in this class of materials in comparison with the respective classical models.

Conclusions
This study introduces a computational framework designed to estimate the thermal conductivity of optical materials, specifically targeting cubic oxides.Our methodology encompasses a series of comprehensive steps utilizing DFT calculations to extract essential structural parameters and elastic constants.These parameters constitute key inputs for refined the Slack and Klemens equations, enhancing the precision of lattice thermal conductivity predictions in these materials.Moreover, the above mentioned parameters can be also obtained from free on line materials database further shortening the time and the costs of the information achievement.Our approach involves the modification of the Slack equation by introducing a corrective factor of the Grüneisen parameter, i.e. ~ 0.3•γ 4 .Also, Klemens equation has been revised by the introduction of a distortion parameter, d, accounting for the impact of point defects on lattice symmetry, a critical factor significantly influencing the thermal conductivity.
The analysis of the outputs of our revised method in comparison with those of the classical Slack and Klemens equations gives rise to a substantial improvement in the accuracy of thermal conductivity prediction for the class of compound analyzed in this work.
Our model is intended to provide a good estimation of thermal conductivity of poorly characterized materials starting from experimental data that are more readily available than thermal conductivity itself, such as lattice parameters and mechanical properties.Thermal properties of the analysed class of materials are particularly important for optical applications and specifically for candidate laser hosts.Further investigations addressed on testing the proposed model also on materials belonging to different space group symmetry are in progress.

Figure 1 .
Figure 1.Fitting (black line) of the ratio between the experimental and the Slack calculated thermal conductivities (red points) as a function of the Grüneisen parameters.The exact expression for the fitting curve is 0.289003 •γ 4.156157 .

Figure 2 .
Figure 2. Left: Schematic representation of lattice perturbation due to the presence of point defects i.e. mass difference (ΔM = |M − m| ), constant harmonic force difference (ΔK), radius difference (ΔR); right: Schematic representation of the symmetry loss of the site.

Figure 3 .
Figure 3. Schematic representation of the method proposed for the prediction of thermal conductivity in both pure and doped cubic oxides.

Figure 4 .
Figure 4. Lattice structure of YAG .Yttrium sites with symmetry D 2 are colored in cyan while Aluminum sites with symmetries C 3i and S 4 in light blue, the oxygen atoms in red.

Figure 5 .
Figure 5. Experimental (black) calculated (blue, red and green) values of thermal conductivity for the entire set of samples.

Figure 6 .
Figure 6.Lattice structure of Y 2 O 3. Yttrium sites with symmetry C 3i and C 2 are colored in dark green.The symmetry of the sites is reported in round brackets.

Figure 7 .
Figure 7. Experimental (black) and calculated (red, blue and green) values of thermal conductivity for the entire set of samples based on 5% Tm 3+ doped mixed Y 2 O 3 /Sc 2 O 3 ceramics.

Table 3 .
Bond distances (after geometry optimization) for the two types of dodecahedral Y-O (Lu-O) bonds which can be identified with subscripts 1 and 2. Al 3+ ions can occupy both octahedral (Al 1 ) and tetrahedral (Al 2 ) sites.

Table 4 .
Mechanical properties of Y 3 Al 5 O 12 and Lu 3 Al 5 O 12 pure crystals expressed in GPa.

Table 7 .
d YAG and d LuAG are the distortions parameters calculated with respect to YAG and LuAG.a is the unit cell side average value, σ is the mean bond length.All the reported quantities are reported in Å.

Table 10 .
Bond lengths X 1 -O, X 2 -O, angle {X 1 OY 2 } and percentage of the relative deviation (Δ%) of the calculated data from the experimental values.More precisely, X 1 -O refers to 3 different types of bond lengths within the C 2 site.All bond lengths are given in Å and the angles in degrees.
. In Table16the experimental and calculated thermal conductivities of the pure samples are reported.Hereon, k 0 and k 100 will label the thermal conductivities of Y 2 O 3 and Sc 2 O 3 as reported in Table16.The ther- mal conductivity characteristics of 5at.%Tm 3+ -doped (Y,Sc) 2 O 3 transparent ceramics present a more intricate scenario compared to those of Lu-doped Yttrium Aluminum Garnet.As a matter of fact, in Sc 12 , Sc 25 and Sc 50

Table 12 .
Side a values calculated with Vegard's law and the mean values of the lattice side calculated with PBEsol functional both compared with experimental data.All values are expressed in Å.

Table 13 .
Average bond lengths Sc-O, X-O and Y-O (X can be Y 3+ or Sc 3+ ).These data are extracted from PBEsol calculations and are expressed in Å.

Table 18
reports the thermal conductivity values calculated by Klemens equation (see Eq. 19) and modified Klemens equation (see Eq. 26) compared with experimental data.With d in round parenthesis, we identify the revised Klemens results.To calculate the Eq. 26 a = 0.027858 and r = −0.113485were set.Mod.k LK (d) refers to revised Klemens equation with Y 2 O 3 and Sc 2 O 3 experimental thermal conductivities and rev.kLK (d) refers to revised Klemens equation with Y 2 O 3 and Sc 2 O 3 thermal conductivities calculated by revised Slack equation.Figure 7 report the comparison between experimental and calculated values of thermal conductivity obtained for 5 at.%Tm 3+ doped mixed Y 2 O 3 /Sc 2 O 3 ceramics.

Table 14 .
Mechanical properties of 5 at.%Tm 3+ -doped Y 2 O 3 , Sc 2 O 3 , Sc 12 , Sc 25 , Sc 50 (see section Ab initio calculations of elastic constants), expressed in GPa, with the exception of υ, which is a dimensionless number.

Table 15 .
Structural and phonon parameters of the materials under study.ρ is expressed in (kg/m 3+ ), v L , v S and v a are expressed in (m/s), θ D and θ a are expressed in K and γ is dimensionless.

Table 17 .
Structural parameters:d Y 2O3 , d Sc2O3 and d Tm2O3 are the distortion parameters with respect to undoped hosts (Y 2 O 3 and Sc 2 O 3 ); a is the unit cell side average value and σ is the mean bond lengths.All the values are reported in Å.

Table 18 .
Experimental and calculated values of thermal conductivity (expressed in Wm −1 K −1 ).