Density functional perturbation theory for lattice dynamics with fully relativistic ultrasoft pseudopotentials: the magnetic case

TL;DR

This paper extends density functional perturbation theory to magnetic materials using fully relativistic ultrasoft pseudopotentials, incorporating time-reversal symmetry and magnetic point group operations for accurate lattice dynamics calculations.

Contribution

It introduces a novel formalism for including magnetic effects and symmetry operations in relativistic lattice dynamics calculations within density functional perturbation theory.

Findings

Validated approach against frozen phonon method in ferromagnetic Ni and Pt wire.

Demonstrated accurate treatment of magnetic symmetry in lattice dynamics.

Extended DFT perturbation theory to fully relativistic magnetic systems.

Abstract

We extend density functional perturbation theory for lattice dynamics with fully relativistic ultrasoft pseudopotentials to magnetic materials. Our approach is based on the application of the time-reversal operator to the Sternheimer linear system and to its self-consistent solutions. Moreover, we discuss how to include in the formalism the symmetry operations of the magnetic point group which require the time-reversal operator. We validate our implementation by comparison with the frozen phonon method in fcc Ni and in a monatomic ferromagnetic Pt wire.

Tables

Table 1. Table 1: Computed FR LDA phonon frequencies at 𝒒 = ( 0 , 2 ​ π / a , 0 ) 𝒒 0 2 𝜋 𝑎 0 \bm{q}=(0,2\pi/a,0) and 𝒒 = ( 0 , 0 , 2 ​ π / a ) 𝒒 0 0 2 𝜋 𝑎 \bm{q}=(0,0,2\pi/a) with DFPT and the frozen phonon method for fcc Ni. The magnetization is oriented along the z 𝑧 z axis. The subscripts indicate the polarization of the phonon modes.

	DFPT		Frozen phonon
	$\nu \left({\text{cm}}^{-1}\right)$		$\nu \left({\text{cm}}^{-1}\right)$
$T_x^{(0,1,0)}$	$232.438$		$232.691\pm 0.006$
$T_z^{(0,1,0)}$	$232.397$		$232.648\pm 0.006$
$T_{\left\{xy\right\}}^{(0,0,1)}$	$232.433$		$232.688\pm 0.006$

Table 2. Table 2: Computed FR LDA phonon frequencies at q = π / a 𝑞 𝜋 𝑎 q=\pi/a and q = π / 2 ​ a 𝑞 𝜋 2 𝑎 q=\pi/2a with DFPT and the frozen phonon method for a monatomic ferromagnetic Pt nanowire. The nanowire is oriented along the z 𝑧 z axis. Results are shown with both 𝒎 ∥ x conditional 𝒎 𝑥 \bm{m}\parallel x and 𝒎 ∥ z conditional 𝒎 𝑧 \bm{m}\parallel z . The subscripts indicate the polarization of the phonon modes.

q			$𝒎\parallel x$			$𝒎\parallel z$
		DFPT		Frozen phonon	DFPT		Frozen phonon
		$\nu \left({\text{cm}}^{-1}\right)$		$\nu \left({\text{cm}}^{-1}\right)$	$\nu \left({\text{cm}}^{-1}\right)$		$\nu \left({\text{cm}}^{-1}\right)$
	$T_x$	$36.51$		$37.02\pm 0.03$	$45.71$		$46.10\pm 0.03$
$\pi /a$	$T_y$	$37.00$		$37.34\pm 0.03$	$45.71$		$46.10\pm 0.03$
	$L$	$113.98$		$114.21\pm 0.03$	$110.30$		$110.51\pm 0.03$
	$T_x$	$25.1$		$25.5\pm 0.1$	$39.17$		$39.63\pm 0.03$
$\pi /2a$	$T_y$	$32.1$		$31.8\pm 0.1$	$39.17$		$39.65\pm 0.03$
	$L$	$54.2$		$53.8\pm 0.1$	$62.66$		$63.14\pm 0.03$

Equations

$$\sum _{{\sigma }^{\prime }}{H}^{\sigma {\sigma }^{\prime }}|{\Psi }_i^{{\sigma }^{\prime }}⟩={ϵ}_i\sum _{{\sigma }^{\prime }}{S}^{\sigma {\sigma }^{\prime }}|{\Psi }_i^{{\sigma }^{\prime }}⟩,$$

$$H^{\sigma {\sigma }^{\prime }}=-\frac{1}{2}{\nabla }^2{\delta }_{\sigma {\sigma }^{\prime }}+V_{KS}^{\sigma {\sigma }^{\prime }}.$$

$$V_{KS}^{\sigma {\sigma }^{\prime }}=V_{NL}^{\sigma {\sigma }^{\prime }}+\sum _{{\sigma }_1{\sigma }_2}\int d^{3}r{V}_{LOC}^{{\sigma }_1{\sigma }_2}(\mathbf{r})K_{\sigma {\sigma }^{\prime }}^{{\sigma }_1{\sigma }_2}(\mathbf{r}),$$

$$\begin{array}{rl}{K}_{\sigma {\sigma }^{\prime }}^{{\sigma }_1{\sigma }_2}(\mathbf{r},{\mathbf{r}}_1,{\mathbf{r}}_2)& =\delta (\mathbf{r}-{\mathbf{r}}_1)\delta (\mathbf{r}-{\mathbf{r}}_2){\delta }_{{\sigma }_1\sigma }{\delta }_{{\sigma }_2{\sigma }^{\prime }}\\ & +\sum _{Imn}\sum _{m_1{n}_1}{Q}_{mn}^I(\mathbf{r})f_{m_{1}m}^{\sigma {\sigma }_1}{\beta }_m^I({\mathbf{r}}_1)\\ & \times f_{n{n}_1}^{{\sigma }_2{\sigma }^{\prime }}{\beta }_{n_1}^{\ast I}({\mathbf{r}}_2),\end{array}$$

$$\mathcal{T}{H}^{[{\mathbf{B}}_{xc}]}{\mathcal{T}}^{†}=H^{[-{\mathbf{B}}_{xc}]}.$$

$$\begin{split}\frac{dn^{\sigma\sigma^{\prime}}(\bm{r})}{d\mu}=&\sum_{i}\sum_{\sigma_{1}\sigma_{2}}\bigg{[}\bra{\Psi_{i}^{\sigma_{1}}}K_{\sigma_{1}\sigma_{2}}^{\sigma\sigma^{\prime}}(\bm{r})\ket{\Delta^{\mu}\Psi_{i}^{\sigma_{2}}}\\ &+\sum_{\sigma_{3}\sigma_{4}}\bra{\left(\mathcal{T}\Psi_{i}\right)^{\sigma_{1}}}\mathcal{T}_{\sigma_{1}\sigma_{3}}\,K_{\sigma_{3}\sigma_{4}}^{\sigma^{\prime}\sigma}(\bm{r})\,\mathcal{T}_{\sigma_{4}\sigma_{2}}^{\dagger}\\ &\times\ket{\left(\mathcal{T}\Delta^{\mu}\Psi_{i}\right)^{\sigma_{2}}}\bigg{]}+\Delta^{\mu}n^{\sigma\sigma^{\prime}}(\bm{r}).\end{split}$$

$$\begin{split}&\frac{d^{2}E_{tot}^{(2)}}{d\mu d\lambda}=\sum_{i}\sum_{\sigma_{1}\sigma_{2}}\langle\Psi_{i}^{\sigma_{1}}\Biggm{\lvert}\frac{\partial V_{KS}^{[\bm{B}]\sigma_{1}\sigma_{2}}}{\partial\lambda}-\epsilon_{i}\frac{\partial S^{\sigma_{1}\sigma_{2}}}{\partial\lambda}\Biggm{\lvert}\Delta^{\mu}\Psi_{i}^{\sigma_{2}}\rangle\\ &+\sum_{i}\sum_{\sigma_{1}\sigma_{2}}\langle\left(\mathcal{T}\Psi_{i}\right)^{\sigma_{1}}\Biggm{\lvert}\frac{\partial V_{KS}^{[-\bm{B}]\sigma_{1}\sigma_{2}}}{\partial\lambda}-\epsilon_{i}\frac{\partial S^{\sigma_{1}\sigma_{2}}}{\partial\lambda}\Biggm{\lvert}\left(\mathcal{T}\Delta^{\mu}\Psi_{i}\right)^{\sigma_{2}}\rangle,\end{split}$$

$$\begin{split}&\sum_{\sigma_{2}}\left[H^{[-\mathbf{B}]\sigma_{1}\sigma_{2}}-\epsilon_{i}S^{\sigma_{1}\sigma_{2}}\right]\Bigm{\lvert}\left(\mathcal{T}\Delta^{\mu}\Psi_{i}\right)^{\sigma_{2}}\rangle=\\ &-\sum_{\sigma_{2}\sigma_{3}}\Pi_{c,i}^{\dagger\sigma_{1}\sigma_{2}}\left[\frac{dV_{KS}^{[-\mathbf{B}]\sigma_{2}\sigma_{3}}}{d\mu}-\epsilon_{i}\frac{\partial S^{\sigma_{2}\sigma_{3}}}{\partial\mu}\right]\Biggm{\lvert}\left(\mathcal{T}\Psi_{i}\right)^{\sigma_{3}}\rangle,\end{split}$$

$$\frac{d{n}^{\sigma {\sigma }^{\prime }}(\mathbf{r})}{d{u}_{s^{\prime }\beta }(\mathbf{q})}$$

$$|{\Delta }^{u_{s^{\prime }\beta }(\mathbf{q})}{\Psi }_{\mathbf{k}v}^{\sigma }⟩$$

$$\begin{split}\frac{dn^{\sigma\sigma^{\prime}}(\bm{r})}{du_{s^{\prime}\beta}(\bm{q})}&=\sum_{\bm{k}v}\sum_{\sigma_{1}\sigma_{2}}\bigg{[}\bra{\Psi_{\bm{k}v}^{\sigma_{1}}}K_{\sigma_{1}\sigma_{2}}^{\sigma\sigma^{\prime}}(\bm{r})\ket{\Delta^{u_{s^{\prime}\beta}(\bm{q})}\Psi_{\bm{k}v}^{\sigma_{2}}}\\ &+\sum_{\sigma_{3}\sigma_{4}}\bra{\left(\mathcal{T}\Psi_{-\bm{k}v}\right)^{\sigma_{1}}}T_{\sigma_{1}\sigma_{3}}K_{\sigma_{3}\sigma_{4}}^{\sigma^{\prime}\sigma}(\bm{r})T_{\sigma_{4}\sigma_{2}}^{\dagger}\\ &\times\ket{\left(\mathcal{T}\Delta^{u_{s^{\prime}\beta}(-\bm{q})}\Psi_{-\bm{k}v}\right)^{\sigma_{2}}}\bigg{]}+\Delta^{u_{s^{\prime}\beta}(\bm{q})}n^{\sigma\sigma^{\prime}}(\bm{r}).\end{split}$$

$$\begin{split}&\sum_{\sigma_{2}}\left[H^{[-\bm{B}]\sigma_{1}\sigma_{2}}-\epsilon_{-\bm{k}v}S^{\sigma_{1}\sigma_{2}}\right]\ket{(\mathcal{T}\Delta^{u_{s^{\prime}\beta}(-\bm{q})}\Psi_{-\bm{k}v})^{\sigma_{2}}}=\\ &-\sum_{\sigma_{2}}\Pi_{c,-\bm{k}v}^{\dagger\sigma_{1}\sigma_{2}}\Bigg{[}\ket{\phi_{T-\bm{k}v}^{u_{s^{\prime}\beta}(\bm{q})[-\bm{B}]\sigma_{2}}}\\ &+\sum_{\sigma_{3}}\sum_{\sigma_{4}\sigma_{5}}\int d^{3}r\frac{dV_{H,xc}^{[-\bm{B}]\sigma_{4}\sigma_{5}}(\bm{r})}{du_{s^{\prime}\beta}(\bm{q})}K^{\sigma_{4}\sigma_{5}}_{\sigma_{2}\sigma_{3}}(\bm{r})\ket{\left(\mathcal{T}\Psi_{-\bm{k}v}\right)^{\sigma_{3}}}\Bigg{]},\end{split}$$

$$|{\varphi }_{T-\mathbf{k}v}^{u_{s^{\prime }\beta }(\mathbf{q})[-\mathbf{B}]{\sigma }_2}⟩=\sum _{{\sigma }_3}\left(\frac{\partial V_{KS}^{[-\mathbf{B}]{\sigma }_2{\sigma }_3}}{\partial u_{s^{\prime }\beta }(\mathbf{q})}-{ϵ}_{-\mathbf{k}v}\frac{\partial S^{{\sigma }_2{\sigma }_3}}{\partial u_{s^{\prime }\beta }(\mathbf{q})}\right)|{\left(\mathcal{T}{\Psi }_{-\mathbf{k}v}\right)}^{{\sigma }_3}⟩,$$

$$\begin{split}\Phi^{(2)}_{\begin{subarray}{c}s\alpha\\ s^{\prime}\beta\end{subarray}}(\bm{q})&=\frac{1}{N}\sum_{\bm{k}v}\sum_{\sigma}\Bigg{[}\langle\phi^{u_{s\alpha}(\bm{q})[\bm{B}]\sigma}_{\bm{k}v}\lvert\Delta^{u_{s^{\prime}\beta}(\bm{q})}\Psi_{\bm{k}v}^{\sigma}\rangle\\ &+\langle\phi^{u_{s\alpha}(\bm{q})[-\bm{B}]\sigma}_{T-\bm{k}v}\lvert\left(\mathcal{T}\Delta^{u_{s^{\prime}\beta}(-\bm{q})}\Psi_{-\bm{k}v}\right)^{\sigma}\rangle\Bigg{]},\end{split}$$

$$\begin{split}\frac{{\widetilde{dn^{\sigma\sigma^{\prime}}}(\bm{r})}}{du_{s^{\prime}\beta}(\bm{q})}&=\sum_{\bm{k}v}\Bigg{[}u_{\bm{k}v}^{*\,\sigma}(\bm{r})\widetilde{\Delta^{u_{s^{\prime}\beta}(\bm{q})}}u_{\bm{k}v}^{\sigma^{\prime}}(\bm{r})+\sum_{\sigma_{1}\sigma_{2}}U_{\sigma^{\prime}\sigma_{1}}\\ &\times(\mathcal{T}u_{-\bm{k}v}(\bm{r}))^{*\,\sigma_{1}}(\mathcal{T}\widetilde{\Delta^{u_{s^{\prime}\beta}(-\bm{q})}}u_{-\bm{k}v}(\bm{r}))^{\sigma_{2}}U_{\sigma_{2}\sigma}^{\dagger}\Bigg{]}\\ &+\sum_{s^{\prime\prime}m\,n}\bigg{(}\widetilde{Q}_{mn}^{s^{\prime\prime}\bm{q}}(\bm{r})\Delta^{u_{s^{\prime}\beta}(\bm{q})}\rho_{mn}^{s^{\prime\prime}\sigma\sigma^{\prime}}\bigg{)}\\ &+\widetilde{\Delta^{u_{s^{\prime}\beta}(\bm{q})}}n^{\sigma\sigma^{\prime}}(\bm{r}),\end{split}$$

$${\tilde{Q}}_{mn}^{s^{\prime \prime }\mathbf{q}}(\mathbf{r})=e^{-ı\mathbf{q}\cdot \mathbf{r}}\sum _{\rho }{e}^{ı\mathbf{q}\cdot {\mathbf{R}}_{\rho }}{Q}_{mn}^I(\mathbf{r}),$$

$$\begin{split}\Delta^{u_{s^{\prime}\beta}(\bm{q})}\rho_{mn}^{s^{\prime\prime}\sigma\sigma^{\prime}}&=\sum_{m_{1}n_{1}}\sum_{\sigma_{1}\sigma_{2}}\sum_{\bm{k}v}\Big{(}\beta_{\bm{k}v}^{*\,s^{\prime\prime}m_{1}\sigma_{1}}f_{m_{1}m}^{\sigma_{1}\sigma}f_{nn_{1}}^{\sigma^{\prime}\sigma_{2}}\Delta^{u_{s^{\prime}\beta}(\bm{q})}\beta_{\bm{k}v}^{s^{\prime\prime}n_{1}\sigma_{2}}\\ &+\sum_{\sigma^{\prime\prime}\sigma^{\prime\prime\prime}}\beta_{T\,-\bm{k}v}^{*\,s^{\prime\prime}m_{1}\sigma_{1}}f_{m_{1}m}^{\sigma_{1}\sigma^{\prime\prime}}U_{\sigma^{\prime\prime}\sigma^{\prime}}U_{\sigma\sigma^{\prime\prime\prime}}^{\dagger}f_{nn_{1}}^{\sigma^{\prime\prime\prime}\sigma_{2}}\Delta^{u_{s^{\prime}\beta}(-\bm{q})}\beta_{T\,-\bm{k}v}^{s^{\prime\prime}n_{1}\sigma_{2}}\Big{)},\end{split}$$

$$\sum _{{\sigma }_3{\sigma }_4}{\mathcal{T}}_{{\sigma }_1{\sigma }_3}{K}_{{\sigma }_3{\sigma }_4}^{{\sigma }^{\prime }\sigma }(\mathbf{r}){\mathcal{T}}_{{\sigma }_4{\sigma }_2}^{†}=\sum _{{\sigma }^{\prime \prime }{\sigma }^{\prime \prime \prime }}{U}_{{\sigma }^{\prime }{\sigma }^{\prime \prime }}{K}_{{\sigma }_1{\sigma }_2}^{{\sigma }^{\prime \prime }{\sigma }^{\prime \prime \prime }}{U}_{{\sigma }^{\prime \prime \prime }\sigma }^{†},$$

$$\frac{\widetilde{dn}(\mathbf{r})}{d{u}_{s\alpha }(\mathbf{q})}$$

$$\frac{{\widetilde{dm}}_{\beta }(\mathbf{r})}{d{u}_{s\alpha }(\mathbf{q})}$$

$$\begin{split}&\sum_{\sigma_{2}}\left({H_{\bm{k}+\bm{q}}^{[-\bm{B}]\sigma_{1}\sigma_{2}}}-\epsilon_{-\bm{k}v}S_{\bm{k}+\bm{q}}^{\sigma_{1}\sigma_{2}}\right)\ket{(\mathcal{T}\widetilde{\Delta^{u_{s^{\prime}\beta}(-\bm{q})}}u_{-\bm{k}v})^{\sigma_{2}}}\\ &=-\sum_{\sigma_{2}}\Pi_{c,-\bm{k}v}^{\dagger\sigma_{1}\sigma_{2},-\bm{k}-\bm{q}}\Bigg{[}\ket{\widetilde{\phi}_{T-\bm{k}v}^{u_{s^{\prime}\beta}(\bm{q})[-\bm{B}]\sigma_{2}}}+\sum_{\sigma_{3}}\frac{\widetilde{dV}_{H,xc}^{[-\bm{B}]\sigma_{2}\sigma_{3}}}{du_{s^{\prime}\beta}(\bm{q})}\\ &\times\ket{(\mathcal{T}u_{\bm{k}v})^{\sigma_{3}}}+\sum_{\sigma_{3}}\sum_{s^{\prime\prime}m_{1}n_{1}}{}^{3}I_{s^{\prime\prime}m_{1}n_{1}}^{u_{s^{\prime}\beta}(\bm{q})[-\bm{B}]\sigma_{2}\sigma_{3}}\ket{\widetilde{\beta}_{m_{1}\bm{k}+\bm{q}}^{s^{\prime\prime}}}\\ &\times\beta_{T\,-\bm{k}v}^{s^{\prime\prime}n_{1}\sigma_{3}}\Bigg{]},\end{split}$$

$$⟨\mathbf{r}|{\tilde{\beta }}_{m_1\mathbf{k}+\mathbf{q}}^{s^{\prime \prime }}⟩=e^{-ı(\mathbf{k}+\mathbf{q})\cdot \mathbf{r}}\sum _{\rho }{e}^{ı(\mathbf{k}+\mathbf{q})\cdot {\mathbf{R}}_{\rho }}{\beta }_{m_1}^{s^{\prime \prime }}(\mathbf{r}-{\mathbf{R}}_I)$$

$${}^3{I}_{s^{\prime \prime }{m}_1{n}_1}^{u_{s^{\prime }\beta }(\mathbf{q})[-\mathbf{B}]{\sigma }_2{\sigma }_3}=\sum _{mn}\sum _{{\sigma }_4{\sigma }_5}{f}_{m_{1}m}^{{\sigma }_2{\sigma }_4}{f}_{n{n}_1}^{{\sigma }_5{\sigma }_3}\sum _{\alpha }{A}_{\alpha }^{{\sigma }_4{\sigma }_5}{}^3{I}_{s^{\prime \prime }mn\alpha }^{u_{s^{\prime }\beta }(\mathbf{q})[-\mathbf{B}]},$$

$${}^3{I}_{s^{\prime \prime }mn\alpha }^{u_{s^{\prime }\beta }(\mathbf{q})[-\mathbf{B}]}=\int Q_{mn}^{\gamma (s^{\prime \prime })}(\mathbf{r}-{\mathbf{d}}_{s^{\prime \prime }})\frac{d{C}_{\alpha }^{[-\mathbf{B}]}(\mathbf{r})}{d{u}_{s^{\prime }\beta }(\mathbf{q})}{d}^{3}r,$$

$$\mathcal{S}\mathbf{q}=\mathbf{q}+{\mathbf{G}}_{\mathcal{S}},$$

$$\mathcal{S}\mathbf{q}=-\mathbf{q}+{\mathbf{G}}_{\mathcal{S}},$$

$$\begin{split}\frac{\widetilde{dn(\bm{r})}}{du_{s^{\prime}\beta}(\bm{q})}=\frac{1}{N_{S}}\sum_{\{\mathcal{S}\lvert\bm{f}\}}\mathcal{O}_{\tau(\mathcal{S})}&\bigg{[}\sum_{\gamma}S_{\gamma\beta}\frac{\widetilde{dn}^{NS}(\{\mathcal{S}\lvert\bm{f}\}\bm{r})}{du_{\bar{s^{\prime}}\gamma}(\bm{q})}\\ &\times e^{\imath\bm{G}_{\mathcal{S}^{-1}}\cdot\bm{r}}e^{-\imath\bm{q}\cdot\bm{R}^{\mathcal{S}}_{s^{\prime}}}\bigg{]},\end{split}$$

$$\begin{split}\frac{\widetilde{dm_{\delta}(\bm{r})}}{du_{s^{\prime}\beta}(\bm{q})}&=\frac{1}{N_{S}}\sum_{\{\mathcal{S}\lvert\bm{f}\}}(-1)^{\tau(\mathcal{S})}\mathcal{O}_{\tau(\mathcal{S})}\Bigg{[}\sum_{\gamma\eta}\tilde{S}^{-1}_{\delta\eta}S_{\gamma\beta}\\ &\times\frac{\widetilde{dm}_{\eta}^{NS}(\{\mathcal{S}\lvert\bm{f}\}\bm{r})}{du_{\bar{s^{\prime}}\gamma}(\bm{q})}e^{\imath\bm{G}_{\mathcal{S}^{-1}}\cdot\bm{r}}e^{-\imath\bm{q}\cdot\bm{R}^{\mathcal{S}}_{s^{\prime}}}\Bigg{]},\end{split}$$

$$\begin{split}\Phi_{\begin{subarray}{c}s\alpha\\ s^{\prime}\beta\end{subarray}}(\bm{q})=\frac{1}{N_{S}}\sum_{\{\mathcal{S}\lvert\bm{f}\}}\mathcal{O}_{\tau(\mathcal{S})}&\Bigg{[}\sum_{\gamma\delta}S_{\gamma\alpha}\,S_{\delta\beta}\,\Phi_{\begin{subarray}{c}\bar{s}\gamma\\ \bar{s^{\prime}}\delta\end{subarray}}^{NS}(\bm{q})\\ &\times e^{\imath\bm{q}\cdot(\bm{R}^{\mathcal{S}}_{s}-\bm{R}^{\mathcal{S}}_{s^{\prime}})}\Bigg{]},\end{split}$$

Full text

I Introduction

Density Functional Perturbation Theory (DFPT) is widely used for the computation of the linear response properties of solids, and in particular for the study of their lattice dynamics. dfpt_review Some years ago, one of us applied DFPT us_fr_dfpt to a scheme based on plane waves and norm conserving (NC) or ultrasoft (US) us_vanderbilt pseudopotentials (PPs), that allow the introduction of spin-orbit effects within a fully relativistic (FR) density functional formalism and can be written in a form very similar to the scalar relativistic (SR) one. However, the theory presented in Ref. us_fr_dfpt, was implemented only for time-reversal invariant systems, and therefore applications that include spin-orbit so far have been limited to non-magnetic solids. Pb_phonon

In this work we extend this theory to the case of magnetic systems, by explicitly considering the presence of an exchange-correlation magnetic field in the Hamiltonian. DFPT equations in presence of a magnetic field have been recently written to calculate magnons with NC PPs in Refs. Giustino, and magnon, . In Ref. Giustino, , the charge density induced by a periodic perturbation was computed by using the response to a perturbation at wave vector $\bm{q}$ and the response to a perturbation at $-\bm{q}$, while in Ref. magnon, the problem at $-\bm{q}$ was not solved, but the time-reversal operator was used to obtain a second Sternheimer equation with a reversed magnetic field. The two formulations are equivalent. We generalize the theory of Ref. magnon, to a phonon perturbation, avoiding the study of the response at $-\bm{q}$, and write it in a form applicable to both NC and US PPs.

In presence of a magnetic field, the solid is invariant upon the symmetry operations of the magnetic space group. Some of these operations require the time-reversal operator. We discuss how to exploit these symmetries for the symmetrization of the induced charge and magnetization densities and for the dynamical matrix.

Finally, we validate our method in ferromagnetic fcc Ni first computing the phonon frequencies at the $X$ point in the Brillouin zone (BZ) and comparing with DFPT results, and then by computing the phonon dispersions. Moreover, we apply our method to a monatomic ferromagnetic Pt nanowire and compare its vibrational properties when the magnetization is parallel or perpendicular to the wire. Also for this case we compare the DFPT results to the frozen phonon method for a phonon wavevector $q=\pi/a$ and $q=\pi/2a$, and then we compute by DFPT the phonon dispersion in the 1D BZ.

II DFPT with fully relativistic US-PPs

In the Density Functional Theory (DFT) with FR US PPs, which accounts for spin-orbit effects, the minimization of the total energy functional leads to the Kohn-Sham KS equations for the two-component spinor wave functions us_fr_dfpt :

[TABLE]

where $S^{\sigma\sigma^{\prime}}$ is the overlap matrix needed in the US scheme, and the Hamiltonian $H^{\sigma\sigma^{\prime}}$ is:

[TABLE]

$V_{KS}^{\sigma\sigma^{\prime}}$ is the total Kohn-Sham potential:

[TABLE]

where:

[TABLE]

where $I=\{\rho,s^{\prime\prime}\}$, while $f_{m_{1}m}^{\sigma_{1}\sigma}$, $f_{nn_{1}}^{\sigma^{\prime}\sigma_{2}}$, as well as the indeces $m$, $n$, $m_{1}$, and $n_{1}$ are defined in Eq. (5) of Ref. us_fr_dfpt, . In particular, in Eq. (3) $V_{LOC}^{\sigma\sigma^{\prime}}=V_{eff}\delta_{\sigma\sigma^{\prime}}-\mu_{B}\bm{B}_{xc}\cdot\bm{\sigma}^{\sigma\sigma^{\prime}}$, and $V_{eff}=V_{loc}+V_{H}+V_{xc}$, is the sum of local, Hartree, and exchange and correlation potential, and $V_{NL}$ is the bare non local potential: in magnetic systems the spin, represented by the Pauli matrices, is coupled to the exchange-correlation magnetic field, $\bm{B}_{xc}$, defined as $\bm{B}_{xc}=-\delta E_{xc}/\delta\bm{m}$. This term breaks the time-reversal symmetry: indeed, introducing the time-reversal operator, $\mathcal{T}=\imath\sigma_{y}\mathcal{K}$, where $\mathcal{K}$ is the complex-conjugation operator and $\sigma_{y}$ is the Pauli matrix, the following relationship holds:

[TABLE]

We first exploit the time-reversal operator to rewrite the induced spin density. Following the notation of Ref. us_fr_dfpt, , we consider both the metallic and the insulating case. The change of the spin density induced by the variation of an external parameter $\mu$ (Eq. (10) of Ref. us_fr_dfpt, ) may be written as:

[TABLE]

$\Delta^{\mu}n^{\sigma\sigma^{\prime}}$ is defined as in the non-magnetic case and corresponds to the last two terms of Eq. (10) of Ref. us_fr_dfpt, .

The same idea can be applied to the second-order derivatives of the total energy. Only the term $d^{2}E_{tot}^{(2)}/d\mu d\lambda$ (Eq. (19) of Ref. us_fr_dfpt, ) needs to be rewritten by using $\mathcal{T}$:

[TABLE]

while the other contributions can be kept in their original form. Both Eq. (6) and Eq. (7) contain two unknown terms, namely $\ket{\Delta^{\mu}\Psi_{i}^{\sigma}}$ and $\ket{(\mathcal{T}\Delta^{\mu}\Psi_{i})^{\sigma}}$. The first can be computed by means of the Sternheimer linear system (Eqs. 13 and 14 of Ref. us_fr_dfpt, ), while the second is obtained by solving the following linear system:

[TABLE]

obtained by applying $\mathcal{T}$ to both sides of the Sternheimer linear system (Eq. (13) of Ref. us_fr_dfpt, ) and using the fact that $\mathcal{T}\,dV_{KS}^{[\bm{B}]}/d\mu\,\mathcal{T}^{\dagger}=dV_{KS}^{[-\bm{B}]}/d\mu$. In particular, here we introduced the time-reversed projector on the conduction manifold, namely $\Pi_{c,i}^{\dagger\sigma_{1}\sigma_{2}}=\sum_{\sigma\sigma^{\prime}}\mathcal{T}_{\sigma_{1}\sigma}P_{c,i}^{\dagger\sigma\sigma^{\prime}}\mathcal{T}^{\dagger}_{\sigma^{\prime}\sigma_{2}}$, similarly to what proposed in Ref. magnon, for the calculation of magnons. Eqs. (6), (7), and (8) are valid for the US PPs scheme, the NC formulation can be obtained as a particular case by writing $K_{\sigma_{1}\sigma_{2}}^{\sigma\sigma^{\prime}}(\bm{r},\bm{r}_{1},\bm{r}_{2})=\delta(\bm{r}-\bm{r}_{1})\,\delta(\bm{r}-\bm{r}_{2})\,\delta_{\sigma\sigma_{1}}\,\delta_{\sigma^{\prime}\sigma_{2}}$ and $S^{\sigma_{1}\sigma_{2}}=\delta_{\sigma_{1}\sigma_{2}}$. Moreover, the insulating case can be dealt with by putting $\widetilde{\theta}_{F,i}=1$ if the state is occupied or [math] if the state is empty (see Ref. us_dfpt, for the definition of $\widetilde{\theta}_{F,i}$).

III Phonons in periodic solids

In this section, we consider a phonon perturbation with a wavevector $\bm{q}$ perturbing a periodic solid, for which the wave functions $\Psi_{i}^{\sigma}(\bm{r})$ may be written in the Bloch form, $\Psi_{\bm{k}v}^{\sigma}(\bm{r})=e^{\imath\bm{k}\cdot\bm{r}}u_{\bm{k}v}^{\sigma}(\bm{r})$, where $u_{\bm{k}v}^{\sigma}(\bm{r})$ is lattice periodic. Following the discussion reported in Appendix A of Ref. us_dfpt, , we introduce the variation of the density and of the wave functions, induced by a phonon perturbation $u_{\nu s^{\prime}\beta}=1/\sqrt{M_{s^{\prime}}}\Re\left(u_{s^{\prime}\beta}(\bm{q})e^{\imath\bm{q}\cdot\bm{R_{\nu}}}\right)$, where $\Re$ indicates the real part, and define:

[TABLE]

Eq. (6) then becomes:

[TABLE]

In Eq. (11) the second term is identical to the first and is not explicitly computed in time-reversal invariant systems. The same holds for (14) below (for the dynamical matrix). Instead, for magnetic systems the two terms must be computed separately. In particular, the time-reversed response of the wave functions can be computed by solving the linear system Eq. (8), which becomes:

[TABLE]

where, similarly to Ref. us_dfpt, we defined:

[TABLE]

in which $\partial V_{KS}^{[-\bm{B}]\sigma_{2}\sigma_{3}}/\partial u_{s^{\prime}\beta}(\bm{q})$ and $\partial S^{\sigma_{2}\sigma_{3}}/\partial u_{s^{\prime}\beta}(\bm{q})$ are defined similarly to Eq. (10). The action of the time-reversal operator on the linear system changes the sign of the exchange and correlation magnetic field, which enters in the Hamiltonian, in $dV_{H,xc}^{[\bm{B}]\sigma_{4}\sigma_{5}}/du_{s^{\prime}\beta}(\bm{q})$, and in $\ket{\phi^{u_{s^{\prime}\beta}(\bm{q})[\bm{B}]\sigma_{2}}_{\bm{k}v}}$ through the third term in Eq. (9) of Ref. us_fr_dfpt, . We can then write the contribution to the dynamical matrix coming from $d^{2}E_{tot}^{(2)}/du_{\mu s\alpha}du_{\nu s^{\prime}\beta}$ in the following way:

[TABLE]

while the other contributions may be kept in their original form, discussed in Ref. us_dfpt, . Here, $N$ is the number of cells in the solid.

Eq. (11) may be further manipulated by writing explicitly $K_{\sigma_{1}\sigma_{2}}^{\sigma\sigma^{\prime}}(\bm{r})$ (Eq. (4)). Introducing the periodic parts of the Bloch functions and of the responses of the wave functions, we obtain the periodic part of the induced spin density (indicated with a tilde ($\sim$)):

[TABLE]

where we defined the quantities $\widetilde{Q}_{mn}^{s^{\prime\prime}\bm{q}}(\bm{r})$ and $\Delta^{u_{s^{\prime}\beta}(\bm{q})}\rho_{mn}^{s^{\prime\prime}\sigma\sigma^{\prime}}$ as:

[TABLE]

[TABLE]

where $\beta_{T\,-\bm{k}v}^{s^{\prime\prime}m_{1}\sigma_{1}}=e^{-\imath\bm{k}\cdot\bm{R}_{\rho}}\,\langle\beta^{I}_{m_{1}}\lvert(\mathcal{T}\Psi_{-\bm{k}v})^{\sigma_{1}}\rangle$, $\Delta^{u_{s^{\prime}\beta}(-\bm{q})}\beta_{T\,-\bm{k}v}^{s^{\prime\prime}n_{1}\sigma_{2}}=e^{-\imath(\bm{k}+\bm{q})\cdot\bm{R}_{\rho}}\langle\beta^{I}_{n_{1}}\lvert(\mathcal{T}\Delta^{u_{s^{\prime}\beta}(-\bm{q})}\Psi_{-\bm{k}v})^{\sigma_{2}}\rangle$, and we used the fact that:

[TABLE]

where $U=\imath\sigma_{y}$ is the unitary part of the time-reversal operator. The induced charge and magnetization densities can be computed from the induced spin density (Eq. (15)) as:

[TABLE]

In particular, for the induced charge density we use the fact that $\sum_{\sigma}U_{\sigma_{2}\sigma}^{\dagger}U_{\sigma\sigma_{1}}=\delta_{\sigma_{1}\sigma_{2}}$ and for the induced magnetization density the fact that $\sum_{\sigma\sigma^{\prime}}U_{\sigma_{2}\sigma}^{\dagger}\sigma_{\alpha}^{\sigma\sigma^{\prime}}U_{\sigma^{\prime}\sigma_{1}}=-\sigma_{\alpha}^{\sigma_{1}\sigma_{2}}$, so that the terms of the induced spin density that contain the time-reversed wave functions are subtracted in Eq. (20).

The linear system (12) may be written in terms of lattice periodic functions:

[TABLE]

where we defined:

[TABLE]

[TABLE]

and:

[TABLE]

where $\alpha=1,\dots,4$, $\bm{A}=(\mathbb{1},\sigma_{x},\sigma_{y},\sigma_{z})$, similarly to Ref. us_fr_pseudo, ($\mathbb{1}$ is the $2\times 2$ identity matrix), and $\bm{C}^{[\bm{B}]}=(V_{H,xc},-\mu_{B}B_{xc,x},-\mu_{B}B_{xc,y},-\mu_{B}B_{xc,z})$.

IV symmetrization

We indicate with $\{\mathcal{S}\lvert\bm{f}\}$ the symmetry operations of the space group of the crystal, where $\mathcal{S}$ is a rotation (proper or improper) and $\bm{f}$ is a translation. In a magnetic crystal, we have to consider also the operations $\mathcal{S}$ such that $\{\mathcal{T}\mathcal{S}\lvert\bm{f}\}$ is a symmetry of the crystal.

Since, for a phonon perturbation, the charge (and magnetization) density response and the dynamical matrix are computed at a given finite wave vector $\bm{q}$, we use as symmetry operations only those $N_{S}$ operations of the antiunitary small space group of $\bm{q}$, the subgroup of the antiunitary space group of the crystal, which contains the symmetry operations $\{\mathcal{S}\lvert\bm{f}\}$ such that:

[TABLE]

if $\{\mathcal{S}\lvert\bm{f}\}$ is a symmetry of the crystal, or:

[TABLE]

if $\{\mathcal{T}\mathcal{S}\lvert\bm{f}\}$ is a symmetry of the crystal. Here, $\bm{G}_{\mathcal{S}}$ is a reciprocal lattice vector that might appear when $\bm{q}$ is at zone border. In order to distinguish the two cases we introduce a variable $\tau(\mathcal{S})$ which may take the values $\tau=0$ or $\tau=1$ if Eq. (25) or Eq. (26) holds, respectively. We compute the unsymmetrized induced spin density by summing over the Irreducible Brillouin Zone (IBZ) in Eqs. (15) and (17), introducing a weight proportional to the number of elements in the star of $\bm{q}$. Then, we calculate the unsymmetrized induced charge and magnetization densities $\widetilde{dn}^{NS}(\bm{r})/du_{s^{\prime}\beta}(\bm{q})$ and $\widetilde{dm}_{\delta}^{NS}(\bm{r})/du_{s^{\prime}\beta}(\bm{q})$ using Eqs. (19) and (20). Finally, the complete responses are obtained through the following relationships:

[TABLE]

[TABLE]

where $\tilde{S}$ is the proper part of $S$, $\mathcal{O}_{\tau(\mathcal{S})}$ is the identity if $\tau(\mathcal{S})=0$, or $\mathcal{O}_{\tau(\mathcal{S})}=\mathcal{K}$ if $\tau(\mathcal{S})=1$. Moreover, $\bm{R}^{\mathcal{S}}_{s^{\prime}}=\mathcal{S}\bm{d}_{s^{\prime}}-\bm{d}_{\bar{s^{\prime}}}$, where $\bm{d}_{s^{\prime}}$ identifies the position of the atom $s^{\prime}$ with respect to the origin of its primitive cell, while $\bm{d}_{\bar{s^{\prime}}}$ is obtained by applying the rotation $\mathcal{S}$ to the atom $s^{\prime}$ ($\{\mathcal{S}\lvert\bm{f}\}(\bm{R}_{\nu}+\bm{d}_{s^{\prime}})=\bm{R}_{\bar{\nu}}+\bm{d}_{\bar{s^{\prime}}}$). Similarly, the dynamical matrix becomes:

[TABLE]

where $\Phi_{\begin{subarray}{c}\bar{s}\gamma\\ \bar{s^{\prime}}\delta\end{subarray}}^{NS}(\bm{q})$ is obtained summing over the IBZ in Eq. (14) and including the terms coming from $d^{2}E_{tot}^{(1)}/du_{\mu s\alpha}du_{\nu s^{\prime}\beta}$, $d^{2}E_{tot}^{(3)}/du_{\mu s\alpha}du_{\nu s^{\prime}\beta}$, and $d^{2}E_{tot}^{(4)}/du_{\mu s\alpha}du_{\nu s^{\prime}\beta}$, defined in Ref. us_fr_dfpt, .

V Applications

In this section we use the theory described above to compute the phonon dispersions of ferromagnetic fcc Ni and of a monatomic ferromagnetic Pt nanowire. We validate the theory by comparing the phonon frequencies obtained by diagonalizing the dynamical matrix (Eq. (29)) with those obtained by the frozen phonon method.

Computational details

First-principle calculations were performed within the Local Density Approximation (LDA) PZ and the Perdew-Burke-Ernzerhof (PBE) PBE schemes, as implemented in the Quantum ESPRESSO QE ; QE_2 and thermo_pw thermo_pw packages. The atoms are described by FR US PPs us_fr_pseudo , with 4$s$ and 3$d$ valence electrons for Ni (PPs Ni.rel-pz-n-rrkjus$\_$psl.0.1.UPF and Ni.rel-pbe-n-rrkjus$\_$psl.0.1.UPF from pslibrary 0.1) and with 6$s$ and 5$d$ valence electrons for Pt (PP Pt.rel-pz-n-rrkjus$\_$psl.1.0.0.UPF from pslibrary 1.0.0 pslibrary ; pslibrary_2 ).

DFPT calculations on ferromagnetic fcc Ni are at the theoretical LDA and PBE lattice constants, $a=6.483$ a.u. and $a=6.658$ a.u., which are $2.6$% and $0.02$% smaller than experiment COD ($a=6.659$ a.u.), respectively. The pseudo wavefunctions (charge density) are expanded in a plane waves basis set with a kinetic energy cut-off of $120$ ($600$) Ry. The BZ integrations were done using a shifted uniform Monkhorst-Pack k_grid $\bm{k}$-point mesh of $28\times 28\times 28$ points for the phonon calculations at a single wave vector $\bm{q}$. The same computational parameters, except the $\bm{k}$-point mesh which has been reduced to $18\times 18\times 18$ points, have been used for the phonon dispersions. The dynamical matrices have been computed by DFPT on a $6\times 6\times 6$ $\bm{q}$-point mesh, and Fourier interpolated to obtain the complete dispersions. Phonon frequencies of ferromagnetic Ni with the frozen phonon method, were calculated with a simple cubic supercell with $4$ Ni atoms. The kinetic energy cut-offs used are the same as for the DFPT calculations, while the BZ integrations were performed on a $\bm{k}$-point mesh of $24\times 24\times 24$ points. The presence of a Fermi surface has been dealt with by the Methfessel-Paxton smearing method MP with a smearing parameter $\sigma=0.02$ Ry.

DFPT calculations on monatomic ferromagnetic Pt nanowire were done at a stretched geometry with interatomic distance $d=4.927$ a.u.. The wire replicas have been separated by a vacuum space of $20$ a.u.. We have checked that by increasing the vacuum space the computed frequencies do not change more than $0.2$ cm*-1*. The system has been studied in a ferromagnetic configuration, with magnetization either parallel or perpendicular to the wire. The kinetic energy cut-off was $60$ ($400$) Ry for the wave functions (charge density). The $\bm{k}$-point mesh is a shifted uniform Monkhorst-Pack mesh of $300$ points. Frozen phonon calculations were performed with supercells with $2$ and $4$ Pt atoms, and Monkhorst-Pack meshes of $150$ and $75$ $\bm{k}$-points, respectively. The smearing parameter was $\sigma=0.002$ Ry.

Fcc Ni

We start our discussion from the computation of the phonon frequencies of ferromagnetic fcc Ni with the magnetization along $[001]$ (and with a magnitude that turns out to be $0.62\,\mu_{B}$ per atom), and compare the DFPT and the frozen phonon method at the $Y$ and $Z$ points. The results obtained are reported in Table 1. The frequencies of the transverse modes at $\bm{q}=(0,0,2\pi/a)$ (Z) are degenerate with both methods, as a consequence of the tetragonal magnetic symmetry ($D_{4h}(C_{4h})$): indeed both transverse modes have atomic displacements perpendicular to the magnetization. Instead, the transverse modes at $\bm{q}=(0,2\pi/a,0)$ (Y) show a small splitting of $0.04\,\text{cm}^{-1}$. The two modes are actually different because the atomic displacements are either parallel or perpendicular to the magnetization. A frequency splitting arises as a consequence of spin-orbit coupling. The DFPT and frozen phonon methods agree within $0.3\,\text{cm}^{-1}$. The DFPT and the frozen phonon method predict the same splitting, which however is small compared to the agreement of the absolute values of the frequencies obtained with the two methods, hence it is not possible to give an accurate quantitative prediction, but only an order of magnitude. With the kinetic energy cut-offs and $\bm{k}$-point mesh used, the frequencies obtained are converged within $5\times 10^{-3}\,\text{cm}^{-1}$, the same order of magnitude as the errorbar reported in Table 1 and due to the fit.

In Fig. 1 we show the complete phonon dispersion of fcc Ni obtained by DFPT. Both LDA and PBE theoretical dispersions are shown, together with inelastic neutron scattering data exp_phonon_Ni . The agreement between the LDA result and the experiment is poor, mainly because LDA underestimates the lattice constant: the highest frequencies of the dispersion (e.g. at the $X$ and $L$ points) are about $30\,\text{cm}^{-1}$ higher than the experiment. On the other hand, the PBE phonon dispersions are in excellent agreement with the experiment. Note however that this agreement is slightly worsened by temperature effects Ni_temp_phonon not included in the present study.

Pt monatomic wire

In this section we consider a monatomic Pt nanowire, a metal with ferromagnetic ordering. It has been shown Pt_nature ; Pt_PRB that at its equilibrium geometry (atomic distance $d=4.441$ a.u.) the system shows a colossal magnetic anisotropy, since the preferred orientation of the magnetization is parallel to the wire and the magnetization vanishes when forced to be perpendicular to the wire. Instead, for stretched geometries with atomic distance higher than $4.913$ a.u. a non-zero magnetization perpendicular to the wire is allowed. Here we consider a stretched geometry with $d=4.927$ a.u. and compute the phonon dispersions with both a magnetization parallel and perpendicular to the wire. In the following the nanowire is along the $z$ direction. In Table 2 we compare the phonon frequencies, at $q=\pi/a$ and $q=\pi/2a$ with $\bm{m}\parallel x$ and $\bm{m}\parallel z$, computed by the DFPT and with the frozen phonon method. With a magnetization $\bm{m}\parallel z$ ($m=0.65\,\mu_{B}$ per atom), the frequencies of the transverse modes are degenerate, while with $\bm{m}\parallel x$ ($m=0.13\,\mu_{B}$ per atom) at $q=\pi/a$ the two transverse modes show a splitting of about $0.5\,\text{cm}^{-1}$, which is of the same order of magnitude as the overall agreement of the two methods. At $q=\pi/2a$ this splitting is about $7\,\text{cm}^{-1}$, one order of magnitude larger than at $q=\pi/a$. In both cases the polarization of the transverse mode with higher frequency is perpendicular to the magnetization. As discussed above for fcc Ni, the two transverse modes are not equivalent due to the presence of the magnetization and of spin-orbit coupling. Pt atoms are heavier than Ni and show a stronger spin-orbit interaction: indeed, the splitting reported for Pt is $1-2$ orders of magnitudes higher than in Ni. The DFPT and frozen phonon results agree within $0.4\,\text{cm}^{-1}$ on average. As before, the errorbars reported in Table 2 come from the linear fit. With the kinetic energy cut-offs and the $\bm{k}$-point mesh used all the frequencies reported are converged within $0.03\,\text{cm}^{-1}$

In Fig. 2 we show the phonon branches along $\Gamma-Z$ for a ferromagnetic wire with magnetization parallel (left panel) or perpendicular to the wire (right panel). The two dispersions show evident differences: at $q=\pi/a$, the longitudinal mode for the wire with $\bm{m}\parallel z$ is lower in frequency than for the wire with $\bm{m}\parallel x$, while the transverse modes are higher in frequency. In the central part of the BZ, around $q=\pi/2a$, the longitudinal mode of the wire with $\bm{m}\parallel z$ has a higher frequency at the Kohn anomaly than the wire with $\bm{m}\parallel x$, while the transverse modes show a Kohn anomaly only for $\bm{m}\parallel x$. We remark that at the stretched geometry studied ($d=4.927$ a.u.) the phonon modes are still stable, but the range of atomic distances at which both modes are stable is quite narrow.

VI Conclusions

We extended the DFPT for lattice dynamics with FR NC and US PPs to the magnetic case. We validated the theory by comparing the DFPT to the frozen phonon method for ferromagnetic fcc Ni and for a monatomic ferromagnetic Pt nanowire. The agreement between the two methods is within $0.5\,\text{cm}^{-1}$. For both systems, we computed by DFPT also the complete phonon dispersions and discussed their features, showing that magnetism together with spin-orbit coupling may lift the degeneracy of some phonon modes. For our systems these splittings range from $10^{-2}\,\text{cm}^{-1}$ (in Ni) to a few $\text{cm}^{-1}$ (in Pt nanowire).

Acknowledgments

Computational facilities have been provided by SISSA through its Linux Cluster and ITCS and by CINECA through the SISSA-CINECA 2018-2019 Agreement.