Density functional perturbation theory within non-collinear magnetism

TL;DR

This paper extends density functional perturbation theory to non-collinear magnetism by developing methods to compute exchange-correlation derivatives in a non-collinear framework, enabling accurate simulations of magnetic materials.

Contribution

The paper introduces two approaches to transform non-collinear exchange-correlation derivatives into a local collinear basis for perturbation theory calculations.

Findings

Methods successfully applied to atomic displacement perturbations.

Methods enable second-order electric field perturbation calculations.

Comparison shows trade-offs between Taylor expansion and explicit evaluation.

Abstract

We extend the density functional perturbation theory formalism to the case of non-collinear magnetism. The main problem comes with the exchange-correlation (XC) potential derivatives, which are the only ones that are affected by the non-collinearity of the system. Most of the present XC functionals are constructed at the collinear level, such that the off-diagonal (containing magnetization densities along $x$ and $y$ directions) derivatives cannot be calculated simply in the non-collinear framework. To solve this problem, we consider here possibilities to transform the non-collinear XC derivatives to a local collinear basis, where the $z$ axis is aligned with the local magnetization at each point. The two methods we explore are i) expanding the spin rotation matrix as a Taylor series, ii) evaluating explicitly the XC for the local density approximation through an analytical expression of the expansion terms. We compare the two methods and describe their practical implementation. We show their application for atomic displacement and electric field perturbations at the second order, within the norm-conserving pseudopotential methods.

Tables

Table 1. Table 1: Atomic displacement tests comparing FD to the three different transformation methods for DFPT: (i) , (i’) , (ii) . Values give the second derivatives of the energy with respect to reduced displacements (values in Ha). SCI is the number of the SCF cycles needed to converge with a tolerance Δ ​ V Δ 𝑉 \Delta V (Ha). FD are performed using 5 discretization points for the Frozen part and 3 points for the Total one. Spin-orbit coupling is not included.

Frozen
RM	$\widehat{𝐦}$	E(2) (FD)	E(2) (DFPT)	$\Delta$V	# SCI
(i)	$\hat{x}$	26330.12639	26330.12833	9 10-9	32
(i’)			26330.12833	5 10-9	40
(ii)			26330.12833	5 10-5	80
(i)	$\hat{y}$	26330.12639	26330.12833	1 10-8	39
(i’)			26330.12833	9 10-9	39
(ii)			26330.12833	1 10-4	80
(i)	$\hat{z}$	26330.12640	26330.12833	7 10-9	30
(i’)			26330.12833	6 10-9	30
(ii)			26330.12833	8 10-7	80
Total
RM	$\widehat{𝐦}$	E(2) (FD)	E(2) (DFPT)	$\Delta$V	# SCI
(i)	$\hat{x}$	16.13695	16.13785	7 10-9	36
(i’)			16.14042	3 10-9	39
(ii)			16.13785	3 10-5	80
(i)	$\hat{y}$	16.13693	16.13785	6 10-9	39
(i’)			16.14042	4 10-9	39
(ii)			16.13785	7 10-5	80
(i)	$\hat{z}$	16.13694	16.13785	4 10-9	32
(i’)			16.14042	1 10-8	40
(ii)			16.13785	3 10-5	80

Table 2. Table 2: Calculated phonon frequencies (cm -1 ) of Cr 2 O 3 for collinear case and non-collinear with magnetic moments along the z 𝑧 z directions (no SOC). We compare the results given by method (i), (i’) and (ii), which all give the same result and are, as expected for collinear antiferromagnet, extremely close to the collinear case.

mode #	Collinear	non-coll. (i)	non-coll. (i’)	non-coll. (ii)
1-2	234	234	234	234
3	236	236	236	236
4	246	246	246	246
5-6	279	280	280	280
7-8	300	300	300	300
9-10	357	357	357	357
11	361	361	361	361
12	382	382	382	382
13-14	414	414	414	414
15	437	437	437	437
16	447	447	447	447
17-18	464	465	465	465
19	466	466	466	466
20-21	473	473	473	473
22-23	522	522	522	522
24	537	537	537	537
25-26	545	545	545	545
27	567	567	567	567

Table 3. Table 3: Calculated phonon frequencies (cm -1 ) of α 𝛼 \alpha -RuCl 3 in its FM phase for collinear case and non-collinear with magnetic moments along the x 𝑥 x and z 𝑧 z directions.

mode #	Label	Collinear	non-collinear ($\hat{x}$)	non-collinear ($\hat{z}$)
1	${\Gamma }_2^{+}$	20$i$	19$i$	21$i$
2-3	${\Gamma }_5^{+}$	123	123	123
4	${\Gamma }_4^{-}$	146	146	147
5-6	${\Gamma }_5^{-}$	148	148	147
7	${\Gamma }_3^{-}$	169	166	158
8-9	${\Gamma }_6^{-}$	176	174	175
10-11	${\Gamma }_6^{+}$	249	249	248
12	${\Gamma }_4^{+}$	279	278	277
13-14	${\Gamma }_5^{-}$	283	283	282
15-16	${\Gamma }_5^{+}$	291	291	291
17	${\Gamma }_2^{-}$	298	298	296
18-19	${\Gamma }_6^{-}$	333	334	335
20	${\Gamma }_3^{-}$	355	347	351
21	${\Gamma }_1^{+}$	355	354	352

Equations

$$\mathrm{E}[\rho ]=\sum _n⟨{\psi }_n⟩T+V_{ext}{\psi }_n+{\mathrm{E}}_{\mathrm{H}}[\rho ]+{\mathrm{E}}_{xc}[\rho ].$$

$$\hat{\rho }(\mathbf{r})=\sum _n{\psi }_n^{\ast }(\mathbf{r}){\psi }_n(\mathbf{r}).$$

$$\mathrm{H}=T+V_{ext}+V_H+V_{xc}.$$

$$V_{xc}=\frac{\delta {\mathrm{E}}_{xc}[\rho ]}{\delta \rho (\mathbf{r})}.$$

$$X(\lambda )=X^{(0)}+\lambda X^{(1)}+{\lambda }^2{X}^{(2)}+{\lambda }^3{X}^{(3)}+...$$

$$X^{(N)}=\dfrac{1}{N!}\dfrac{d^{N}X}{d\lambda^{N}}\Bigr{|}_{\lambda=0}\quad.$$

$$V_{xc}^{(i)}=\left(\dfrac{\delta E_{xc}[\rho(\lambda)]}{\delta\rho({\bf r})}\right)^{(i)}=\dfrac{1}{i!}\left(\dfrac{d^{(i)}}{d\lambda^{(i)}}\dfrac{\delta E_{xc}[\rho(\lambda)]}{{\delta\rho({\bf r})}}\right)\Bigr{|}_{\lambda=0}\ .$$

$$\left[\mathrm{H}(\lambda )-E_n(\lambda )\right]|{\psi }_n(\lambda )⟩=0,$$

$$P_c(H-E_n{)}^{(0)}{P}_c|{\psi }_n^{(1)}⟩+P_c(H-E_n{)}^{(1)}|{\psi }_n^{(0)}⟩=0;$$

$${\mathrm{H}}^{(1)}$$

$$⟨{\psi }_n^{(0)}⟩{\psi }_n^{(1)}+⟨{\psi }_n^{(1)}⟩{\psi }_n^{(0)}=0.$$

$$\left(\begin{array}[]{cc}H^{\alpha\alpha}&V_{xc}^{\alpha\beta}\\[15.05624pt] V_{xc}^{\beta\alpha}&H^{\beta\beta}\end{array}\right)\ \left(\begin{array}[]{c}\ket{\psi_{n}^{\alpha}}\\[16.06pt] \ket{\psi_{n}^{\beta}}\end{array}\right)\ =E_{n}\left(\begin{array}[]{c}\ket{\psi_{n}^{\alpha}}\\[16.06pt] \ket{\psi_{n}^{\beta}}\end{array}\right),$$

$$V_{xc}^{\alpha \beta }\left[\hat{\rho }(\mathbf{r})\right]=\frac{\delta E_{xc}\left[\hat{\rho }(\mathbf{r})\right]}{\delta {\rho }^{\alpha \beta }(\mathbf{r})}.$$

$$V_{xc}^{\alpha \beta }$$

$$+$$

$${\rho }_{11}$$

$${\rho }_{22}$$

$${\rho }_{12}$$

$$\hat{\rho }=|{\psi }^{(0)}⟩⟨{\psi }^{(0)}|.$$

$${\rho }_{\alpha \beta }=⟨\alpha ⟩{\psi }^{(0)}⟨{\psi }^{(0)}⟩\beta ,$$

$$\rho$$

$$=$$

$${\rho }_{\alpha \beta }$$

$$=$$

$$n$$

$$m_x$$

$$m_y$$

$$m_z$$

$$n^{(0)}(\mathbf{r})=\frac{1}{(2\pi {)}^3}{\int }_{\text{BZ}}\sum _i^{occ}\sum _{s=1}^2{\psi }_{i\mathbf{k}s}^{(0)\ast }(\mathbf{r}){\psi }_{i\mathbf{k}s}^{(0)}(\mathbf{r})d\mathbf{k},$$

$$N_1$$

$$N_2$$

$$D_1$$

$$D_2$$

$$\sum _{\alpha \beta }{U}_{i\alpha }^{†}{\hat{\rho }}^{\alpha \beta }{U}_{\beta j}={\rho }_{ij}{\delta }_{ij},$$

$$\hat{U}=\left(\begin{array}[]{cc}\dfrac{m_{z}+m}{D_{1}}&\dfrac{m_{z}-m}{D_{2}}\\[15.05624pt] \dfrac{m_{x}+i\ m_{y}}{D_{1}}&\dfrac{m_{x}+i\ m_{y}}{D_{2}}\end{array}\right).$$

$$\hat{U}^{-1}=\left(\begin{array}[]{cc}\dfrac{D_{1}}{2m}&\dfrac{\left(m-m_{z}\right)D_{1}}{2m\left(m_{x}+i\ m_{y}\right)}\\[15.05624pt] -\dfrac{D_{2}}{2m}&\dfrac{\left(m+m_{z}\right)D_{2}}{2m\left(m_{x}+i\ m_{y}\right)}\end{array}\right).$$

Full text

I Introduction

The density functional perturbation theory (DFPT) is the application of linear response theory to the density functional theory (DFT) formalism. Its development traces back to the birth of many DFT codes during the 80’s Baroni et al. (2001); Gonze (1995); Baroni et al. (1987). The first applications were for the calculations of phonon spectra, dielectric constants, Born effective charges and related properties Gonze and Lee (1997). These responses require to expand the DFT as a function of atomic displacements and electric field perturbations. Later on, strain perturbations were added to the list, which allows the calculation of elastic and piezoelectric tensors Hamann et al. (2005); Wu et al. (2005). DFPT represented a breakthrough in the calculation of many physical properties since, once the equations are implemented, the resulting work-flow is very simple for the end user; this also simplifies high-throughput calculations Petousis et al. (2016); Hautier et al. (2012).

However, in spite of the past decade’s increase of interest for non-collinear magnetic materials, DFPT has not yet been generalized for non-collinear magnets. This is due to the fact that non-collinear magnetism within DFT is much more computationally demanding (the density has four components instead of two for collinear magnetism) and it requires to include the spin-orbit coupling (SOC), which, in present implementations of DFT codes, often reduces the number of usable symmetries Peralta et al. (2007); Hobbs et al. (2000); Kubler et al. (1988). Additionally, non-collinear magnetism with SOC works in very low energy ranges (from 1 to hundreds of $\mu$eV), which require highly accurate simulations Bousquet and Spaldin (2010). Non-collinear ground state (GS) calculations (spin and lattice) have been first developed in the 80’s Kubler et al. (1988) but became widespread only recently thanks to increased supercomputer power and the improvement of code scalability Naveh and Kronik (2009); Peralta et al. (2007); Hobbs et al. (2000); Verstraete et al. (2008); Dal Corso (2010). It is thus timely to extend DFPT to non-collinear magnetism, for phonon spectra, elasticity or dielectric responses.

In this paper, we derive the formalism of DFPT in the presence of non-collinear magnetism. The main difficulty is represented by the XC term, for which the functionals have been developed using a collinear framework. Only a few attempts have been made to go beyond this and parameterize directly a non-collinear functional Eich and Gross (2013); Eich et al. (2013); Bulik et al. (2013); Scalmani and Frisch (2012). We show here two different approaches to perform the energy derivatives in the non-collinear regime starting with the collinear $E_{xc}$ functionals. We start with a rapid reminder of collinear and of non-collinear DFT formalisms, followed by the description of two formally equivalent methods to derive the XC potential from collinear $E_{xc}$ functionals. We also show the specific derivations (using norm-conserving pseudo potentials), for atomic displacement and electric field perturbations. In each case, we show how the two methods compare in precision and efficiency.

II Basics of DFPT

The framework of DFTKohn and Sham (1965); Hohenberg and Kohn (1964) allows to solve a many body problem as a single particle in an effective self-consistent potential. In this formalism Kohn (1999), the reference Ground State (GS) energy is estimated by minimizing with respect to the Kohn-Sham (KS) wavefunctions $\psi_{\alpha}$ the following expression:

[TABLE]

Here, $n$ runs over all the occupied states, $T$ is the kinetic energy operator, $v_{ext}$ is the external potential (the sum of ionic potentials, which can be non-local, and eventual applied fields), $E_{\rm H}$ and $E_{xc}$ are the interaction energy terms, namely the Hartree ($H$) and XC functionals. These two last terms contain the approximations to the many body electronic interactions, as functionals of the density operator:

[TABLE]

The Hamiltonian related to Eq. (1), needed to self-consistently solve the KS equation, is:

[TABLE]

The solution of this equation depends on the way the last term is approximated in the so-called XC potential $V_{xc}$. Interestingly, $V_{xc}$ can be written as the functional derivative of the XC energy with respect to the density:

[TABLE]

When the system is subject to a perturbation, its effects can be related to a (small) parameter $\lambda$ on which any generic observable $X(\lambda)$ depends. Expanding explicitly the perturbation series we have:

[TABLE]

where the expansion coefficients are related to the $N^{\rm th}$ order derivative of $X(\lambda)$ with respect to $\lambda$ by a scalar coefficient:

[TABLE]

Our objective is to evaluate the modification of the GS energy E or density $\rho({\bf r})$ caused by the presence of the external perturbation, through the computation of this expansion.

At the $i^{\rm th}$ generic order in $\lambda$, the XC potential can be written as:

[TABLE]

Truncating the expansion at the first order, the equation that has to be self-consistently solved is obtained expanding the single particle KS equation Gonze (1995); Gonze and Lee (1997):

[TABLE]

giving the so-called Sternheimer equation Sternheimer (1954):

[TABLE]

here, $P_{c}$ is the projector upon the empty states (conduction bands), $H^{(0)}$, $E_{\alpha}^{(0)}$ and $\psi^{(0)}$ are, respectively, the GS Hamiltonian, eigenvalues and wavefunctions and $H^{(1)}$ is the first-order Hamiltonian:

[TABLE]

The orthonormalization condition that holds for the GS wavefunctions $\braket{\psi_{n}}{\psi_{n}}=1$, becomes to first order:

[TABLE]

Once the framework has been specified (DFT for the search of the GS and DFPT to explore the effects of the perturbations), it is very useful to write the operators in the spinorial form: each operator belonging to the KS equation is represented as a 2$\times$2 matrix. In this way, we can rewrite the Eq. (8) as:

[TABLE]

here, the indices $\alpha$ and $\beta$ run over the up and down spin components (or major and minor Pauli spinors) and $n$ over the bands. As we can see, the off-diagonal part of the hamiltonian comes only from the XC potential: all other terms of H (kinetic energy, Hartree and pseudo-potential components) do not cross-couple up and down spin channels. At the zeroth order, the general $\alpha\beta$ element of the XC potential is obtained simply re-writing Eq. (4):

[TABLE]

and since the (spin) density matrix is diagonal ($\rho_{\alpha\beta}({\bf r})=\hat{\rho}^{\alpha\beta}\delta_{\alpha\beta}$) only if the magnetization is directed along the quantization axis $z$ (then the magnetization ${\bf m}({\bf r})$ = m${}_{z}({\bf r})$), it is useful to diagonalize the 2$\times$2 matrix in the previous section, finding the transformation $\hat{U}$ (often referred to as the “spin-$1/2$ rotation matrix”). Using the new basis, the XC potential becomes:

[TABLE]

where $\rho^{\uparrow(\downarrow)}({\bf r})=1/2\left(n\pm m\right)$, where $n$ is the charge density; $\delta_{\alpha\beta}$ is the Kronecker delta, $\hat{{\bf m}}={\bf m}/m$ is the magnetization versor and $m$ the norm of ${\bf m}$; $\boldsymbol{\sigma}_{\alpha\beta}(\sigma^{x}_{\alpha\beta},\sigma^{y}_{\alpha\beta},\sigma^{z}_{\alpha\beta})$ are the Pauli matrices.

In the following, we extend the DFPT formalism to non-collinear magnetism using the spinorial wavefunctions. With respect to collinear magnetism, i. e., where the direction of all the spins is along a well defined quantization axis, there are at least two issues to be addressed: 1) the calculation of the derivative of the XC potential $V^{(1)}_{xc\ \alpha\beta}$ and 2) the computation of the density derivatives $\rho^{(1)}_{\alpha\beta}$ from the derivative of the KS wavefunctions $\psi^{(1)}_{n\ \alpha}$.

III Non-collinear DFT framework

Our goal is now to generalize the DFPT equations with spinorial wavefunctions, to take into account the NCM and thus a generic magnetization direction. The main concern of such a generalization is to treat the derivatives of the XC term (see above). A first problem arises from the fact that in spin-dependent DFT the XC functionals are constructed on a collinear basis where the magnetization is considered to be have the same orientation at each point in space. The standard way to handle non-collinear DFT with collinear XC functionals is to locally rotate the magnetization density with respect the quantization axis at each point in space, so that the derivatives of the XC term can be treated as locally collinear Hobbs et al. (2000), and then rotated back to the “laboratory frame”. To avoid this “local collinearity approximation”, some attempts have been made to build specific non-collinear XC functionals Eich and Gross (2013); Eich et al. (2013); Bulik et al. (2013); Scalmani and Frisch (2012). In this paper, we will restrict ourselves to the case of “local collinearity approximation” and use the common collinear XC functionals.

In a collinear magnetic system, the magnetization at each point r in space is considered to be directed along a fixed quantization axis (usually $z$). In this picture, the density matrix has two components only, which are often specified using $\alpha,\beta=1,2$, to indicate the spin up $\ket{1}$ and down $\ket{2}$ states:

[TABLE]

In these expressions and the following ones, we omit the (0) or (1) to indicate the order of the perturbation, since the form of the density matrix is generally valid.

In the case of NCM, the density matrix operator can be written in terms of the GS spinorial wavefunctions $\psi^{(0)}$ as follows:

[TABLE]

In the representation of the spin, projecting the spinors on the orthonormal spin basis $\ket{\alpha}$ and $\ket{\beta}$, $\hat{\rho}_{\alpha\beta}$ can be evaluated explicitly in terms of spin components:

[TABLE]

or in matrix form:

[TABLE]

where $\braket{\alpha}{\psi^{(0)}}=\psi_{\alpha}^{(0)}$ and $\braket{\psi^{(0)}}{\beta}=\psi_{\beta}^{(0)*}$, $m_{i}$ is the $i=x,y,z$ component of the spin vector, or magnetization density, and $n$ is the electronic density. This expression for the density matrix is defined by its relation with the Pauli matrices $\boldsymbol{\sigma}_{\alpha\beta}$ Wang and Ziegler (2006), which gives the following expression:

[TABLE]

where the electronic and magnetization components can be written using Eq. (22):

[TABLE]

The expression of the GS charge density $n$ in Eq. (28) can be written as specified in Ref. Gonze (1995):

[TABLE]

where $i$ runs over all the occupied states, $s$ is the spin component (if the system is non-magnetic $s=1$) and the integral is performed on the whole Brillouin Zone (BZ).

IV DFPT within non-collinear magnetism

In order to address a NCM problem, we want to map the overall 2$\times$2 density matrix form (Eq. (25)) onto a local collinear case. In the following we develop two possibilities to do this mapping. Both techniques rely on the idea of first evaluating the perturbed XC potential ($V^{(1)}_{xc}({\bf r})$) in a local reference frame, whose $z$-axis is aligned with the direction of local perturbed magnetization at each point r of the real space. Such a choice of local coordinates reduces the problem to that of collinear magnetism, and after the evaluation of the corresponding “collinear” $V^{(1)}_{xc}({\bf r})$ values, we can restore its non-collinear form applying the backwards rotation that aligns the local reference frame with the globally defined Cartesian axes. The first order potential can be then evaluated based on either (i) the $\lambda$ expansion of the rotation matrices and corresponding series of the collinear $V^{(1)}_{xc}$ expression or (ii) by direct evaluation of the $\lambda$ series of the final non-collinear $V^{(1)}_{xc}$ expression. Essentially, we develop either the expansion of the general expression (case (i)) or of the final result (case (ii)). Although both methods give the same results, the first method (i) is more general: the procedure to evaluate $V^{(1)}_{xc}$ can be applied without any modification of the XC potential (even a non-diagonal one), while the second method requires the explicit re-evaluation of the $\lambda$ series for each particular choice of XC flavour (LDA or GGA). For the sake of completeness we provide below the development of both schemes (see sections IV.1 and IV.1.2 for methods (i) (and an alternative one that we call (i’)) and (ii), respectively). Furthermore, implementation of both methods allows us to validate the correctness of the implemented code by contrasting results obtained using the two different techniques.

IV.1 (i) Expansion of the rotation matrix method

IV.1.1 Magnetization density local transformation

To project the non-collinear density onto the local axis of magnetization, we need to determine the direction of the local spin-quantization axis to build a diagonal collinear density. In the case of a general GS density matrix $\hat{\rho}$ of the form Eq. (25), the eigenvalues are given by $\rho_{\uparrow,\downarrow}=1/2\left(n\pm|m|\right)$ and the two corresponding normalized and complex eigenvectors are:

[TABLE]

where $m=|{\bf m}|=\sqrt{m_{x}^{2}+m_{y}^{2}+m_{z}^{2}}$. The two $N_{1}$ and $N_{2}$ eigenvectors define the spin-1/2 unitary transformation matrix $\hat{U}$ diagonalizing $\hat{\rho}$. Then, the transformation can be detailed as:

[TABLE]

here $\alpha,\beta,i,j$ are indices running on the two possible spin states $\ket{1}$ or $\ket{2}$ (as previously defined). The expression of $\hat{U}$ is thus:

[TABLE]

The inverse of this matrix, such that $\hat{U}^{-1}\hat{U}={\bf I}$, with ${\bf I}$ the identity matrix, is:

[TABLE]

Moreover, at all orders in $\lambda$, we have the unitarity property: $\hat{U}^{\dagger}=\hat{U}^{-1}$. Importantly, in the last two expressions the presence of the ${\bf m}$ components in the denominator does not cause a divergence if ${\bf m}=0$, since the dependence on the magnetization of the numerator is of the same order. To the first order in the perturbation, the unitarity condition continues to hold:

[TABLE]

and can be expanded:

[TABLE]

Developing this expression, neglecting the second order perturbations ($\lambda^{2}$) and higher terms, we obtain the following relationship:

[TABLE]

All these properties will be exploited in order to obtain the expression of the perturbed XC potential and density matrix in the following part of the manuscript.

IV.1.2 XC potential and density matrix

As we have seen for a general observable $X$ in Sec. II, the perturbed XC potential (in the following we will omit the “$xc$” subscript in $V_{xc}$) can be expanded as a Taylor series in $\lambda$ and truncated at the first order:

[TABLE]

At each point in space the equation to locally diagonalize the total XC potential is:

[TABLE]

Here, ${U}$ is a general transformation matrix to be determined. Substituting now the expression in Eq. (43), truncated at the first order, into the last equation, we obtain the complete expression needed to rotate the reference frame along the local direction where the XC potential is diagonal:

[TABLE]

In passing we note that, for general non-collinear XC functionals, the relation between diagonalizing $\rho$ and $V$ might be more complex. It is interesting to separate the terms belonging to the last expression, as a function of their order in $\lambda$:

[TABLE]

From Eq. (IV.1.2) we can extract two useful properties of the zeroth order transformation equation:

[TABLE]

We can rewrite the Eq. (47) using the unitary properties (42) and the previous relations (49) and (50) as follows:

[TABLE]

and after some algebraic manipulation we get:

[TABLE]

The last equation can be represented symbolically:

[TABLE]

As we can see, the first term on the left-hand side is a purely off-diagonal matrix which is needed to compensate the off-diagonal terms present in the second term to guarantee that the right-hand side term is diagonal. Interestingly, this implies that the first order transformation matrix $\hat{U}^{(1)}$ can never be neglected since its presence ensures the unitarity of the transformation. Determining the elements of $\hat{U}^{(1)}$ requires some algebraic effort, which can be circumvented by exploiting the analogous of Eq. (IV.1.2) for the density matrix $\hat{\rho}$:

[TABLE]

In fact, the same mathematical process from Eq. (43) to (IV.1.2) is valid for the density operator as well. In this way, it holds for both Eq.s (IV.1.2) and (IV.1.2) that the quantity

[TABLE]

can be written explicitly from the density matrix Eq. (IV.1.2) as:

[TABLE]

where we have defined

[TABLE]

The knowledge of this quantity fully solves the system of Eq.s (IV.1.2) and (IV.1.2), indeed, to obtain the transformed XC potential in the Cartesian system of reference we shift the second term from the left-hand to the right-hand side of Eq. (IV.1.2) and apply the inverse transformation. Similarly to Eq. (57), we can define an analogous expression for XC potential:

[TABLE]

Now, we can rewrite Eq. (IV.1.2) and find the expression for the perturbed potential in the local reference system:

[TABLE]

In this way, the final expression for the matrix elements is:

[TABLE]

and its expression in the global Cartesian reference requires to apply the inverse transformation to Eq. (60):

[TABLE]

This method and the expression for the perturbed XC potential is completely general and independent of its flavour (LSDA or GGA or others, for example).

It is possible, moreover, to find the analytical expressions of the $U^{(0)}$ and $U^{(1)}$ transformation matrices in terms of the Pauli matrices ${\boldsymbol{\sigma}_{\alpha\beta}}$, i. e., the generators of the Lie group SU(2). Specifically, we can write:

[TABLE]

where ${\hat{n}}$ denotes the rotation axis direction and $\theta$ the rotation angle (in the following, the $(0)$ superscript will be omitted for 0th order quantities). In particular, these two quantities can be written in terms of the magnetization components

[TABLE]

[TABLE]

thus defining the rotation operation from the local spin-coordinate frame (with the local $z$ axis aligned within the local magnetization direction) to the global spin-coordinate reference system. The use of Eq.s (62), (63) and (64) allows to obtain an analytical expression for the $U^{(1)}$ and $U^{\dagger(1)}$ matrices:

[TABLE]

[TABLE]

Here, ${\hat{n}}^{(1)}$ and $\theta^{(1)}$ are the derivatives of ${\hat{n}}$ and $\theta$, respectively, with respect to the perturbation and precisely:

[TABLE]

[TABLE]

Combining Eq.s (65), (66), (67), (68) and using the normalization conditions for the rotation axis components ($n_{x}^{2}+n_{y}^{2}=1$ and $n_{x}n_{x}^{(1)}+n_{y}n_{y}^{(1)}=0$), we recognize that:

[TABLE]

and

[TABLE]

The last two expressions constitute the final equations to treat the LSDA approximation of the XC functional. These equations, (69) and (70), have been implemented as an alternative to the method (i) presented at the previous subsection, and which we identify as (i’).

IV.2 (ii) Explicit evaluation of the first order XC potential

Alternatively to the expansion of the rotation matrix, an explicit expression of the first order XC potential can be obtained by performing a Taylor expansion of the right hand side of Eq. (15). Introducing the LSDA definitions of the XC electrostatic potential $v_{xc}$ and the XC magnetic field magnitude $B_{xc}$,

[TABLE]

one obtains the following expression for $V_{xc}^{(1)}$ relative to the global spin quantization axis (omitting again the subscript “$xc$”):

[TABLE]

where ${\bf I}$ denotes the identity matrix. In the above expression, the first term corresponds to the change of the direction of the XC magnetic field due to rotation of magnetization induced by an external perturbation (e.g. atomic displacements). The second term describes the change of electric scalar potential induced by changes of $n$ and $m$, while the last term stems from the magnitude change of $B_{xc}$ due to the variation of the electronic density $n$ and the magnetic moment magnitude $m$. Hence, it can be readily recognized that the combination of the second and the third terms corresponds to a rotation of the collinear LSDA potential using the $\hat{U}^{(0)}$ transformation matrix, while the first term, related to the change of magnetization direction, contains contributions involving the $\hat{U}^{(1)}$ transformation matrix introduced above. The second derivation is comparatively simple and uses more intuitive ingredients like the xc effective magnetic field.

Within the LSDA approximation, the spin diagonal $V^{(1)}_{xc}$ and $B^{(1)}_{xc}$ can be expressed as

[TABLE]

In these expressions, we use $\rho_{\uparrow}$ and $\rho_{\downarrow}$ which are defined in Eq.s (16) and (17). Although the expressions (61) and (72) should be equivalent, Eq. (61) is more generic since, in principle, it can be applied to any combination of $V_{xc}^{(0)}$ and $V_{xc}^{(1)}$.

In the following, we will test and compare the two methods to verify that we obtain the same result, and see if one method is more efficient than the other.

IV.3 Implementation for atomic and electric field perturbations

We have implemented the aforementioned treatments of DFPT within non-collinear magnetism in the abinit code Gonze et al. (2016) (available from version 8.8.3 onwards). This implementation is currently only available in the norm-conserving pseudopotential formalism; the generalization to the Projector Augmented-Wave (PAW) approach is in progress and should be straightforward: the spin rotation should be applied to the PAW on-site XC potential expressed in terms of spherical harmonics. In the case of an atomic displacement perturbation, one has to take care of the so-called “frozen” part of the second energy derivative (see Eq. (13) in Ref. Gonze and Lee (1997)), which is not affected by the non-collinearity of the density. Another delicate term is the so-called “non-linear exchange-correlation core correction” (NLCC)Bechstedt (2015), which is a correction to the XC energy to compensate the error caused by the frozen core approximation, and the separation between core and valence charge densities in the pseudopotential formalism, accounting for possible overlap between them. Even if this correction only concerns the electronic charge density $n$, it should be accounted for correctly because it is not invariant under spinor rotation transformations. The treatment of electric field perturbation involves the derivative with respect to the k-point (ddk), which does not involve any frozen part or NLCC, such that its generalization into non-collinear magnetism is straightforward. The implementation works also in the case where the spin-orbit coupling is included.

The different methods of magnetization rotation can be set through the ixcrot abinit input flag. The default value is method (i) (ixcrot=1). The method (ii) or (i’) can be used by setting ixcrot to 2 or 3 respectively.

V Test cases

V.1 Accuracy and convergence of the implementation: Cr2O3

In this section we show how the different methods to treat the non-collinear DFPT perform for phonons with the LDA XC functional. We start by performing tests on Cr2O3, which is a collinear antiferromagnetic insulator. Cr2O3 is also known because it was the first proposed and measured magnetoelectric crystal Dzyaloshinskii (1959); Astrov (1960, 1961) and became a reference system to study the microscopic origin of bulk magnetoelectricity from DFT Íñiguez (2008); Bousquet et al. (2011); Malashevich et al. (2012); Tillack et al. (2016), which requires non-collinear treatment of the magnetism (to go beyond the exchange-striction effectMostovoy et al. (2010)).

We perform our calculations with the abinit package and the Pseudodojo LDA norm-conserving pseudopotentials van Setten et al. (2018) (including spin-orbit coupling). Cr2O3 keeps its insulating state even without Hubbard correction for the d orbitals of the Cr atom. To test our implementation of the DFPT within non-collinear magnetism, we first compare the second energy derivatives (E${}^{(2)}=\frac{1}{2}\frac{\partial^{2}E}{\partial\lambda^{2}}$) obtained from DFPT and finite difference (FD) calculations for an atomic displacement perturbation in one direction (Cr atom along the $x$ direction displaced by four amplitudes - two positive and two negative: -2$\tau$, $-\tau$, $+\tau$ and $+2\tau$, with $\tau=0.003$ Å - ; analogous results hold for the other two displacement directions that are not shown for simplicity). This 5-point FD scheme allows to extract the pure harmonic contribution to the energy variation and thus guarantees the validity of the comparison with the linear response calculations at the second order. The results are shown in Tab. 1 where we report the comparison of the frozen contribution alone and the total energy second derivatives for both DFPT (within its three (i), (i’) and (ii) methods) and FD. In Tab. 1 we also include the number of self-consistent iterations (SCI) needed to reach a given residual of the potential to stop the SCI. We also put a maximum number of iterations of 80 such that if the SCI did not reach the required precision the calculation will stop and we report the potential residual that the calculation could reach at the 80${}^{\text{th}}$ step. As can be seen the three approaches agree to within machine precision, and the agreement with FD is within 5-7 significant digits. We also observe that within 80 SCI steps, the method (ii) does not reach the same precision as methods (i) or (i’) do. Method (ii) seems then to give slower converge. We also remark that method (i) converges with a bit less number of SCI steps than method (i’). We note that this slower convergence of method (ii) is observed even with changing the mixing or preconditioning parameters. However, this characteristics may vary depending on the case study and method (ii) could appear better in other crystals than Cr2O3.

In Tab.2 we report the calculated TO phonon frequencies of Cr2O3 for the collinear case and for the non-collinear (magnetic moments along the $z$ direction, the $x$ and $y$ directions give the same results and are not shown for clarity) case without spin-orbit coupling. We can see that the method (i), (i’) and (ii) give strictly the same phonon frequencies within less than 1 cm*-1*. The comparison is also extremely good with the collinear case, which is what we expect for a non-frustrated collinear antiferroagnet as Cr2O3. These very good agreement confirm the robustness of our implementation.

V.2 Strong spin-orbit interaction case: RuCl3

$\alpha$-RuCl3 is a layered material, with space group $P6_{3}/mcm$ and a honeycomb arrangement of the Ru atoms, where the 4$d$ electrons drive unusual electronic and magnetic behavior due to strong spin-orbit interaction Plumb et al. (2014). The magnetic anisotropy of RuCl3 is very strong, where the magnetic moments lie perpendicular to the stacking plane and where a magnetic field of about 14 T is necessary to flip the moments in-plane Majumder et al. (2015). The magnetic ground state is known to be the so-called ZZ antiferromagnetism (ZZ-AFM) but for our test purpose we will simulate the crystal within its ferromagnetic phase. As for Cr2O3, we performed our calculations with the abinit package and the Pseudodojo LDA norm-conserving pseudopotentials van Setten et al. (2018) including spin-orbit coupling. The calculations were done with a grid of 8$\times$8$\times$8 k-points and a cutoff energy of 50 Ha for the plane wave expansion.

We performed internal atomic coordinate relaxations at fixed cell parameters ($a=$6.561 and $c=$5.709 Å) for both collinear and non-collinear cases (with a residual on the forces of 2 10*-5* Ha/Bohr) before performing the DFPT calculations of the phonons (using method (i)). For the non-collinear simulations, we treated both in-plane ($x$-direction) and out-of-plane ($z$-direction) alignment of the magnetic moments. Since no Hubbard correction is done for the Ru-$d$ orbital correlations, we find that the ground state is metallic with or without the spin-orbit coupling (as reported in Ref. Majumder et al. (2015)). We obtain a magnetocrystalline anisotropy of about 780 $\mu$eV, which is consistent with previous calculations Iyikanat et al. (2018) and the strong anisotropy reported experimentally.

We report the results in Tab.3. We can observe that the differences between collinear and non-collinear cases are small. Frequency differences for different orientations are within $<2$ cm*-1*, except the $\Gamma_{3}^{-}$ modes. A deviation of 9 cm*-1* is observed for mode number 7 when spins are along the $z$ directions and of 3 cm*-1* when spins are along the $x$-direction (when comparing to the collinear calculation). A deviation of 9 cm*-1* and 4 -1 is observed for mode number 20 when spins are aligned along $x$ and $z$ respectively. We also remark a small frequency shift of 1 (spins along $x$ direction) and 3 cm*-1* (spins along $z$ direction) of mode number 21 with label $\Gamma_{1}^{+}$. The modes with $\Gamma_{3}^{-}$ label, which exhibit the strong non-collinear spin-phonon coupling, are the only two modes that involves motions of the Ru against each other along the $z$ direction. The associated Cl motions are along the in-plane direction only and are such that they come closer to Ru going away from their $xy$ plan and they go away from Ru going closer to their plan. This spin-orbit shift of phonon bands denotes an appreciable spin-electron-phonon coupling, which would require further studies to understand the exact microscopic origin. We also observe a small unstable mode for all the cases, showing that the FM phase at the fixed cell parameters we used is structurally unstable (which is not surprising as it is not the ground state).

The results on RuCl3 validate our implementation, in that performing explicitly non-collinear DFPT calculations is important for systems where the non-collinearity is strongly coupled to the lattice properties of the crystal.

VI Conclusion

In this paper we have shown how to treat non-collinear magnetism within the DFPT. The main problem is that most of the XC functionals, like LDA or GGA, are derived at the collinear level, thus missing the non-diagonal terms of the density containing the x and y components of the magnetization. This prevents to use them as such in DFPT, where perturbations can change the orientation of the magnetization. We have derived several possibilities to treat this problem, which all are based on reducing the non-collinear derivatives to their collinear ones by aligning the z coordinate axis to the local magnetization direction (before or after perturbation). At each point the derivative of the XC potential can thus be performed in the usual collinear framework, the knack being to perform the back and forth transforms between the local collinear system and the global non-collinear one. We propose two possibilities to perform these changes of coordinate system for the XC functional, either through a Taylor expansion of the rotation matrix between the two coordinate systems, or using an analytical expression of the XC functional (only done for the LDA here). These different methods produce identical energy results in the test cases we explored, with method (ii) being somewhat slower and harder to converge. In pathological cases it is possible that one method could be qualitatively better. We have done tests on the atomic displacement and electric field perturbations in Cr2O3 and shown that in RuCl3 with strong spin-orbit coupling and magnetic anisotropy the phonons are affected by the non-collinear spin directions.

The immediate perspectives of our work are to extend this formalism to the Projector Augmented-Wave (PAW) approach Audouze et al. (2006, 2008) and to GGA and hybrid XC functionals, and to combine it with the linear response of crystals to magnetic fieldsSavrasov (1998); Cao et al. (2018).

Acknowledgements.

• SP and FR contributed equally to this work. The authors would like to thank X. Gonze for the fruitful discussions. EB MJV and FR acknowledge the FRS-FNRS for support through PDR projects MaRePeThe (GA 19528980) and “Transport in novel VDW heterostructures” (GA T.1077.15-1/7). SP acknowledge the Marie Curie COFUND postdoctoral fellow at the University of Liege and with the support of the European Union. MV, FR, SP and EB acknowledge the ARC AIMED project, the PRACE project TheDeNoMo and the CECI facilities funded by F.R.S-FNRS (Grant No. 2.5020.1) and Tier-1 supercomputer of the Fédération Wallonie-Bruxelles funded by the Walloon Region (Grant No. 1117545). MT was granted access to the french HPC resources of TGCC under the allocation 2018-AP010910358 attributed by GENCI (Grand Equipement National de Calcul Intensif).

Appendix A Building the density matrix

A.1 DFPT: first order

In order to write the expression for the first order density matrix $\hat{\rho}^{(1)}$ we need to derive Eq. (18) with respect to a given perturbation $\lambda$:

[TABLE]

The generic element of this operator, projected on the $\alpha,\beta$ spin states, is:

[TABLE]

At this point, we can write the $\hat{\rho}^{(1)}$ general expression, in terms of the wavefuctions analogously to the GS (see Eq. (22)):

[TABLE]

The last equation has been obtained as a derivative of the Eq. (25) and, interestingly, we note that the shape of the density matrix does not change with respect a given order of perturbation. Following the analogous procedure used for the GS case (Eq. (28) and following), we obtain the expression of the charge and magnetization density components in terms of spinors:

[TABLE]

Separating the real and imaginary parts of the wavefunctions to rewrite the last equations we remark that all the quantities are real, as we expected:

[TABLE]

Appendix B Local potential integral

In addition to the perturbed density matrix elements, that have to be built for a non-collinear magnetic system, we need to specify the terms involved in this complicated formalism, as well as the so-called “local potential”. This term (see Eq. (91) in Gonze (1995)), for a given band, is defined as follows:

[TABLE]

The terms to be integrated then can be explicitly written performing the matrix product between spinors and the local potential matrix:

[TABLE]