# Lattice dynamics of palladium in the presence of electronic correlations

###### Abstract

We compute the phonon dispersion, density of states, and the Grüneisen parameters of bulk palladium in the combined Density Functional Theory and Dynamical Mean-Field Theory. When electronic correlations are taken into account we find a good agreement with the experimentally measured phonon spectra, and ground state properties (equilibrium lattice parameter and bulk modulus) in comparison with density functional theory results. At the same time we demonstrate that the previously predicted softening of the phonon mode along the direction is absent in the presence of electronic correlations.

## I Introduction

The vibrational spectrum is a key quantity for understanding the interatomic interactions in a material. For solids, various inelastic scattering techniques exist to measure the phonon dispersion relation. These results can be readily compared with calculations within Density Functional Theory (DFT) ho.ko.64; kohn.99; jo.gu.89; jone.15. Within the DFT, a frequently used method to compute the phonon dispersion relation is the frozen phonon approach based on total energy calculations. The same formalism is also used to compute ground state properties such as equilibrium lattice parameters and bulk modulus. In the last decade it has been demonstrated that within the combined DFT and Dynamical Mean-Field Theory ge.ko.96; go.ko.92; me.vo.89; ko.vo.04, the so called LDA+DMFT approach held.07; ko.sa.06, many of the electronic ground state properties of -metal elements, their alloys and compounds can be properly described ko.sa.06; held.07; ka.ir.08; ku.le.17. This is also true in the case of palladium for which the LDA calculation of the bulk modulus turns out to overestimate the experimental value by more than . Indeed, using LDA+DMFT and employing a fluctuation exchange (FLEX) approximation as impurity solver in the Matsubara domain at K, with a local Coulomb interaction eV, and a local Hund’s rule coupling eV it was shown oe.ap.16 that the experimental volume a.u. is reproduced and the deviation in the bulk modulus from the experimental value is significantly reduced (). In addition, the results for the spectral function demonstrate the existence of a high-energy satellite formation in agreement with experiment ch.81. In our previous study we also discussed the importance of local and nonlocal correlation effects oe.ap.16 by comparing the LDA+DMFT results with the self-consistent quasiparticle GW approach sc.06; ko.07. More recently, a LDA+DMFT study sa.re.18 using a quantum Monte Carlo impurity solver found comparable electronic effective mass as in Ref. oe.ap.16. In addition in Ref. sa.re.18 employing a lattice (-dependent) FLEX solver resulted in a self-energy with weak -dependence. These results show the importance of including local, dynamic correlation effects in computing the ground state properties of Pd.

The goal of this paper is to calculate the changes in the phonon dispersions and phonon density of states of elemental solid palladium under pressure when electronic correlation effects are included. In addition, we investigate the volume dependence of the phonon frequencies, i.e. the mode Grüneisen parameters, since they provide information on the crystal anharmonicity. We employ the quasi-harmonic approximation (QHA), according to which the harmonic approximation holds at every volume, and the harmonic frequencies are replaced by renormalized volume dependent frequencies do.93; ho.74. It has been shown in the case of Pd that the QHA describes thermal dilation effects accurately and that higher order corrections are three orders of magnitude smaller for temperatures up to the melting temperature zo2.90.

Early experimental studies of the phonon spectra in palladium were performed by
Miiller and Brockhouse mii.bro.68; mii.bro.71; mi.75 using neutron scattering.
Results were reported for room temperature (K) and a Kohn anomaly
at a vector around was proposed.
The analysis by Freeman et al. fr.wa.79 showed a pronounced
enhancement of the
generalized susceptibility (Lindhard function) at a -vector of
around .
Later, Savrasov and Savrasov sa.sa.96 and Takezawa ta.na.05
reported numerical studies and found no signature of phonon anomalies in Pd.
Stewart st.08 reported a softening in the phonon dispersion at around
using DFT, which still somewhat
underestimated the experimental value. He showed that the electronic degrees
of freedom are responsible for the softening of the phonon mode.
However, Liu et al. li.ya.11 pointed out that the observed
Kohn anomaly by Stewart depends sensitively on the technical details of the
simulation. They observed that the Kohn anomaly vanishes when a
calculation is performed with a denser -mesh.
Recently, Dal Corso da.13 investigated the influence of various
exchange correlation functionals on the phonon spectrum in transition and noble
metals including the zero-point anharmonic expansion
(ZPAE) ^{1}^{1}1see for example chapter 7 Ref. gr.99 and thermal expansion of the lattice constant, which also affects the equilibrium bulk modulus gr.hi.07; st.04.
He showed that the LDA calculations overestimate the experimentally measured phonon frequencies by while the generalized gradient approximation (GGA) pbe1 underestimates by about at the high symmetry .

Until now, the numerical modeling of phonons in Pd did not include correlation effects. Following up our previous work, in this article we compare the results of the lattice dynamics obtained by LDA and LDA+DMFT. In Sec. II we discuss the details of the calculation and in Sec. III we briefly describe the methodology of frozen phonon calculations within the LDA and LDA+DMFT method. In Sec. III.1 the phonon dispersion curves are computed within the Born-Oppenheimer (BO) adiabatic approximation including DMFT corrections for the electrons by assuming the system to be perfectly harmonic. In Sec. III.2 we calculate the mode Grüneisen parameters which can be used to assess the anharmonicity within the QHA. This allows us to estimate the lineshifts originating from the thermal expansion and pressure zo.90; ma.62; ma.62_1. The article is closed with a brief conclusion in Sec. IV.

## Ii Computational Details

Pd crystallizes in the face-centered cubic structure ( space group) with a one-atom basis. The palladium atom is located at the Wyckoff position , i.e. . We performed a lattice relaxation within LDA with the parametrization of Perdew and Wang pe.wa.92 for the exchange-correlation functional. In the present study we employ the full-potential linearized muffin-tin orbitals (FPLMTO) method, as implemented in the RSPt code rsptbook; gr.ma.12. We used three kinetic energy tails with corresponding energies , , and Ry. The -mesh grid is for the equations of state. The muffin-tin radius was set to a.u. and was kept constant for all unit-cell volumes. For the charge density and the angular decomposition of the potential inside the muffin-tin spheres a maximum angular momentum was set. The calculations include spin-orbit coupling and scalar-relativistic terms. The so-called thermal broadening with K using the Fermi-Dirac function as integration kernel was employed for the Brillouin-zone sampling. To go beyond the LDA and include electronic correlations the local interaction term is added to the Hamiltonian, where . Here the operators annihilate (create) electrons on the orbital with the spin -projection . The summation over the set of basis functions () is restricted to the correlated subset (the -orbital subspace in our case). The interaction matrix elements are components of a rotational invariant fourth-rank tensor. These are computed from the on-site Slater-Koster integrals using atomic wave functions co.du.16. The averaged local Hubbard interaction is varied in the range from eV to eV, while the Hund’s rule coupling is kept fixed at eV. A double-counting correction is included in order to cancel the contributions due to the total energy functional originating from the electron-electron interaction captured within the exchange-correlation functional of DFT. This correction is applied to the self-energy gr.ma.12: with , since this is known to be the appropriate correction for metals, which ensures Fermi-liquid behavior pe.ma.03. We note that the frozen-phonon approach is not limited to a specific choice of double-counting. The Spin-polarized T-matrix Fluctuation Exchange Approximation (Sp-TFLEX) li.ka.97 on the Matsubara domain, with imaginary frequencies, was used as an impurity solver for the DMFT problem. We chose the same temperature (K) for the imaginary-frequency Matsubara mesh in DMFT as for the Brillouin-zone integration. The LDA+DMFT results do not change significantly in the temperature range from 275 to 325 K, which allows for a comparison between our results (K) and the experimental data measured at room temperature.

For the frozen phonon superlattice we used the Brillouin-zone integration by Fourier quadrature, which permits us to keep the -mesh density (namely -points per fcc Brillouin zone) consistent in all calculations fr.89. This minimizes systematic numerical errors when total energies are compared. We estimated errors due to -mesh sampling for the LDA study by computing all phonon frequencies using a denser -mesh of within the fcc Brillouin-zone for the frozen phonon calculations. We found that they depend on the considered -point and branch index, but remain below meV. Furthermore, the statistical error bar, estimated from the curve fit of the total energies as a function of phonon mode amplitude, has been found to be less than meV. In the presented Figures, we indicate only the systematic error-bars since they are the dominant source of error. The same error estimates have been performed for the LDA+DMFT calculations for the zone-boundary phonon frequencies. The LDA+DMFT error estimates are found to be smaller than their corresponding LDA values. A fine -mesh grid of is used to compute the phonon density of states (DOS) employing the tetrahedron integration method la.vi.84 after fitting to a five nearest-neighbor Born-von Karman (BvK) force constant model.

## Iii Phonon calculations from first principles

The frozen phonon approach ku.81; yi.82 is frequently used to calculate the phonon spectra. In this method the eigenvector of the given phonon is obtained from symmetry considerations when the appropriate atomic displacements are frozen in. The associated lowering of the symmetry corresponds to supercells of dimensions that match the phonon wavelength. In this way the method is not limited to the linear response of the energy functional with respect to lattice vibrations. The shortcoming of the frozen-phonon method is the supercell size, which may become impractically large, especially when long-wavelength phonons are considered. This difficulty can be circumvented by employing Density Functional Perturbation Theory (DFPT) (see Ref. ba.01 and references therein) in which the Hellmann-Feynman theorem is used hel.39; fe.39. In contrast to the frozen-phonon approach, DFPT enables the calculation of vibrational frequencies and eigenmodes at arbitrary wave vectors. This improvement in the computations is possible since the exchange and correlation kernel in DFPT can be parametrized and it takes a particularly simple form in the LDA ba.01; fel.17. There are more sophisticated versions of the DFPT which also include a Hartree-Fock mean-field correction (DFT+U) fl.gi.11, but they are not considered in the present study.

The linear response method in combination with LDA+DMFT was first applied on the Mott insulators, NiO and MnO sa.ko.03, and is currently being further extended le.an.14; ha.pa.16; ha.18. In combination with the many-body (LDA/GGA)+DMFT approach the frozen phonon technique was used to compute the phonon spectrum of the elemental metals, namely Fe (Ref. le.po.12) and Pu (Ref. da.sa.03), as well as the compound KCuF (Ref. le.bi.08). Below we take a similar approach and present results for the phonon band structure and the corresponding phonon density of states of palladium using the LDA and LDA+DMFT approach.

### iii.1 Phonon dispersions and density of states

In the BO approximation the nuclei are considered immobile during the characteristic electronic time scales. Within the harmonic approximation the potential energy surface is expanded up to second order in the ionic displacement. The equation of motion for the ionic subsystem is determined by the BO energy surface, obtained from the total energy within the LDA and LDA+DMFT. Its solution specifies the eigenmodes of a vibration with the characteristic frequencies , where the index enumerates the three linear independent phonon branches. The relation between the symmetry adapted basis and the Cartesian basis is explained in Appendix A.

method | (K) | (a.u.) | (GPa) | (GPa) | Refs | |||||
---|---|---|---|---|---|---|---|---|---|---|

LDA | 0 | 95.94 | 0 | 226.6 | T/S | [oe.ap.16] | 30.2 | 13.6 | 30.0 | 20.7 |

LDA | 0 | 99.30 | -10.3 | 186.5 | T/S | [oe.ap.16] | 27.9 | 12.8 | 27.7 | 19.3 |

LDA+DMFT | 95.94 | 4.5 | 225.7 | S | [oe.ap.16] | 30.9 | 13.7 | 30.3 | 20.8 | |

LDA+DMFT | 99.30 | -6.0 | 186.0 | S | [oe.ap.16] | 28.5 | 13.0 | 28.1 | 19.4 | |

Exp | 296 | 99.30 | 0 | 195 | S | [mii.bro.68,si.79] | 28.2 | 13.0 | 27.5 | 18.8 |

Exp | 298 | 99.30 | 0 | 189 | T | [gs.64] | - | - | - | - |

Exp | 296 | 99.30 | 0 | 193 | T | [so.93] | - | - | - | - |

The values in the column correspond either to the adiabatic () value, the isothermal () value, or to both when and coincide in the formalism. The difference between the reciprocal values (Ref. mii.bro.68) and (Ref. gs.64) is , where is the heat capacity and constant pressure, and is the volumetric thermal expansion coefficient. The bulk modulus is computed from the elastic constants (see Table I in Ref. si.79).

In Figs. 1a and 1c we present the phonon dispersion curves computed at the LDA equilibrium volume and at the experimental equilibrium volume, respectively, in comparison with the experimental measurements mii.bro.68. We find an good qualitative agreement with the experimentally measured dispersion curves. For the LDA equilibrium volume () the phonon frequencies are overestimated for most of the -points. This overestimation can be corrected considering the experimental equilibrium volume which is larger than oe.ap.16. Hence the phonon frequencies computed with correspond to the values which would be obtained for a negative pressure of GPa (see also Table 1). We find that by including electronic interactions of eV and does not improve the LDA results further. The bulk modulus and lattice constant computed within LDA+DMFT, however, show a better agreement with experiment than the LDA result oe.ap.16. One significant difference is observed along the -direction (see Figs. 1a and 1c), namely the mode softening obtained within the LDA at in the branch (see Appendix A for definition of the branch labeling) is reduced within LDA+DMFT. A similar reduction is also found in the experiment where increasing temperature (in the range of K up to K) has the effect of smoothening out the alleged Kohn anomaly, which is in agreement with the presented LDA+DMFT result mi.75. Whether this is a result of temperature alone or a combined effect of temperature and correlations, cannot be decided at this moment.

In Figs. 1b and 1d we present the phonon DOS computed at the LDA equilibrium volume and at the experimental equilibrium volume, respectively, and compare the LDA and LDA+DMFT phonon DOS with the experimental results. The main features of the experimental DOS are the two peaks at around 20 and 25 meV. At a.u. we find that both within the LDA and LDA+DMFT the peak positions are shifted up compared to the experimental phonon DOS. At a.u. the lower peak (around 20 meV) is captured by both the LDA and LDA+DMFT method, while the higher peak (around 25 meV) is shifted up in both methods. In LDA+DMFT the peak position (around 25 meV) in the DOS agrees well with the experiment, since it is only shifted by about meV.

In Fig. 2 we show the comparison of phonon spectra computed for different values of the local Hubbard interaction . For eV and eV the overall fit to the experimental phonon dispersion is good and a similar deviation is found as in the case of the LDA. For eV and eV we find an increase of the phonon frequencies relative to the LDA result. The softening of the phonon dispersion along the -line () is not seen in the case of LDA+DMFT for all investigated interaction strengths. For eV we find that the dispersion curve around even has the opposite curvature as compared to the experimentally measured phonon dispersion.

### iii.2 Phonon dispersion at different volumes

method | (K) | (K) | ||||||
---|---|---|---|---|---|---|---|---|

LDA: | 2.210 | 2.216 | 2.217 | |||||

LDA+DMFT: | 2.019 | 2.020 | 2.020 | |||||

Difference |

In the first attempt to analyze the volume dependence of the phonon frequencies on a quantitative level, it is useful to consider the Grüneisen parameter . It represents the relative change of with isotropic change in volume for a particular phonon branch at the symmetry point , and is computed according to:

(1) |

Usually the values of are positive and lie in the range gr.99, where a larger value suggests that anharmonic effects are important for the particular phonon mode of interest. Using the values provided by Eq. 1, we computed the average Grüneisen parameter defined as . We also computed the thermodynamic Grüneisen parameter which within the QHA takes the form: . The specific heat per mode has the form , with . We find that is weakly temperature dependent and essentially equivalent to the average value (see Tab. 2). The main effect of electronic correlations is a decrease of by about . This shows that the electronic correlations have the tendency to reduce the anharmonicity of the lattice degrees of freedom for typical phonon modes, i.e. the lattice becomes stiffer with increasing electronic correlations. One reaches a similar conclusion by analyzing the bulk modulus data.

In Fig. 3 we present the phonon frequencies at the zone boundaries as a function of the unit-cell volume in LDA and LDA+DMFT. Red and blue decorated lines correspond to transversal and longitudinal modes, respectively. Upon lattice expansion (increasing volume) phonon modes are seen to become softer. Note that in the studied range the phonon frequencies decrease almost linearly with the volume, and within a good approximation the linear curves are shifted rigidly by an amount (see also Tab. 2). These results indicate that the correlation effects lead to a stiffening of the lattice. As seen in Fig. 3, the frequency shifts due to volume expansion are significant in the considered volume range in comparison with the error estimates for the -modes, while this is not the case for the -mode.

In Tab. 2 we present the values of the for the modes , and in LDA and LDA+DMFT. The results show that is larger for the longitudinal modes and than for the transversal modes and . The overall effect of electronic correlations is the significant reduction of for the longitudinal modes, while for the transversal modes remains almost unchanged (-mode) or even increases (-mode).

## Iv Conclusion

We performed LDA and LDA+DMFT calculations and determined the phonon dispersion , the phonon density of states, the thermodynamic Grüneisen parameter , and the mode Grüneisen parameters for zone-boundary phonon modes. The phonon dispersion was computed using the frozen-phonon method which relies on reliable calculation of total energies and on the Born-Oppenheimer and the harmonic approximations. In our LDA+DMFT calculations the static part of the electronic response to lattice oscillations is included, which reduces the phonon-softening. The change in the thermodynamic Grüneisen parameter, , as a function of the strength of the Coulomb and exchange interaction parameters, demonstrates that in the presence of electron correlations anharmonic effects become less important.

## Acknowledgement

We are grateful to Ferdi Aryasetiawan, Igor Di Marco, Milos Radonjić and Weiwei Sun for stimulating discussions. WHA thanks Levente Vitos for the kind hospitality at KTH Stockholm. Part of the computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at Beskow. Financial support offered by the Augsburg Center for Innovative Technologies, and by the Deutsche Forschungsgemeinschaft (through TRR 80/F6) is gratefully acknowledged. IL is grateful for the financial support provided by the Russian Science Foundation (project Nr. 18-12-00492).

## Appendix A Symmetry Adapted Dynamical Matrix

In order to study lattice vibrations in the harmonic approximation the potential energy surface is expanded in second order in the ionic displacement

where are the deviations of the ionic positions from their equilibrium positions, and denote the components of a vector in the 3D Cartesian space, is the BO energy surface, and is the matrix of interatomic force constants in the Cartesian representation. We introduce the dynamical matrix

(3) |

were is the ionic mass. The Hamiltonian of the effective phonon system is quadratic in the displacement and anharmonic terms are neglected. The normal modes of the phonon system are determined by solving the following equation of motion:

(4) |

where are components of the normalized vector in 3D Cartesian vector space. Non-trivial solutions are obtained by solving the secular equation . The solutions are denoted as and the corresponding eigenmodes are given by the nullspace of the matrix :

(5) |

The index enumerates the 3 linear independent phonon branches in Pd. The solution of the secular equation can be simplified by making use of the symmetry of the system. It turns out that the eigenvalue problem in Eq. 4 can be decomposed into one-dimensional blocks along the certain high symmetry lines in the irreducible Brillouin zone. This allows us to determine the normal modes by applying group theoretical techniques without explicitly obtaining the solution of the secular equation.

The phonon dispersion curves for a five nearest neighbor force-constant model is shown in Fig 4. The irreducible representations (irreps) are indicated in the figure. When the dynamical matrix is written in the symmetry adapted basis then the eigenvalue problem in Eq. 4 decomposes into blocks according to the generalized Wigner-Eckart theorem lu.12; wi.93:

here we introduced the multi-index (also the branch index in our case) which runs over all irreps , the row of the irrep and the the multiplicity of an irrep in the basis . The unitary transformation from the Cartesian to the symmetry adapted basis is denoted as . The symmetry-adapted basis is obtained employing the projecting technique of Ref. lu.12. Due to the highly symmetric space-group and the simple lattice basis, consisting only of a single atom per unit cell, the dynamical matrix decomposes into blocks of dimension . Therefore the phonon frequencies can be computed efficiently in the frozen phonon approach since only a single symmetry mode needs to be set up in order to determine the mode frequency. The supercell of the subgroup has to be commensurate with the wave-length and ions located at low symmetric Wyckoff positions are displaced according to the symmetry mode of interest.

The construction of the symmetry adapted basis functions is a cumbersome task since one needs to perform the projection of the basis onto space group irreps for various -points in the irreducible Brillouin zone. It is useful to interpret the symmetry modes as order parameter directions of a structural second-order phase transition, and the potential energy as a function of the ionic displacement as the Landau free energy expansion landau1968statistical; landau1937theory. The advantage of this approach is that the order parameters of a phase transition can be found without directly referring to the atomic positions and instead focusing on specific irreps and their invariant subspaces. The invariant subspaces determine subgroups, the so called isotropy groups, of the high symmetric parent structure.

Stokes and Hatch have developed a software package which allows one to obtain information about the normal modes of oscillations in a crystal by employing the concept of isotropy subgroups st.ha.84; ha.sto.87; st.hat.91; st.ha.88; st.95. By listing the number of irreps of the space-group at a given -point one can determine the symmetry modes systematically. Using one irrep at a time we project out modes which transform like basis functions of that irrep st.hat.91; st.95.