Lattice dynamics of palladium in the presence of electronic correlations

TL;DR

This study combines DFT and DMFT to accurately compute phonon properties of palladium, revealing that electronic correlations cause low-temperature phonon mode softening near the Kohn anomaly.

Contribution

It introduces a combined DFT+DMFT approach to study lattice dynamics in palladium, highlighting the impact of electronic correlations on phonon behavior at low temperatures.

Findings

Good agreement with experimental phonon spectra and ground state properties.

Phonon frequency near the Kohn anomaly decreases significantly below 20 K.

Correlation effects weaken ionic restoring forces, causing mode softening.

Abstract

We compute the phonon dispersion, density of states, and the Gr\"uneisen parameters of bulk palladium in the combined density functional theory (DFT) and dynamical mean-field theory (DMFT). We find good agreement with experimental results for ground state properties (equilibrium lattice parameter and bulk modulus) and the experimentally measured phonon spectra. We demonstrate that at temperatures $T \lesssim 20~K$ the phonon frequency in the vicinity of the Kohn anomaly, $\omega_{T1}({\bf q}_{K})$, strongly decreases. This is in contrast to DFT where this frequency remains essentially constant in the whole temperature range. Apparently correlation effects reduce the restoring force of the ionic displacements at low temperatures, leading to a mode softening.

Tables

Table 1. Table 1: Left part: Volume V 𝑉 V for various methods (upper part) and experiments (lower part). Middle part: Computed or experimentally measured pressure P 𝑃 P , and isothermal/adiabatic bulk modulus K x = T / S subscript 𝐾 𝑥 𝑇 𝑆 K_{x=T/S} in the upper part and lower part, respectively. The computed values for the bulk modulus K x subscript 𝐾 𝑥 K_{x} includes a correction term for the pressurized sample. The computations where done at T = 316 𝑇 316 T=316\, K and the experiments have been performed at room temperature as well. Right part: Zone boundary frequencies computed under adiabatic conditions are compared with the experimental phonon frequencies  Miiller and Brockhouse ( 1968 ) .

method	$V$(a.u.3)	$P$(GPa)	$K_x$(GPa)a	$x$	Refs	${\omega }_{\nu }\left(𝐪\right)$
						$(L,L_2^{-})$	$(L,L_3^{-})$	$(X,X_3^{-})$	$(X,X_5^{-})$
LDA	95.94	0	226.6	T/S	[Östlin et al., 2016]	30.2	13.6	30.0	20.7
LDA	99.30	-10.3	186.5	T/S	[Östlin et al., 2016]	27.9	12.8	27.7	19.3
LDA+DMFT	95.94	4.5	225.7	S	[Östlin et al., 2016]	30.9	13.7	30.3	20.8
LDA+DMFT	99.30	-6.0	186.0	S	[Östlin et al., 2016]	28.5	13.0	28.1	19.4
Exp	99.30	0	195	S	[Miiller and Brockhouse, 1968,Singh and Hemkar, 1979] b	28.2	13.0	27.5	18.8
Exp	99.30	0	189	T	[Gschneidner Jr, 1964]	-	-	-	-
Exp	99.30	0	193	T	[Söderlind et al., 1993]	-	-	-	-

Table 2. Table 2: Left part: Grüneisen parameter γ ν G ​ ( 𝐪 ) superscript subscript 𝛾 𝜈 𝐺 𝐪 \gamma_{\nu}^{G}({\mathbf{q}}) computed in LDA and LDA+DMFT, the difference between frequencies calculated in LDA and LDA+DMFT δ ​ ω ν ​ ( 𝐪 ) = ω ν ​ ( 𝐪 ) DMFT − ω ν ​ ( 𝐪 ) LDA 𝛿 subscript 𝜔 𝜈 𝐪 subscript 𝜔 𝜈 superscript 𝐪 DMFT subscript 𝜔 𝜈 superscript 𝐪 LDA \delta\omega_{\nu}(\mathbf{q})=\omega_{\nu}(\mathbf{q})^{\mathrm{DMFT}}-\omega_{\nu}(\mathbf{q})^{\mathrm{LDA}} computed for the branches ν = L 2 − , L 3 − , X 3 − , X 5 − 𝜈 superscript subscript 𝐿 2 superscript subscript 𝐿 3 superscript subscript 𝑋 3 superscript subscript 𝑋 5 \nu=L_{2}^{-},L_{3}^{-},X_{3}^{-},X_{5}^{-} . Right part: Average Grüneisen parameter ⟨ γ G ⟩ delimited-⟨⟩ superscript 𝛾 𝐺 \langle\gamma^{G}\rangle , and thermodynamic Grüneisen parameter γ t ​ d G superscript subscript 𝛾 𝑡 𝑑 𝐺 \gamma_{td}^{G} computed in both methods.

method		$(L,L_2^{-})$	$(L,L_3^{-})$	$(X,X_3^{-})$	$(X,X_5^{-})$	$⟨{\gamma }^G\left(𝐪\right)⟩$	${\gamma }_{td}^G$ ($T=298$K)	${\gamma }_{td}^G$ ($T=316$K)
LDA:	${\gamma }_{\nu }^G\left(𝐪\right)$	$2.34$	$1.68$	$2.29$	$2.01$	2.210	2.216	2.217
LDA+DMFT:	${\gamma }_{\nu }^G\left(𝐪\right)$	$2.28$	$1.74$	$2.24$	$2.01$	2.019	2.020	2.020
Difference	$\delta {\omega }_{\nu }\left(𝐪\right)$	$0.45$	$-0.15$	$0.39$	$0.13$

Equations

$${\gamma }_{\nu }^G(\mathbf{q})$$

$$U_p(\{\mathbf{R}\})$$

$$D_{{\alpha }_1{\alpha }_2}(\mathbf{q})=\frac{1}{M}\sum _{\mathbf{R}}{C}_{{\alpha }_1{\alpha }_2}(\mathbf{R})e^{i\mathbf{q}\cdot \mathbf{R}},$$

$${\omega }^2{e}_{{\alpha }_1}(\mathbf{q})$$

$${\omega }_{\nu }^2(\mathbf{q})e_{{\alpha }_1,\nu }(\mathbf{q})$$

$$D_{{\alpha }_1{\alpha }_2}(\mathbf{q})$$

Full text

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) Hohenberg and Kohn (1964); Kohn (1999); Jones and Gunnarsson (1989); Jones (2015). 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 (DMFT) Metzner and Vollhardt (1989); Georges and Kotliar (1992); Georges et al. (1996); Kotliar and Vollhardt (2004), the so called DFT+DMFT approach Held (2007); Kotliar et al. (2006), many of the electronic ground state properties of $d$-metal elements, their alloys and compounds can be properly described Kotliar et al. (2006); Held (2007); Katsnelson et al. (2008); Kuneš et al. (2017). This is also true in the case of palladium where results of the local density approximation (LDA) for the bulk modulus turns out to overestimate the experimental value by more than $10\%$. Indeed, using LDA+DMFT it was shown Östlin et al. (2016) that the experimental volume $V=99.3\,$a.u.3 is reproduced and that the difference in the bulk modulus from the experimental value is reduced to less than $1\,\%$. In addition, the results for the spectral function demonstrate the existence of a high-energy satellite formation in agreement with experiment Chandesris et al. (1981). In our previous study we also discussed the importance of local and nonlocal correlation effects Östlin et al. (2016) by comparing the LDA+DMFT results with the self-consistent quasiparticle GW approach van Schilfgaarde et al. (2006); Kotani et al. (2007). More recently, a LDA+DMFT study Savrasov et al. (2018) using a quantum Monte Carlo impurity solver found comparable electronic effective mass as in Ref. Östlin et al., 2016. In addition in Ref. Savrasov et al., 2018 employing a lattice (${\bf k}$-dependent) FLEX solver resulted in a self-energy with weak ${\bf k}$-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. 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 Dove (1993); Horton and Maradudin (1974). 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 Zoli (1990).

Early experimental studies of the phonon spectra in palladium were performed by Miiller and Brockhouse Miiller and Brockhouse (1968, 1971); Miiller (1975) using neutron scattering, and a Kohn anomaly at a vector around ${\mathbf{q}}=\frac{2\pi}{a}[0.35,0.35,0]$ was proposed. The analysis by Freeman et al. Freeman et al. (1979) showed a pronounced enhancement of the generalized susceptibility (Lindhard function) at a $\mathbf{q}$-vector of around $\frac{2\pi}{a}[0.325,0.325,0]$. Later, Savrasov and Savrasov Savrasov and Savrasov (1996) and Takezawa Takezawa et al. (2005) reported numerical studies where they found no signature of phonon anomalies in Pd. Stewart Stewart (2008) reported a softening in the phonon dispersion at around ${\mathbf{q}}=\frac{2\pi}{a}[0.3,0.3,0]$ 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. Liu et al. (2011) 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 $\mathbf{k}$-mesh. Recently, Dal Corso Dal Corso (2013) investigated the influence of various exchange correlation functionals on the phonon spectrum in transition and noble metals. This included the zero-point anharmonic expansion (ZPAE) 111see for example chapter 7 Ref. Grimvall, 1999 and the thermal expansion of the lattice constant, which also affects the equilibrium bulk modulus Grabowski et al. (2007); Staroverov et al. (2004). He showed that the LDA calculations overestimate the experimentally measured phonon frequencies by $5\,\%$ while the generalized gradient approximation (GGA) Perdew et al. (1996) underestimates them by about $7\,\%$ at the high symmetry point $L=\frac{\pi}{a}[1,1,1]$.

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 employed to assess the anharmonicity on the level of the QHA. Namely, although the QHA neglects cubic or quartic terms in the expansion of the Born-Oppenheimer energy surface, it includes the anharmonic strain- or volume-dependence of the frequencies. Sec. III.3 addresses the finite temperature behaviour of a specific mode close to the Kohn anomaly. The phonon frequencies computed within LDA+DMFT show a strong temperature dependence which is not observed in DFT(LDA). The article is closed with a brief conclusion in Sec. IV.

II Computational Details

Pd crystallizes in the face-centered cubic structure ($Fm\overline{3}m$ space group) with a one-atom basis. The palladium atom is located at the Wyckoff position $4a$, i.e. $(0,0,0)$. We performed a lattice relaxation within LDA with the parametrization of Perdew and Wang Perdew and Wang (1992) 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 Wills et al. (2010); Grånäs et al. (2012); Di Marco et al. (2009). We used three kinetic energy tails with corresponding energies $-0.3$, $-2.3$, and $-1.5\,$Ry. The muffin-tin radius was set to $2.45\,$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 $l_{\text{max}}=8$ was set. The calculations include spin-orbit coupling and scalar-relativistic terms. Brillouin-zone integration was performed using a Fermi-Dirac function (thermal broadening). The ${\mathbf{k}}$-mesh grid is $16\times 16\times 16$ for the equations of state. To go beyond the LDA and include electronic correlations the local interaction term $1/2\sum_{\{\lambda_{i}\}}U_{\lambda_{1}\lambda_{2}\lambda_{4}\lambda_{3}}c^{\dagger}_{\lambda_{1}}c^{\dagger}_{\lambda_{2}}c_{\lambda_{3}}c_{\lambda_{4}}$ is added to the Hamiltonian, where $\lambda_{i}=(m_{i},\sigma_{i})$. Here the operators $c^{(\dagger)}_{\lambda_{i}}$ annihilate (create) electrons on the orbital $m_{i}$ with the spin $z$-projection $\sigma_{i}=\uparrow,\downarrow$. The summation over the set of basis functions $\{m_{i}=-l,\dots,l\}$ ($l=2$) is restricted to the correlated subset (the $4d$-orbital subspace in our case). We used the so-called muffin-tin-only correlated subspace which corresponds to basis functions which are zero outside the muffin tin region; see Sec. II in Ref. Grechnev et al., 2007 and Ref. Grånäs et al., 2012 for additional information about the basis-function dependence. The interaction matrix elements $U_{\lambda_{1}\lambda_{2}\lambda_{4}\lambda_{3}}$ are components of a rotational invariant fourth-rank tensor. These are computed from the on-site Slater-Koster integrals using atomic wave functions Coury et al. (2016). The averaged local Hubbard interaction $U$ is varied in the range from $0.75\,$eV to $1.30$ eV, while the Hund’s rule coupling is kept fixed at $J=0.3\,$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ånäs et al. (2012): $\Sigma_{\lambda_{1}\lambda_{2}}(i\omega_{n})\to\Sigma_{\lambda_{1}\lambda_{2}}(i\omega_{n})-\Sigma_{\lambda_{1}\lambda_{2}}^{DC}(0)$ with $\Sigma_{\lambda_{1}\lambda_{2}}^{DC}(0)=\delta_{m_{1}m_{2}}/(2l+1)\sum_{m_{3}}\Sigma_{m_{3}\sigma_{1},m_{3}\sigma_{2}}(0)$, since this is known to be the appropriate correction for metals, which ensures Fermi-liquid behavior Petukhov et al. (2003). 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) Lichtenstein and Katsnelson (1998) on the Matsubara domain was used as an impurity solver for the DMFT problem. Throughout this paper finite temperatures are only considered for the electronic subsystem, where the temperature enters in the Matsubara frequencies $\omega_{n}=(2n+1)\pi T$, ($n\in\mathcal{Z}$). We chose the same temperature for the imaginary-frequency Matsubara mesh in DMFT as for the Brillouin-zone integration. The electronic charge was updated self-consistently in all LDA and LDA+DMFT calculations Grånäs et al. (2012). 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 ($T=316\,$K) and the experimental data measured at room temperature. In the search of a mode softening in the vicinity of the Kohn anomaly we also performed a finite temperature study in an extended temperature range down to $15\,$K.

For the frozen phonon superlattice we used the Brillouin-zone integration by Fourier quadrature, which permits us to keep the $\mathbf{k}$-mesh density (namely $16\times 16\times 16$ $\mathbf{k}$-points per fcc Brillouin zone) consistent in all calculations Froyen (1989). This minimizes systematic numerical errors when total energies are compared. We estimated errors due to ${\mathbf{k}}$-mesh sampling for the LDA study by computing all phonon frequencies using a denser ${\mathbf{k}}$-mesh of $32\times 32\times 32$ within the fcc Brillouin-zone for the frozen phonon calculations. The error estimates for the phonon frequencies were computed as the differences $\Delta_{\nu}(\mathbf{q})=|\omega_{\nu}(\mathbf{q})_{\text{fine}}-\omega_{\nu}(\mathbf{q})_{\text{sparse}}|$, where $\omega_{\nu}(\mathbf{q})_{\text{fine}}$ and $\omega_{\nu}(\mathbf{q})_{\text{sparse}}$ denote the phonon frequencies which are computed using the $32\times 32\times 32$ and $16\times 16\times 16$ ${\mathbf{k}}$-mesh, respectively. We found that $\Delta_{\nu}(\mathbf{q})$ is smaller than $3\,$meV. In the presented Figures, we indicate the errors due to finite k-mesh sampling. 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. The reduction of the sensitivity of the LDA+DMFT results as compared to the LDA results is expected to originate from the broadening of the one-electron spectra due to many-body effects. A fine ${\mathbf{q}}$-mesh grid of $24\times 24\times 24$ is used to compute the phonon density of states (DOS) employing the tetrahedron integration method Lambin and Vigneron (1984) after fitting to a five nearest-neighbor Born-von Karman (BvK) force constant model.

III Phonon calculations from first principles

The frozen phonon approach Kunc and Martin (1981); Yin and Cohen (1982) 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. A shortcoming of the frozen-phonon method is the supercell size, which may become impractically large, especially when long-wavelength phonons are considered.

The linear response method in combination with LDA+DMFT was first applied on the Mott insulators, NiO and MnO Savrasov and Kotliar (2003), and is currently being further extended Leonov et al. (2014); Haule and Pascut (2016); Haule (2018). 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. Leonov et al., 2012) and Pu (Ref. Dai et al., 2003), as well as the compound KCuF3 (Ref. Leonov et al., 2008). Below we take the frozen-phonon 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 $\omega_{\nu}(\mathbf{q})$, where the index $\nu=1,2,3$ enumerates the three linear independent phonon branches. The relation between the symmetry adapted basis $|\nu\rangle$ and the Cartesian basis is explained in Appendix A.

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 Miiller and Brockhouse (1968). We find an good qualitative agreement with the experimentally measured dispersion curves (see also Table 1). For the LDA equilibrium volume ($V_{\rm LDA}$) the phonon frequencies are overestimated for most of the ${\mathbf{q}}$-points. This overestimation can be corrected considering the experimental equilibrium volume $V_{\rm Exp}$ which is larger than $V_{\rm LDA}$ Östlin et al. (2016). Hence the phonon frequencies computed with $V_{\rm Exp}$ correspond to the values which would be obtained for a negative pressure of $P=-10.3\,$GPa (see also Table 1). We find that by including electronic interactions of $U=0.75\,$eV and $J=0.3\,$eV 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 Östlin et al. (2016). One significant difference is observed along the $[110]$-direction (see Figs. 1a and 1c), namely the mode softening obtained within the LDA at ${\bf q}=\frac{2\pi}{a}[0.32,0.32,0]$ in the $T_{1}$ branch (see Appendix A for definition of the branch labeling) is reduced within LDA+DMFT. Hence the presented result, computed at high temperatures, is not in line with the experiment. A detailed study on the temperature dependence of the mode in the vicinity of ${\bf q}=\frac{2\pi}{a}[0.32,0.32,0]$ will be presented in Sec. III.3. There it is shown that at lower temperatures ($T\leq 15\,$K) the mode softening is well described within LDA+DMFT.

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 $V=96.0\,$a.u.3 we find that both within the LDA and LDA+DMFT the peak positions are shifted up compared to the experimental phonon DOS. At $V=99.3\,$a.u.3 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 $1\,$meV.

In Fig. 2 we compare phonon spectra computed at $316\,$K for different values of the local Hubbard interaction $U$ with the experimental results measured at room temperature Miiller and Brockhouse (1968). For $U=0.75\,$eV and $J=0.3\,$eV the overall agreement with the experimental phonon dispersion is good, differences being similar to those in the case of the LDA results. For $U=1.0\,$eV and $U=1.3\,$eV at $T=316$ K we find that the phonon frequencies are higher than those obtained within the LDA. Furthermore, the phonon frequencies are found to increase with increasing Coulomb interaction strength, i.e., the modes become stiffer.

The softening of the phonon dispersion along the $\Sigma$-line ($\mathbf{q}_{\text{nest}}=\frac{2\pi}{a}[0.32,0.32,0]$) is not described by the LDA+DMFT results for any of the investigated interaction strengths. This will be further discussed in Sec. III.3, where we address the finite temperature behavior of the phonon frequencies.

III.2 Phonon dispersion at different volumes

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 $\gamma_{\nu}^{G}({\mathbf{q}})$. It represents the relative change of $\omega_{\nu}(\mathbf{q})$ with isotropic change in volume for a particular phonon branch $\nu$ at the symmetry point ${\mathbf{q}}$, and is computed according to:

[TABLE]

Usually the values of $\gamma_{\nu}^{G}({\mathbf{q}})$ are positive and lie in the range $1.5\pm 1.0$ (Ref. Grimvall, 1999), where a larger value suggests that anharmonic effects are important for the particular phonon mode of interest. This is only valid when anharmonicity is weak to start with. Using the values provided by Eq. 1, we computed the average Grüneisen parameter $\langle\gamma^{G}\rangle$ defined as $\langle\gamma^{G}\rangle=\sum_{\nu}1/3V_{BZ}\int d\mathbf{q}\gamma_{\nu}^{G}({\mathbf{q}})$. We also computed the thermodynamic Grüneisen parameter which within the QHA takes the form: $\gamma_{td}^{G}=\sum_{\mathbf{q},\nu}C_{V}(\mathbf{q},\nu)\gamma^{G}_{\nu}(\mathbf{q})/\sum_{\mathbf{q},\nu}C_{V}(\mathbf{q},\nu)$. The specific heat per mode has the form $C_{V}(\mathbf{q},\nu)=k_{B}x^{2}e^{x}/(e^{x}-1)^{2}$, with $x=\hbar\beta\omega_{\nu}(\mathbf{q})$. We find that $\gamma_{td}^{G}$ is essentially equivalent to the average value $\langle\gamma^{G}\rangle$ (see Tab. 2). The main effect of electronic correlations is a decrease of $\gamma_{td}^{G}$ by about $10\%$. We provide these results with the following interpretation: We know from previous studiesÖstlin et al. (2016) that electronic correlations lead to a softening of the lattice it can be expected that the coupling between phonon degrees of freedom and the volume are reduced. Hence the Grüneisen parameter is expected to be reduced by correlation and that is what is seen.

In Fig. 3 we present the phonon frequencies at the zone boundaries as a function of the unit-cell volume $V$ 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 $\delta\omega_{\nu}({\bf q}=X,L)=\omega_{\nu}(\mathbf{q})^{\mathrm{DMFT}}-\omega_{\nu}(\mathbf{q})^{\mathrm{LDA}}$ (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 $L_{2}^{-},X_{3}^{-},X_{5}^{-}$-modes, while this is not the case for the $L_{3}^{-}$-mode.

In Tab. 2 we present the values of the $\gamma_{\nu}^{G}({\mathbf{q}})$ for the modes $X_{5}^{-},L_{3}^{-},X_{3}^{-}$, and $L_{2}^{-}$ in LDA and LDA+DMFT. The results show that $\gamma_{\nu}^{G}({\mathbf{q}},\nu)$ is larger for the longitudinal modes $X_{3}^{-}$ and $L_{2}^{-}$ than for the transversal modes $L_{3}^{-}$ and $X_{5}^{-}$. The overall effect of electronic correlations is the significant reduction of $\gamma_{\nu}^{G}({\mathbf{q}},\nu)$ for the longitudinal modes, while for the transversal modes $\gamma_{\nu}^{G}({\mathbf{q}},\nu)$ remains almost unchanged ($X_{5}^{-}$-mode) or even increases ($L_{3}^{-}$-mode).

III.3 Phonon dispersion at different temperatures

In the previous sections we discussed LDA+DMFT results for phonon spectra calculated at $316\,$K. For the wave vector ${\mathbf{q}}=\frac{2\pi}{a}[0.375,0.375,0]$, which is close to the nesting vector $\mathbf{q}_{\text{nest}}$, the phonon modes were found to become stiffer with increasing Coulomb interaction strength.

In this section we present results of LDA+DMFT calculations performed in a wide temperature range for interaction parameters $U=1.00\,$eV and $J=0.30\,$eV. We used the harmonic approximation and, consequently, the renormalization of phonon frequencies is only due to the electronic correlations at finite temperatures while the anharmonic (phonon-phonon) terms are ignored.

Within the LDA it is easy to compute the phonon spectra at finite temperatures since this requires only a multiplication with the Fermi function when only the electronic temperature is considered, as is the case in our investigation. By contrast, for low temperature calculations in the LDA+DMFT approach a large number of Matsubara frequencies must be included in the impurity solver to reach convergence. This increases the computational time very significantly. The LDA+DMFT calculation of phonons were performed using both a fine and a sparse Matsubara mesh of 32768 and 4096 frequencies, respectively. We found that for temperatures $T$ larger than about $150\,$K the difference in the phonon frequencies computed with the finer Matsubara mesh becomes negligible in comparison with the sparse mesh. Hence, we will use the sparse mesh in the following to keep the computational effort at a manageable level.

In Fig. 4a we compare the experimental and the computed phonon dispersions along the direction $\Sigma$ close to room temperature. One clearly observes that in the vicinity of the possible Kohn anomaly at ${\bf q}_{nest}$ the DMFT results lie above the LDA results and the experimental data Miiller and Brockhouse (1968). To understand the frequency behavior near ${\bf q}_{nest}$ better we computed the temperature dependence down to low temperatures ($T=15\,$K), both within LDA and LDA+DMFT. The results are shown in Fig. 4b. For temperatures $T\leq 20\,$K the phonon frequencies obtained within LDA+DMFT are found to be significantly reduced compared to room temperature. Not only are they now lower than the LDA results but, in fact, they approach the experimental values. Clearly the LDA+DMFT approach performs better at lower temperatures where quantum effects are more relevant. The strong temperature dependence of the phonon frequencies at wave-vectors close to the nesting vector (Kohn anomaly) in the presence of electronic correlations may be interpreted as follows. When the temperature is lowered the screening of the atomic forces becomes more efficient, which leads to significant phonon softening. This effect cannot be captured within LDA. The remaining discrepancies may be partly attributed to phonon-phonon scattering not considered in the present calculations.

IV Conclusion

We performed LDA and LDA+DMFT calculations and determined the phonon dispersion $\omega_{\nu}(\mathbf{q})$, the phonon density of states, the thermodynamic Grüneisen parameter $\gamma_{td}^{G}$, and the mode Grüneisen parameters $\gamma_{\nu}^{G}(\mathbf{q})$ for zone-boundary phonon modes. The phonon dispersion was computed using the frozen-phonon method which relies on the Born-Oppenheimer and the harmonic approximations and requires accurate total energies.

The change in the thermodynamic Grüneisen parameter, $\gamma_{td}^{G}$, as a function of the strength of the Coulomb $U$ and exchange $J$ interaction parameters, demonstrates that in the presence of electron correlations anharmonic effects on the level of the QHA become less important.

We also studied the phonon dispersion at finite temperatures in the vicinity of the possible Kohn anomaly. In LDA+DMFT low temperature calculations require considerable computational effort. We find that at $T\sim 20$ K the phonon frequency obtained within LDA+DMFT has strongly decreased relative to the value at room temperature and approaches the experimental result. We interpret this softening as being a consequence of the correlation-induced reduction of the restoring force of the ionic displacement in the harmonic approximation.

Acknowledgement

We are grateful to Ferdi Aryasetiawan, Milos Radonjić and Weiwei Sun for stimulating discussions. WHA thanks Levente Vitos for the kind hospitality at KTH Stockholm. I.D.M. acknowledges support by the appointment to the JRG program at the APCTP through the Science and Technology Promotion Fund and Lottery Fund of the Korean Government, as well as support by the Korean Local Governments, Gyeongsangbuk-do Province and Pohang City. 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

[TABLE]

where $u_{\alpha}(\mathbf{R})=R_{\alpha}-R^{0}_{\alpha}$ are the deviations of the ionic positions from their equilibrium positions, and $\alpha_{1},\alpha_{2}=x,y,z$ denote the components of a vector in the 3D Cartesian space, $U_{p}(\{\mathbf{R}_{n}\})$ is the BO energy surface, and $C_{\alpha_{1}\alpha_{2}}(\mathbf{R}-\mathbf{R}^{\prime})$ is the matrix of interatomic force constants in the Cartesian representation. We introduce the dynamical matrix

[TABLE]

were $M$ is the ionic mass. The Hamiltonian of the effective phonon system $H_{p}=-\sum_{{\mathbf{R}}}\frac{1}{2M}\nabla^{2}_{{\mathbf{R}}}+U_{p}(\{\mathbf{R}\})$ 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:

[TABLE]

where $e_{\alpha}(\mathbf{q})$ are components of the normalized vector in 3D Cartesian vector space. Non-trivial solutions are obtained by solving the secular equation ${\rm det}\left(D_{\alpha_{1}\alpha_{2}}({\mathbf{q}})-\delta_{\alpha_{1}\alpha_{2}}\omega^{2}\right)=0$. The solutions are denoted as $\omega_{\nu}^{2}({\mathbf{q}})$ and the corresponding eigenmodes $e_{\alpha,\nu}(\mathbf{q})$ are given by the nullspace of the matrix $D_{\alpha_{1},\alpha_{2}}(\mathbf{q})-\omega_{\nu}^{2}({\mathbf{q}})\delta_{\alpha_{1},\alpha_{2}}$:

[TABLE]

The index $\nu=1,2,3$ 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 $e_{\alpha,\nu}(\mathbf{q})$ 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 5. The irreducible representations (irreps) $\gamma$ are indicated in the figure. When the dynamical matrix $D_{\alpha_{1},\alpha_{2}}(\mathbf{q})$ is written in the symmetry adapted basis then the eigenvalue problem in Eq. 4 decomposes into $m_{\alpha}$ blocks according to the generalized Wigner-Eckart theorem Ludwig and Falter (1988); Wigner (1993):

[TABLE]

here we introduced the multi-index $\nu=\gamma,s_{\gamma},l$ (also the branch index in our case) which runs over all irreps $\gamma$, the row $l$ of the irrep and the the multiplicity $s_{\gamma}=1,\dots,m_{\gamma}$ of an irrep $\gamma$ in the basis $\alpha$. The unitary transformation from the Cartesian to the symmetry adapted basis is denoted as $S_{\nu\alpha}=\langle\nu|\alpha\rangle$. The symmetry-adapted basis is obtained employing the projecting technique of Ref. Ludwig and Falter, 1988. 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 $m_{\alpha}=1$. 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 ${\mathbf{q}}$ 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 $\alpha$ onto space group irreps for various ${\mathbf{q}}$-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 Landau and Lifshitz (1968); Landau (1937). 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 Stokes and Hatch (1984); Hatch et al. (1987); Stokes et al. (1991); Stokes and Hatch (1988); Stokes (1995). By listing the number of irreps of the space-group $Fm\overline{3}m$ at a given ${\mathbf{q}}$-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 Stokes et al. (1991); Stokes (1995). The supercells for the frozen-phonon calculations can be generated using the ISOTROPY Software Suite Stokes (2019). In the supercell geometries the ionic positions are shifted rigidly, where we used three displacement amplitudes of $0.02$, $0.06$, and $0.1$ atomic units.