Mechanism of Negative Thermal Expansion in Monoclinic Cu2P2O7 from First Principles

Negative thermal expansion (NTE) materials generally have high-symmetry space groups, large average atomic volumes, and corner-sharing octahedral and tetrahedral coordination structures. By contrast, monoclinic α-Cu2P2O7, which has a small average atomic volume and edge-sharing structure, has been reported to exhibit NTE, the detailed mechanism of which is unclear. In this study, we investigate the A2B2O7 polymorphs and analyze the NTE behavior of α-Cu2P2O7 using first-principles lattice-dynamics calculations. From the polymorphism investigation in 20 A2B2O7 compounds using 6 representative crystal structures, small A and B cationic radii are found to stabilize the α-Cu2P2O7-type structure. We then analyze the NTE behavior of α-Cu2P2O7 using quasi-harmonic approximation. Our calculated thermal expansion coefficients and anisotropic atomic displacement parameters were in good agreement with those of the experimental reports at low temperatures. From the mode-Grüneisen parameter distribution plotted over the entire first-Brillouin zone, we found that the phonon contributing most significantly to NTE emerges not into the special points but between them. In this phonon mode, the O connecting two PO4 tetrahedra rotates, and the Cu and O vibrate perpendicular to the bottom of the CuO5 pyramidal unit, which folds the ac lattice plane. This vibration behavior can explain the experimentally reported anisotropic NTE behavior of α-Cu2P2O7. Our results demonstrate that the most negative mode-Grüneisen parameter contributing to NTE behavior is not always located on high-symmetry special points, indicating the importance of lattice vibration analyses for the entire first-Brillouin zone.

many computational studies of lattice dynamics and molecular dynamics have reported the calculation results within the PBEsol functional [S5-S8] or the PBE functional [S9-S13].

PAW data sets used in the first-principles calculations
In the present calculations, the PAW data sets enumerated in Table S2 were used with a cutoff energy of 550 eV for the plane-wave basis.

Antiferromagnetic structure for Cu2P2O7
As shown in Figure S1, we show the antiferromagnetic configuration used for Cu2P2O7.Note that the antiferromagnetic configuration was the most stable among the three possible antiferromagnetic configurations in the primitive Cu2P2O7, which is consistent with the previous calculation results [S4, S14].

The correlation effect owing to +U correction for the d electrons in Cu 2+ ions
The formal charge of Cu in Cu2P2O7 is 2+, which gives a d 9 electronic configuration for Cu ions.The electrons correlation effect of d 9 in Cu 2+ has been extensively studied since the discovery of high-temperature cuprate superconductors.In the field of strongly correlated electrons systems, it is widely known that a strong electronic correlation generally gives a localization effect to the electrons, leading a system to a Mott insulator and/or a magnetic ordering [S15].However, our calculation results of Cu2P2O7 have clearly shown that the magnetic moments and antiferromagnetic structure remain even without +U correction, indicating that the +U correction in 3d electrons in Cu 2+ ions is not necessarily essential.This is probably due to the superexchange interaction of Cu2O6 magnetic dimers [S16].In short, considering the magnetic structure is much more important than the +U correction.To show the validity of our insistence, we also calculated the anisotropic atomic displacement parameters of oxygen atoms at 200 K with and without the +U correction, as shown in Figure S2.
Here, the Dudarev formulation was adopted for the +U correction [S17].From these results, little change can be observed by applying the +U correction.In addition, we also show the calculated phonon bands with and without the +U correction for the 3d electrons in Cu, as illustrated in Figure S3.The little change in the low-frequency region owing to the +U correction can be observed, while the increase of phonon frequency in the high-frequency region can also be observed.The increase of phonon frequency should stem from the shrinkage of lattice constants due to the +U correction (compare the calculated lattice constants of PBEsol and PBEsol+U5 in Table S1).Furthermore, the atomic displacement parameters are largely affected mainly by the low-frequency phonons, which can explain the reason why the +U correction gives little change in them (Figure S2).S6

The determination of the most stable phase in the comparison of polymorphs
For the comparison of polymorphs for A2B2O7 (A = Cu, Zn, Mg, Sr, Sc, Y; B = Ge, Sn, Ti, Zr, P, V, Ta), the total energies of the C2/c and C2/m phases for Mg2P2O7, Sc2Ti2O7, and Mg2Ta2O7 are comparable.Therefore, we determined the most stable phase of these three compounds as follows: as the most stable phase, we adopted the C2/m phase for Mg2P2O7, which is experimentally reported [S18], and the C2/c (C2/m) phase for Mg2Ta2O7 (Sc2Ti2O7) because the total energy was slightly lower than that of the C2/m (C2/c) phase.
We also present the phonon bands of the Cu2P2O7-type C2/c phases for Cu2P2O7, Cu2V2O7, Zn2V2O7, Sc2Zr2O7, Cu2Ta2O7, Zn2Ta2O7, and Mg2Ta2O7 as shown in Figure S4.We can see that the three C2/c phases of Cu2P2O7, Cu2V2O7, and Zn2V2O7 are found to be dynamically stable [Figure S4a-c  S7

The cationic radii for the comparison of polymorphs
In the map of structural stability (Figure 3), we evaluated the cationic radii of A and B by using Shannon's ionic radii [S19].We present the coordination number dependence of ionic radius as shown in Figure S3.Here, we linearly approximated the correlation between the coordination number and ionic radius, and then we evaluated the cationic radii A and B for the relevant polymorphs by using the linear regression model and the effective coordination number extracted by using VESTA [S20].S8

Comparison of calculated phonon frequencies with the experimental reports
To show the validity of the calculated phonon frequencies in the present study, we compared them with the experimental reports.The calculated phonon frequencies at Γ point and experimentally reported Raman frequencies are enumerated in Table S3.There are 3 acoustic and 63 optic modes; the phonon modes that are active to Raman and IR spectroscopy are 16 Ag + 17 Bg and 15 Au + 15 Bu, respectively.One can see that our calculation results are compatible with the experimental reports.The little difference in phonon frequencies between the calculation results and experimental reports is owing to the difference in lattice constants (see Table S1).

The distribution of mode-Grüneisen parameter on a different band path
As shown in Figure S6, we show the phonon band and the distribution of mode-Grüneisen parameters for the C2/c phases of Cu2P2O7 along with the band path suggested by Hinuma et al [S22].In Figure 5, the band path passes the reciprocal point (1/4, 1/4, 1/4), whereas the band path suggested by Hinuma et al does not (Figure S6).These results indicate that the selection of band paths might miss the significant phonon modes for the NTE behavior.K. Here, we adopted the equilibrium volumes at the relevant finite temperatures obtained by minimizing free energy via QHA.In the finite temperatures, the phonons in the range between 0 and 7 THz are mainly excited, which have negative mode-Grüneisen parameters (Figure 5a).

Figure S1 .
Figure S1.Antiferromagnetic configuration of Cu2P2O7 visualized with (a) the primitive cell and (b) the conventional cell.Yellow arrows indicate the direction of the magnetic moment in Cu.

Figure S2 .
Figure S2.Comparison of calculated and experimental anisotropic atomic displacement parameters of oxygen atoms at 200 K (a) without and (b) with +U correction.
], whereas the other four C2/c phases of Sc2Zr2O7, Cu2Ta2O7, Zn2Ta2O7, and Mg2Ta2O7 are found to have imaginary phonon modes, indicating that they are dynamically unstable [FigureS4d-g].

Figure S5 .
Figure S5.Coordination number dependences of ionic radii for cations A, B and O.The dashed lines are the least-squares linear fittings, and R 2 values denote their coefficients of determination.

Figure S6 .
Figure S6.(a) Phonon band structure with respect to the values of the mode-Grüneisen parameters for Cu2P2O7.Blue and red indicate negative and positive values of mode-Grüneisen parameters, respectively.(b) Mode-Grüneisen parameters along the path of first-Brillouin zone of C2/c space group, which are colored according to the calculated phonon frequency.(c) First Brillouin zone shape and k-paths generated from the Hinuma et al [S22].

Figure S7 .
Figure S7.Schematic of Cu2P2O7 with calculated thermal ellipsoids at 200 K representing a 98% probability of containing the relevant atoms above the b-axis.PO4 tetrahedral units and CuO4 quadrilaterals are illustrated.

Figure S8 .
Figure S8.(a) Calculated total phonon densities of states and (b) the Bose-Einstein distributions at 100 K, 200 K, and 300 K.

Table S2 .
Valence electrons and cutoff radii of the PAW data sets used in the present study.