A first-principles investigation of the linear thermal expansion coefficients of BeF2: giant thermal expansion

We present the results of a theoretical investigation of the linear thermal expansion coefficients (TECs) of BeF2, within a direct Grüneisen formalism where symmetry-preserving deformations are employed. The required physical quantities such as the optimized crystal structures, elastic constants, mode Grüneisen parameters, and phonon density of states are calculated from first-principles. BeF2 shows an extensive polymorphism at low pressures, and the lowest energy phases [α-cristobalite with space group (SG) P41212 and its similar phase with SG P43212] are considered in addition to the experimentally observed α-quartz phase. For benchmarking purposes, similar calculations are performed for the rutile phase of ZnF2, where the volumetric TEC (αv), derived from the calculated linear TECs along the a (αa) and c (αc) directions, is in very good agreement with experimental data and previous theoretical results. For the considered phases of BeF2, we do not find any negative thermal expansion (NTE). However we observe diverse thermal properties for the distinct phases. The linear TECs are very large, especially αc of the α-cristobalite phase and its similar phase, leading to giant αv (∼175 × 10−6 K−1 at 300 K). The giant αv arises from large Grüneisen parameters of low-frequency phonon modes, and the C13 elastic constant that is negatively signed and large in magnitude for the α-cristobalite phase. The elastic constants, high-frequency dielectric constants, Born effective charge tensors, and thermal properties of the above phases of BeF2 are reported for the first time and hence serve as predictions.


Introduction
Beryllium uoride (BeF 2 ) is known to exist in glass and crystalline phases and has a variety of technological applications. BeF 2 glass has a large bandgap of about 13.8 eV, the lowest refractive index and highest Abbe number of any inorganic material, and exceptional resistance to damage. These properties have enabled the manufacturing of special glasses (from BeF 2 and its mixtures with uorides and other diuorides) that have excellent transmittance in the UV region 1,2 and for use in high-power laser systems. 3 The LiF-BeF 2 mixture is a primary coolant and fuel solvent in molten salt nuclear reactors. 4 In protein crystallography, BeF 2 is used to restrict protein motion to facilitate the crystallography process. 5 Very recently, crystalline BeF 2 is predicted to be a better neutron lter than MgF 2 , which has been considered an effective neutron lter candidate. 6 The main aim of this work is to investigate, for the rst time, the linear thermal expansion coefficients (TECs) of a few low-energy crystalline phases of BeF 2 . We also consider a benchmark system ZnF 2 that has exceptional electric and optical properties, and interesting technological applications ranging from catalysis to spectroscopy and laser applications. 7 Single crystal BeF 2 has been grown and found to have a crystal structure remarkably similar to that of the a-quartz (SiO 2 ) structure, 8 which has a trigonal symmetry with space group (SG) P3 1 21 (#152). A recent rst-principles study has revealed that BeF 2 shows extensive polymorphism at low pressures. 9 Interestingly, three crystal phases [namely, (i) the acristobalite phase that has a tetragonal symmetry with SG P4 1 2 1 2 (#92), (ii) a similar phase to the a-cristobalite phase (hereaer referred to as the a 0 -cristobalite phase) with SG P4 3 2 1 2 (#96), and (iii) the C2/c-2 Â 4 phase with SG C12/c1 (#15)] are predicted to be energetically more stable than a-quartz. However, these phases have a very small stability pressure range (less than 0.4 GPa), and the a-quartz phase transforms to the coesite-I phase SG C2/c at 3.1 GPa. The high-pressure phases of BeF 2 have been the subject of other rst-principles calculations. 10 Very recently, rst-principles calculations have also been employed to construct the P-T phase diagram of BeF 2 . 11 The HSE06 optical bandgap of the a-quartz structure is found to be about 10.6 eV, and increases by increasing the applied pressure. 9 The lattice vibrations, inelastic scattering crosssections, and neutron transmission of BeF 2 have been thoroughly investigated 6 using rst-principles calculations and compared to those of MgF 2 .
The benchmark system ZnF 2 crystallizes in the tetragonal rutile structure with SG P4 2 /mnm (#136). Very recently, Raman scattering measurements with the use of the diamond anvil cell have been employed to investigate the structural phase transformations of ZnF 2 under high pressures. 12 This experimental work is supplemented by rst-principles calculations. In addition to the structural stability and pressure variation of the Raman active phonon modes, the electronic bandgap of the considered phases as a function of pressure has been investigated at the HSE06 level. Neutron diffraction has been employed to study the temperature dependence of the lattice parameters and unit cell volume of ZnF 2 , 13 and NTE has been observed in a small temperature range (below 75 K). This NTE behavior has been supported by rst-principles calculations. 13,14 However, only the volumetric TEC has been theoretically investigated.
In the present work, the linear TECs of BeF 2 and ZnF 2 are investigated by employing the recently introduced direct approach 15 in which the symmetry of the deformed structures could be preserved. Since this approach has not been applied to systems with tetragonal symmetry, ZnF 2 is thus chosen as a suitable benchmark system because of its tetragonal crystal structure, in addition to the existence of experimental and previous theoretical results of its volumetric thermal expansion. The elastic constants and phonon frequencies required to compute linear TECs are calculated from rst-principles. For BeF 2 , the a-quartz, a-cristobalite and a 0 -cristobalite phases will be considered. Moreover, the relative stability of the above three phases of BeF 2 are also investigated using different levels of approximation of the exchange-correlation potential.

Methodology
The linear TECs of the considered phases of ZnF 2 and BeF 2 are calculated within the Grüneisen formalism following the procedure described in ref. 16-22. To compute the mode Grüneisen parameters we considered two types of symmetrypreserving deformations obtained by changing the in-plane (a) or out-of-plane (c) lattice parameters by AE0.5%. These deformations allow for the full utilization of the tetragonal or trigonal point-group symmetry 23 of the considered systems, which minimizes the required number of independent atomic displacements (i.e., number of supercells) to calculate the phonon frequencies within the direct method. [24][25][26][27][28] The amplitude of atomic displacements, from the corresponding equilibrium positions, is 0.015Å. The supercell sizes are of 2 Â 2 Â 3 for ZnF 2 , and 2 Â 2 Â 2 for the a-quartz, a-cristobalite, and a 0cristobalite phases of BeF 2 . The adequacy of these supercells have been checked by considering larger ones for each of these systems, and we found that such actions do not alter appreciably our main results and conclusions. The determination of the linear TECs also requires the elastic constants that may be obtained through ttings of energy versus strain curves. 29,30 Specically, these TECs at temperature T in the a (a a ) and c (a c ) directions of the above systems are given by " a a ðTÞ a c ðTÞ where C ij are the elastic constants, D ¼ (C 11 + C 12 )C 33 À 2C 13 2 and U is the volume of the primitive cell. The phonon density of states (PDOS) weighted by the Grüneisen parameters are with the Grüneisen parameter g i,lq ¼ Àn lq À1 vn lq /v3 i for the deformation of type i (with i ¼ a for in-plane and i ¼ c for out-ofplane deformations). The derivative vn lq /v3 i measures the change of the frequency n lq with respect to the strain parameter 3 i . 15 The summation is over all frequencies n lq for the phonon band index l and q vector in the Brillouin zone (BZ). The I i (T) are calculated from where c(n,T) ¼ k B (r/sinh r) 2 is the contribution of the phonon modes with frequency n to the specic heat. Here, r ¼ hn/2k B T, and h and k B are the Planck and Boltzmann constants, respectively. The maximum (minimum) frequency is denoted by n max (n min ). The DFT calculations of the optimized structural parameters, phonon frequencies, and elastic constants are performed by employing the projector augmented wave (PAW) method, as implemented in the Vienna Ab Initio Simulation Package (VASP). A relatively high cutoff energy of 600 eV is used for the plane-wave basis. Geometry optimization is stopped when the maximum force on each atom is less than 10 À3 eVÅ À1 . We nd that phonon frequency shis are more consistent when we use the local density approximation (LDA) for BeF 2 , and the PBE_sol functional of the generalized gradient approximation (GGA) for ZnF 2 . Therefore, for the linear TECs only the results of these calculations are reported.

Structural properties
We start with the benchmark system, ZnF 2 , which at ambient conditions crystallizes in the rutile structure, shown in Fig. 1(a). This crystal structure has a tetragonal symmetry and six atoms per primitive unit cell. The two Zn atoms occupy the Wyckoff 2a(0, 0, 0) sites and the four F atoms are located at the 4f(x, x, 0) sites. Therefore, this structure is characterized by three crystallographic parameters: two lattice parameters (a and c) and an internal parameter for the coordinates of the four F atoms (x). The structural parameters obtained using the PBE_sol functional are (a, c, x) ¼ (4.7194Å, 3.1376Å, 0.3037), while the corresponding LDA results are (4.6373Å, 3.0990Å, 0.3033). These calculated values are in very good agreement with the experimental data (4.7038Å, 3.1336Å, 0.3035) 32 and (4.7034Å, 3.1335Å, 0.303). 13 As expected, the LDA lattice parameters are underestimated while those of the PBE_sol are slightly overestimated.
For BeF 2 , the three considered crystal structures are a-quartz [ Fig. 1 and a 0 -cristobalite. The aquartz phase has trigonal symmetry and nine atoms in the conventional hexagonal unit cell. The three Be atoms occupy the Wyckoff 3a(x 1 , x 1 , 0) sites and the six F atoms occupy the 6c(x 2 , y 2 , z 2 ) sites. On the other hand, the a-cristobalite phase has a tetragonal symmetry and twelve atoms per primitive unit cell. The four Be atoms occupy the Wyckoff 4a(x 1 , x 1 , 0) sites and the eight F atoms occupy the 8c(x 2 , y 2 , z 2 ) sites. Therefore, each of these two structures has six crystallographic parameters: two lattice parameters (a and c), and four internal parameters (denoted as x 1 , x 2 , y 2 and z 2 ). The atomic coordinates of the a 0cristobalite phase can be obtained from those of the a-cristobalite by mirror-image transformation (x, y, z) / (Ày, Àx, z), and the lattice parameters of the two structures are identical. Therefore, only the structure parameters of the rst two crystal structures are reported. Our LDA, PBE_sol, and PBE results, shown in Table 1, are in good agreement with available experimental data and other theoretical calculations. The PBE_sol results lie between the corresponding LDA and PBE results and show the best agreement with the experimental data for the a-quartz.

Elastic properties and stability
We show in Table 1 the LDA, PBE_sol and PBE relative energies (DE) of the a-cristobalite phase with respect to those of a-quartz. According to the PBE and PBE_sol calculations the latter phase is slightly more stable, in accordance with the Nelson et al. 9 GGA calculations. However, the LDA calculations lead to an opposite conclusion. Similar conclusions have recently been reached by Masoumi, 11 using both LDA and GGA calculations. These results show that these two phases have extremely close cohesive energies.
The a-quartz phase of BeF 2 with a trigonal crystal symmetry has six independent elastic constants. 34 On the other hand, the rutile phase of ZnF 2 and a-cristobalite phase of BeF 2 (both have Fig. 1 The crystal structures of (a) rutile ZnF 2 , (b) a-quartz BeF 2 , and (c) a-cristobalite BeF 2 . The a 0 -cristobalite BeF 2 structure which is similar to the a-cristobalite BeF 2 structure is not shown. The small gray balls represent the F atoms. The c axis of the three crystal structures is along the vertical direction. Table 1 Calculated lattice constants and internal parameters of the a-quartz and a-cristobalite phases of BeF 2 . Also shown are the relative energy (DE) of the a-cristobalite phase with respect to that of a-quartz, and the available experimental data (measured at 100 K) and other theoretical results

Phase Approach
Lattice constants Internal parameters  31 4.695 6.318 LDA 11 4.684 6.373 PBE 11 4.960 6.910 a tetragonal crystal symmetry) also have six independent elastic constants. 34 The elastic constants of these phases, obtained by using the LDA and PBE_sol functionals are shown in Table 2.
There are two features to note from this The elastic constants could be used to investigate the mechanical stability of the crystal structure. For a-quartz structure, the Born stability criteria 34 are and (C 11 À C 12 )C 44 À 2C 14 2 > 0.
Note that the same expression of D appears in eqn (1). As for the rutile and a-cristobalite phases (of tetragonal (I) class), 35 the Born stability criteria are: C 11 > jC 12 j, D > 0, C 44 > 0 and C 66 > 0. The elastic constants reported in Table 2 show that these criteria are satised, and therefore the considered phases of ZnF 2 and BeF 2 are mechanically stable. The mechanical stability of these crystal structures can also be inferred from phonon dispersion relations, discussed below.

Phonon dispersion relations
Since the considered crystals are polar in character we perform non-analytical correction (NAC) to their dynamical matrices. To do that, we have calculated the high-frequency dielectric constant and Born effective charge tensors, and the results are listed in Table 3. The features to note from this table are the following. (i) The calculated results have a very weak dependence on the used exchange-correlation functional. (ii) The effective charges in the ZnF 2 are larger than in BeF 2 , which shows that the ionicity of Zn-F bond is larger than that in Be-F. This is consistent with the corrected Allred-Rochow electronegativity values 39 (larger for Be). (iii) Our calculated xx and yy components of the dielectric constant are in very good agreement with available experimental data, 36 while that of zz is larger than the measured one. However, the comparable values of diagonal components of the dielectric constant in our calculations are consistent with the experimental and calculated values for other metal uorides crystallizing in the rutile structure (such as MgF 2 and FeF 2 ). 40 Fig. 2 shows the calculated phonon dispersion relations and PDOS of the rutile ZnF 2 , with and without NAC. Also shown are  Fig. 2. For the designation of these phonon modes see ref. 40. Fig. 2 shows that these experimental data agree reasonably well with our rst-principles results. Fig. 3 shows the phonon dispersion relations of the a-quartz and a-cristobalite phases of BeF 2 , taking into account the NAC. Also shown are the PDOS, and Be and F projected PDOS of the acristobalite phase. The results of the a 0 -cristobalite phase are very similar to those of a-cristobalite and hence are not shown. The features to note from this gure are the following. (i) The very wide frequency range of the phonon modes in these systems, compared to that of ZnF 2 . This can be understood as a consequence of the rather large mass difference between Be and Zn atoms. (ii) The frequency range of both phases of BeF 2 can be separated, according to the character of the phonon modes, into three sub-regions. (a) The lower frequency region between 0 and about 700 cm À1 , where the phonon modes are mainly due to the vibrations of F atoms. The contribution of the Be atoms becomes appreciable above 300 cm À1 . It is worth noting that in the case of ZnF 2 the dominance of the vibrations of the F atoms occurs in the upper part of the frequency range because the Zn atom is heavier than the F atom. (b) A narrow intermediate region at about 770 cm À1 , where the rather localized phonon modes originate from vibrations involving both Be and F atoms. (c) The upper-frequency region, where the phonon modes originate mainly from vibrations of Be atoms. (iii) The opening of two frequency gaps, between (a) and (b), and between (b) and (c) sub-regions. These frequency gaps can be understood as a consequence of the localization of phonon modes in the (b) sub-region.

Thermal expansion
The calculated linear and volumetric TECs of the considered phases of ZnF 2 and BeF 2 , according to the procedure described in Sec. 2, are depicted in Fig. 4 and 5, respectively. It is worth mentioning that, as expected, the NACs to the dynamical matrices have negligible effects on the calculated linear TECs.
We will rst look at the TECs of ZnF 2 . The important features to note from Fig. 4 are the following. (i) The NTE at low temperatures is mostly due to a a . The negative values of a c are smaller (in magnitude) than those of a a and lie in a considerably shorter T-range. This is clear from the magnitude and the location of the minimum values: a a $ À1.05 Â 10 À6 K À1 at 55 K, and a c $ À0.5 Â 10 À6 K À1 at 40 K. These results are consistent with observed T-variations of the a and c lattice parameters at low temperatures (see Fig. 3 of ref. 13). (ii) The calculated a v from the linear TECs (i.e., a v ¼ 2a a + a c ) are in good agreement with the previous direct theoretical calculations, 13,14 and the results of all these theoretical calculations are in a qualitative agreement with experimental data. 13 This nding reects the accuracy and reliability of our calculated linear TECs. (iii) a c is systematically and appreciably larger than a a . For example, at 300 K the calculated value of a c (of 8.9 Â 10 À6 K À1 ) is about 60% larger than that of a a (of 5.6 Â 10 À6 K À1 ).
As for the thermal expansion of the considered phases of BeF 2 , the features to note from Fig. 5 are the following. (i) Unlike ZnF 2 , the calculated values of both a a and a c are always positive for all of the considered phases of BeF 2 . (ii) Both a c and a a of the a-cristobalite structure are very close to those of a 0 -cristobalite, and hence only those of the former phase will be discussed below. (iii) In the considered T-range, a a (T) of the a-quartz structure is slightly larger than a c (T), whereas a a (T) of the acristobalite phase is much smaller than a c (T). (iv) The large a c and a a lead to very large a v for both phases of BeF 2 . For example, at 300 K, the values of a v are of 77.6 and 169.9 Â 10 À6 K À1 respectively for the above two phases of BeF 2 , compared to that of 20.0 Â 10 À6 K À1 for ZnF 2 . Our largest calculated linear TEC is of $95 Â 10 À6 K À1 at 300 K for a c of the a-cristobalite phase. This is indeed large compared to the experimental linear TECs at 300 K of four uorites, i.e., CaF 2 , SrF 2 , BaF 2 , and PbF 2 (ref. 41) that range between 18.1 and 29 Â 10 À6 K À1 , but still is somewhat smaller than the measured linear TEC value of 163.9 Â 10 À6 K À1 of an Ti-Nb alloy. 42 The key physically insightful quantity for the interpretation of the above results is the PDOS weighted by the Grüneisen parameters, G i (n), dened in eqn (2). Fig. 6 shows G i (n) of ZnF 2 , and the a-quartz and a-cristobalite phases of BeF 2 . The important features to note from this gure are the following. (i) For ZnF 2 , the low-frequency modes (n < 150 cm À1 ) have negative Grüneisen parameters, which lead to negative G i (n) in this n range. Since low-frequency modes are easily thermally excited, this nding explains the observed NTE in ZnF 2 . Moreover, by inspecting the differences between G a (n) and G c (n) one can easily understand why a a is always lower than a c . (ii) The G i (n) of the considered phases of BeF 2 are always positive, which reects the dominance of positive mode Grüneisen parameters in these phases. This explains the lack of NTE in the considered phases of BeF 2 . (iii) The G a (n) and G c (n) of a-quartz BeF 2 have comparable magnitudes, with G c (n) being smaller than G a (n) for n < 100 cm À1 , which explain the comparable magnitudes and ordering of its a c and a a . (iv) The peak in G c (n) of the a-cristobalite phase around n $ 34 cm À1 is much higher than that of the G a (n), which results in a large a c compared to a a . This nding means that, in this n-range, the positive mode Grüneisen parameters associated with the out-of-plane deformation are signicantly larger than those associated with the in-plane deformation. The large Grüneisen parameters can be viewed as a manifestation of strong anharmonic effects in the a-cristobalite and a 0 -cristobalite structures of BeF 2 . However, it should be noted that large Grüneisen parameters are not the only factor that is responsible for the giant a v of the a-  Calculated linear and volumetric TECs of ZnF 2 using the LDA, compared to the available experimental data 13 and previous theoretical results. 13,14 cristobalite: the elastic property via the negative (and with a large magnitude) C 13 elastic constant (see Table 2 and eqn (1)) plays also a major role. The above ndings explain the much larger volumetric TEC of the a-cristobalite phase of BeF 2 , compared to that of the a-quartz phase, which, in turn, is larger than that of ZnF 2 .

Summary
First-principles calculations are performed to investigate the structural, elastic, and vibrational properties of the rutile structure of ZnF 2 and three crystal structures of BeF 2 (a-quartz, a-cristobalite and its similar phase with space group P4 3 2 1 2). The so-obtained phonon density of states, mode Grüneisen parameters, and elastic constants are used to study the linear thermal expansion coefficients (TECs) of the compounds mentioned above, within a Grüneisen formalism. We have used deformations that preserve the symmetry of the crystal to obtain the Grüneisen parameters. The considered crystal structures of both systems are found to be mechanically stable. The calculated physical quantities for both systems are in very good agreement with the available experimental data and previous theoretical results. For ZnF 2 , the calculated linear TECs, a a and a c , along the a and c directions are consistent with the experimental T-variations of the corresponding lattice parameters, respectively. The volumetric TEC a v computed from these linear TECs is in qualitative agreement with experiment at low temperatures, including negative thermal expansion (NTE) behavior. The considered phases of BeF 2 are not NTE materials, and their linear TECs are much higher than those of ZnF 2 , especially for the a-cristobalite phase. The elastic constants, high-frequency dielectric constants, Born effective charge tensors, and TECs of the considered phases of BeF 2 are reported in this work for the rst time and could serve as predictions.

Conflicts of interest
There are no conicts to declare.