# Exact Formulation of the Transverse Dynamic Spin Susceptibility as an Initial-Value Problem

Authors: Arno Schindlmayr
Journal: Advances in Mathematical Physics
 (2018)

---

## Abstract

The transverse dynamic spin susceptibility is a correlation function that yields exact information about spin excitations in systems with a collinear magnetic ground state, including collective spin-wave modes. In anab initio context, it may be calculated within many-body perturbation theory or time-dependent density-functional theory, but the quantitative accuracy is currently limited by the available functionals for exchange and correlation in dynamically evolving systems. To circumvent this limitation, the spin susceptibility is here alternatively formulated as the solution of an initial-value problem. In this way, the challenge of accurately describing exchange and correlation in many-electron systems is shifted to the stationary initial state, which is much better understood. The proposed scheme further requires the choice of an auxiliary basis set, which determines the speed of convergence but always allows systematic convergence in practical implementations.


---

## Body

## 1. Introduction

Accurately predicting the excitation spectra of interacting quantum-mechanical many-electron systems from first principles is a notoriously difficult problem, because the variational principle, which underlies such successful ground-state schemes as density-functional theory [1] or quantum Monte Carlo methods [2], cannot be exploited in this case. Furthermore, strategies based on an expansion of the many-electron wave function are limited to very small systems due to the rapidly increasing number of possible configurations. Under these circumstances, the most successful computational approaches currently rely on correlation functions, which contain the complete and, in principle, exact information about specific classes of excited states and can be directly related to experimental spectroscopies. Well-known examples include the single-particle Green function, which describes the propagation of individual quasiparticles, that is, injected electrons or holes, and the dynamic charge susceptibility or linear density-response function, which gives information about plasmons, excitons, and other charge-neutral optical excitations [3]. Another example, which will be the focus of this paper, is the transverse dynamic spin susceptibility, which describes spin-flip excitations, including single-particle Stoner excitations as well as collective spin-wave modes, in systems with a collinear magnetic ground state. Among the properties that can be directly obtained from these correlation functions are the energies, lifetimes, and oscillator strengths of the excitations.In the context of first-principles calculations, the transverse dynamic spin susceptibility may be constructed using either of two methods. Within many-body perturbation theory, it can be expanded in terms of single-particle Green functions and a kernel that couples the dynamics of electrons and holes with opposite spin [4, 5]. In actual implementations, the latter must be drastically simplified and is usually replaced by the static zero-frequency limit of the screened Coulomb interaction, often in combination with additional approximations, such as a Yukawa-potential model [6] or the neglect of interatomic screening in a localized Wannier basis [7, 8]. Alternatively, in time-dependent density-functional theory, the spin susceptibility of independent electrons and holes is renormalized with a dynamic exchange-correlation kernel [9–12]. The former can be expressed in terms of the Kohn-Sham orbitals and eigenvalues; accurate static density functionals exist for this purpose. The functional form of the latter is much more elusive, however. In practice, adiabatic kernels derived from static density functionals are almost exclusively used, but this choice, necessitated by a lack of viable alternatives, entails significant errors. For example, the local-density approximation yields the exact exchange-correlation potential and hence the exact Kohn-Sham orbitals and eigenvalues, for a homogeneous electron gas, but the adiabatic kernel derived from it neglects both the spatial and temporal nonlocality of the exact exchange-correlation kernel [13]. In both frameworks, therefore, the lack of accurate models for time-dependent exchange and correlation effects in nonstationary systems currently constitutes the main impediment for systematic quantitative improvements. Moreover, inconsistencies between the prevalent static and relatively crude approximations for dynamic two-particle correlation on the one hand and the accurate description of many-body effects on the quasiparticle electronic structure of stationary systems on the other hand give rise to qualitative errors, such as the numerical violation of Goldstone’s theorem for the spin-wave dispersion in implementations of many-body perturbation theory [7, 14].In this paper, I develop a different approach that holds the promise of circumventing the above-mentioned problems. Instead of calculating the transverse dynamic spin susceptibility from a Dyson-type integral equation with a dynamic exchange-correlation kernel, as is usually done, the key idea advanced here is to reformulate it as the solution of an initial-value problem in the time domain. To this end, I derive an exact differential equation that is of first order in the time variable, using the equations of motion of the field operators in second quantization. Furthermore, I show that the initial value can be expressed in terms of the spin-resolved one-particle reduced density matrix. As a consequence, this scheme requires an accurate description of nonlocal correlation in the stationary ground state, which is achievable even without explicit wave functions, but no dynamic exchange-correlation kernel. To simplify the notation, Hartree atomic units are used throughout.


## 2. Derivation

### 2.1. Definition of the Spin Susceptibility

In collinear magnetic systems, the direction of the magnetization vector is identical at all points in space and varies only in magnitude and sign. In such circumstances, spin is a good quantum number, as each electron has a definite spin orientation with respect to the global quantization axis. The retarded transverse dynamic spin susceptibility is then defined as(1)χ+-r,r′;t-t′=iS^+r,t,S^-r′,t′Θt-t′,where the angular brackets indicate the expectation value with respect to the ground-state many-electron wave function in second quantization and(2)S^+r,t,S^-r′,t′=S^+r,tS^-r′,t′-S^-r′,t′S^+r,tis the commutator of two quantum-mechanical operators. The spin-raising and spin-lowering operators are given by(3)S^+r,t=ψ^↑†r,tψ^↓r,tand by(4)S^-r′,t′=ψ^↓†r′,t′ψ^↑r′,t′,respectively, in terms of the time-dependent annihilation operator ψ^σ(r,t) for an electron with spin σ∈{↑,↓} and its adjoint, the creation operator ψ^σ†(r,t), in the Heisenberg picture. The fermionic field operators with equal time arguments satisfy the anticommutation relations [3](5)ψ^σr,t,ψ^σ′†r′,t=δσσ′δr-r′as well as(6)ψ^σr,t,ψ^σ′r′,t=ψ^σ†r,t,ψ^σ′†r′,t=0,where the curly brackets denote the anticommutator, defined analogous to (2) but with a plus sign instead of a minus sign on the right-hand side. Finally, the Heaviside step function Θ(t-t′) in (1) reflects causality and ensures that the spin density at time t is only influenced by the magnetic field at times t′<t, if the dynamic spin susceptibility is interpreted as the linear spin-density response function by means of Kubo’s formula.


### 2.2. Equation of Motion

The time derivative of the retarded transverse dynamic spin susceptibility, which will eventually lead to the desired differential equation, follows immediately from definition (1) and is given by(7)∂∂tχ+-r,r′;t-t′=i∂∂tS^+r,t,S^-r′,t′Θt-t′+iS^+r,t,S^-r′,tδt-t′.The derivative of the spin-raising operator (3) is in turn given by(8)i∂∂tS^+r,t=ψ^↑†r,ti∂∂tψ^↓r,t--i∂∂tψ^↑†r,tψ^↓r,t.For a general interacting many-electron system governed by the Hamiltonian(9)H^=∑σ∫ψ^σ†r,thrψ^σr,td3r+12∑σ,σ′∬ψ^σ†r,tψ^σ′†r′,tvr-r′×ψ^σ′r′,tψ^σr,td3r′d3r, where h(r)=-∇2/2+V(r) comprises the kinetic energy and the ionic potential and v(r-r′) is the pairwise Coulomb interaction, the field operators obey the equations of motion [3](10a)i∂∂tψ^↓r,t=hrψ^↓r,t+∑σ′∫ψ^σ′†r′,tvr-r′ψ^σ′r′,tψ^↓r,td3r′,(10b)-i∂∂tψ^↑†r,t=hrψ^↑†r,t+∑σ′∫ψ^↑†r,tψ^σ′†r′,tvr-r′ψ^σ′r′,td3r′.When applied to the single-particle Green function, it is well known that this procedure generates an infinite hierarchy of equations of motion where each Green function is coupled to a higher-order Green function due to the occurrence of terms involving products of three field operators on the right-hand sides of the above equations [3]. However, if (10a) and (10b) are inserted into (8), one finds that the terms involving the Coulomb interaction cancel in this case, leading to the remarkably simple result(11)i∂∂tS^+r,t=ψ^↑†r,thrψ^↓r,t-hrψ^↑†r,tψ^↓r,t.In fact, only the kinetic energy in h(r) yields a nonvanishing contribution, whereas the external potential commutes with the field operators and has no influence. As a consequence, the time evolution of S^+(r,t) turns out to be to some degree universal, as long as there are no spin-dependent magnetic fields that couple directly to the spin density and give rise to extra terms in (11). In particular, it would be identical if h(r) were replaced by another Hamiltonian h~(r) with a different spin-independent local potential V~(r) instead of V(r). This will be turned to an advantage below. It should be noted that although the time evolution of the operator S^+(r,t) is to some degree independent of the specific system, physical observables derived from it are not, of course, because the expectation values are still formed with a specific system-dependent wave function. The fact that the time derivative of an operator has a generic value is also not unusual in itself; for example, all operators that commute with the Hamiltonian have a vanishing time derivative, irrespective of the specific system. Another example of practical relevance is the continuity equation for the spin-density operator, to which the present result is related.Turning next to the second term on the right-hand side of (7), the commutator can be simplified to(12)S^+r,t,S^-r′,t=ψ^↑†r,tψ^↑r′,tδr-r′-ψ^↓†r′,tψ^↓r,tδr-r′by exploiting the anticommutation relation (5). The expectation values of the operators on the right-hand side correspond to the components of the spin-resolved one-particle reduced density matrix(13)nσr,r′=ψ^σ†r,tψ^σr′,tand do not depend on the time variable t for a stationary system in the ground state.If all terms are collected, one thus obtains the equation of motion(14)∂∂tχ+-r,r′;t-t′=ψ^↑†r,th~rψ^↓r,t-h~rψ^↑†r,tψ^↓r,t,S^-r′,t′Θt-t′+in↑r,r′-n↓r′,rδr-r′δt-t′,which takes the form of a differential equation that is of first order in the time variable. In order to turn it into a practical scheme for calculating the transverse dynamic spin susceptibility, two further problems must still be addressed, however: First, the expectation value involving products of four field operators on the right-hand side must be expressed in terms of χ+-(r,r′;t-t′) or other, known quantities. Second, an initial value from which the time propagation can start must be established.


### 2.3. Expansion in an Auxiliary Basis

In order to achieve the first of the two outlined tasks, a set of single-particle orbitalsφq(r) is chosen, with the sole requirement that they are eigenfunctions(15)h~rφqr=ϵqφqrof a Hamiltonian h~(r)=-∇2/2+V~(r) with some spin-independent local potential V~(r). This ensures that the orbitals form a complete orthonormal set, so that other quantities can be expanded in this basis. In particular, the fermionic field operators can be written as(16a)ψ^σ†r,t=∑pφp∗rc^pσ†t,(16b)ψ^σr,t=∑qφqrc^qσt, where the annihilation operator c^qσ(t) and the creation operator c^pσ†(t) now refer to the chosen basis. Inserting these expressions into the definition of the transverse dynamic spin susceptibility (1) yields the representation(17)χ+-r,r′;t-t′=∑p,q,r,sφp∗rφqrχpqrs+-t-t′φr∗r′φsr′with the coefficients(18)χpqrs+-t-t′=ic^p↑†tc^q↓t,c^r↓†t′c^s↑t′Θt-t′. Crucially, if one exploits the relations(19a)h~rψ^σ†r,t=∑pϵpφp∗rc^pσ†t,(19b)h~rψ^σr,t=∑qϵqφqrc^qσt,then the hitherto elusive expectation value on the right-hand side of (14) can be rewritten in terms of the same coefficients. With the expansion of the one-particle reduced density matrix(20)nσr,r′=∑p,sφp∗rnpsσφsr′and the resolution of the identity(21)δr-r′=∑q,rφqrφr∗r′, all terms in (14) can thus be cast into an identical form and written as sums over products of four orbitals. Although the products of the orbitals exhibit linear degeneracies and do not form an orthogonal basis themselves, a sufficient condition for the equation of motion to be fulfilled is that the coefficients inside the sums on both sides of the equation agree, which eventually yields(22)∂∂tχpqrs+-t-t′=-iϵq-ϵpχpqrs+-t-t′+inps↑-nrq↓δt-t′.


### 2.4. Initial Value

From the Heaviside step function contained in the definition of the transverse dynamic spin susceptibility (1), it already follows that χ+-(r,r′;t-t′)=0 and hence(23)χpqrs+-t-t′=0fort<t′.The discontinuity at t=t′ is described by the second term on the right-hand side of (22), which establishes that the limit from the above equals(24)limt↘t′χpqrs+-t-t′=inps↑-nrq↓.The same result can also be directly derived from (1), which implies(25)limt↘t′χ+-r,r′;t-t′=iS^+r,t,S^-r′,t. It then follows by inserting the already established simplified expression for the commutator (12) in combination with the definition of the one-particle reduced density matrix (13) and its expansion in the chosen basis set, as described in the previous section.Starting from the initial value (24), the further forward evolution of the transverse dynamic spin susceptibility is continuous and governed by(26)∂∂tχpqrs+-t-t′=-iϵq-ϵpχpqrs+-t-t′fort>t′,which only retains the first term on the right-hand side of (22). Taken together, (23), (24), and (26) thus define the initial-value problem that lies at the center of this work. Due to the simple linear form of (26), it can in fact be integrated analytically, yielding the explicit solution(27)χpqrs+-t-t′=inps↑-nrq↓e-iϵq-ϵpt-t′Θt-t′.Due to the inclusion of the Heaviside step function, which ensures that the solution fulfills (23) and thus respects causality, this formula is valid for all values of the time argument. Its Fourier transformation to frequency space can also be carried out analytically and is given by(28)χpqrs+-ω+iη=∫-∞+∞χpqrs+-τeiω+iητdτ=nps↑-nrq↓ω-ϵq-ϵp+iηwith η>0. The location of the poles in the lower half of the complex frequency plane is a characteristic feature of retarded (or causal) correlation functions.


### 2.5. Discussion

The formal solution of the initial-value problem derived in the previous section takes the form of an expansion in products of orbitals from a complete orthonormal basis, given by (17), for the spatial variation of the transverse dynamic spin susceptibility and a superposition of harmonic oscillations, given by (27), for its time dependence. Despite the very suggestive appearance, especially of the Fourier transform (28), it is important to recall that the parameters ϵq do not constitute any actual energy levels of the physical system in question but are instead the eigenvalues of a different, auxiliary problem (15) and thus serve a purely mathematical purpose. For the same reason, the sums in (17) also run over the entire set of basis functions and are not restricted to transitions between occupied and unoccupied states. Even for systems with a finite number of electrons, the sums are thus infinite, necessitating truncations in practical applications. The description of the time dependence as a superposition of an infinite number of harmonic contributions is akin to a Fourier representation.As pointed out above, the basis set can be chosen freely with the sole requirement that it corresponds to the eigenfunctions of an auxiliary single-particle Hamiltonianh~(r) with a spin-independent local potential V~(r). In practical implementations, where truncations of infinite sums are unavoidable, this degree of freedom may be exploited in order to achieve good convergence even with a finite number of terms. The auxiliary Hamiltonian must then be chosen in such a way that the expansion (20) of the one-particle reduced density matrix of the actual physical system is numerically accurate with a manageable, finite number of coefficients. Consequently, it must reflect the characteristics of the specific physical system. Obvious choices are not admissible, however: The noninteracting Kohn-Sham system of spin-density-functional theory is ruled out because, for collinear magnetic systems, the effective potential contains an explicit spin-dependent exchange-correlation field that forces the breaking of the spin symmetry [15], which arises spontaneously even without magnetic fields in the real physical system due to the Coulomb interaction. Likewise, the so-called natural orbitals [16] that diagonalize the one-particle reduced density matrix and thereby enable its most efficient representation cannot be employed; although the natural orbitals constitute a complete orthonormal set, they do not, in general, correspond to the eigenfunctions of a single-particle Hamiltonian with known eigenvalues, and for a collinear magnetic system they are also spin dependent. While generic basis sets like plane waves are always an option, of course, a weighted average of the local effective potentials of the two spin channels, or the self-consistent effective potential obtained from constrained density-functional theory with a non-spin-polarized electron configuration, might be preferable in actual implementations, because fewer basis functions are required for an accurate representation. In any case, although the choice of an optimal auxiliary basis set remains a practical challenge, it ultimately affects only the speed of convergence, whereas the converged results themselves are independent of the basis set.A key feature of the constructed initial-value problem is that the characteristics of exchange and correlation are contained in the initial value in form of the one-particle reduced density matrix. This is plausible from a fundamental point of view, because the dynamic spin susceptibility equals the linear spin-density response function, which is a property of the stationary ground state. As the one-particle reduced density matrix is not obtained within the present scheme, it must be constructed independently. For practical purposes, it is important that it can in fact be derived self-consistently without explicit recourse to many-electron wave functions by means of reduced density-matrix functional theory [17], a nonlocal extension of conventional density-functional theory. Although the design of reliable quantitative energy functionals is still very much in progress, reduced density-matrix functional theory is already a practically useful scheme [18]. Alternatively, it can be derived from the single-particle Green function, which is obtainable from implementations of many-body perturbation theory, using suitable approximations [19].


## 2.1. Definition of the Spin Susceptibility

In collinear magnetic systems, the direction of the magnetization vector is identical at all points in space and varies only in magnitude and sign. In such circumstances, spin is a good quantum number, as each electron has a definite spin orientation with respect to the global quantization axis. The retarded transverse dynamic spin susceptibility is then defined as(1)χ+-r,r′;t-t′=iS^+r,t,S^-r′,t′Θt-t′,where the angular brackets indicate the expectation value with respect to the ground-state many-electron wave function in second quantization and(2)S^+r,t,S^-r′,t′=S^+r,tS^-r′,t′-S^-r′,t′S^+r,tis the commutator of two quantum-mechanical operators. The spin-raising and spin-lowering operators are given by(3)S^+r,t=ψ^↑†r,tψ^↓r,tand by(4)S^-r′,t′=ψ^↓†r′,t′ψ^↑r′,t′,respectively, in terms of the time-dependent annihilation operator ψ^σ(r,t) for an electron with spin σ∈{↑,↓} and its adjoint, the creation operator ψ^σ†(r,t), in the Heisenberg picture. The fermionic field operators with equal time arguments satisfy the anticommutation relations [3](5)ψ^σr,t,ψ^σ′†r′,t=δσσ′δr-r′as well as(6)ψ^σr,t,ψ^σ′r′,t=ψ^σ†r,t,ψ^σ′†r′,t=0,where the curly brackets denote the anticommutator, defined analogous to (2) but with a plus sign instead of a minus sign on the right-hand side. Finally, the Heaviside step function Θ(t-t′) in (1) reflects causality and ensures that the spin density at time t is only influenced by the magnetic field at times t′<t, if the dynamic spin susceptibility is interpreted as the linear spin-density response function by means of Kubo’s formula.


## 2.2. Equation of Motion

The time derivative of the retarded transverse dynamic spin susceptibility, which will eventually lead to the desired differential equation, follows immediately from definition (1) and is given by(7)∂∂tχ+-r,r′;t-t′=i∂∂tS^+r,t,S^-r′,t′Θt-t′+iS^+r,t,S^-r′,tδt-t′.The derivative of the spin-raising operator (3) is in turn given by(8)i∂∂tS^+r,t=ψ^↑†r,ti∂∂tψ^↓r,t--i∂∂tψ^↑†r,tψ^↓r,t.For a general interacting many-electron system governed by the Hamiltonian(9)H^=∑σ∫ψ^σ†r,thrψ^σr,td3r+12∑σ,σ′∬ψ^σ†r,tψ^σ′†r′,tvr-r′×ψ^σ′r′,tψ^σr,td3r′d3r, where h(r)=-∇2/2+V(r) comprises the kinetic energy and the ionic potential and v(r-r′) is the pairwise Coulomb interaction, the field operators obey the equations of motion [3](10a)i∂∂tψ^↓r,t=hrψ^↓r,t+∑σ′∫ψ^σ′†r′,tvr-r′ψ^σ′r′,tψ^↓r,td3r′,(10b)-i∂∂tψ^↑†r,t=hrψ^↑†r,t+∑σ′∫ψ^↑†r,tψ^σ′†r′,tvr-r′ψ^σ′r′,td3r′.When applied to the single-particle Green function, it is well known that this procedure generates an infinite hierarchy of equations of motion where each Green function is coupled to a higher-order Green function due to the occurrence of terms involving products of three field operators on the right-hand sides of the above equations [3]. However, if (10a) and (10b) are inserted into (8), one finds that the terms involving the Coulomb interaction cancel in this case, leading to the remarkably simple result(11)i∂∂tS^+r,t=ψ^↑†r,thrψ^↓r,t-hrψ^↑†r,tψ^↓r,t.In fact, only the kinetic energy in h(r) yields a nonvanishing contribution, whereas the external potential commutes with the field operators and has no influence. As a consequence, the time evolution of S^+(r,t) turns out to be to some degree universal, as long as there are no spin-dependent magnetic fields that couple directly to the spin density and give rise to extra terms in (11). In particular, it would be identical if h(r) were replaced by another Hamiltonian h~(r) with a different spin-independent local potential V~(r) instead of V(r). This will be turned to an advantage below. It should be noted that although the time evolution of the operator S^+(r,t) is to some degree independent of the specific system, physical observables derived from it are not, of course, because the expectation values are still formed with a specific system-dependent wave function. The fact that the time derivative of an operator has a generic value is also not unusual in itself; for example, all operators that commute with the Hamiltonian have a vanishing time derivative, irrespective of the specific system. Another example of practical relevance is the continuity equation for the spin-density operator, to which the present result is related.Turning next to the second term on the right-hand side of (7), the commutator can be simplified to(12)S^+r,t,S^-r′,t=ψ^↑†r,tψ^↑r′,tδr-r′-ψ^↓†r′,tψ^↓r,tδr-r′by exploiting the anticommutation relation (5). The expectation values of the operators on the right-hand side correspond to the components of the spin-resolved one-particle reduced density matrix(13)nσr,r′=ψ^σ†r,tψ^σr′,tand do not depend on the time variable t for a stationary system in the ground state.If all terms are collected, one thus obtains the equation of motion(14)∂∂tχ+-r,r′;t-t′=ψ^↑†r,th~rψ^↓r,t-h~rψ^↑†r,tψ^↓r,t,S^-r′,t′Θt-t′+in↑r,r′-n↓r′,rδr-r′δt-t′,which takes the form of a differential equation that is of first order in the time variable. In order to turn it into a practical scheme for calculating the transverse dynamic spin susceptibility, two further problems must still be addressed, however: First, the expectation value involving products of four field operators on the right-hand side must be expressed in terms of χ+-(r,r′;t-t′) or other, known quantities. Second, an initial value from which the time propagation can start must be established.


## 2.3. Expansion in an Auxiliary Basis

In order to achieve the first of the two outlined tasks, a set of single-particle orbitalsφq(r) is chosen, with the sole requirement that they are eigenfunctions(15)h~rφqr=ϵqφqrof a Hamiltonian h~(r)=-∇2/2+V~(r) with some spin-independent local potential V~(r). This ensures that the orbitals form a complete orthonormal set, so that other quantities can be expanded in this basis. In particular, the fermionic field operators can be written as(16a)ψ^σ†r,t=∑pφp∗rc^pσ†t,(16b)ψ^σr,t=∑qφqrc^qσt, where the annihilation operator c^qσ(t) and the creation operator c^pσ†(t) now refer to the chosen basis. Inserting these expressions into the definition of the transverse dynamic spin susceptibility (1) yields the representation(17)χ+-r,r′;t-t′=∑p,q,r,sφp∗rφqrχpqrs+-t-t′φr∗r′φsr′with the coefficients(18)χpqrs+-t-t′=ic^p↑†tc^q↓t,c^r↓†t′c^s↑t′Θt-t′. Crucially, if one exploits the relations(19a)h~rψ^σ†r,t=∑pϵpφp∗rc^pσ†t,(19b)h~rψ^σr,t=∑qϵqφqrc^qσt,then the hitherto elusive expectation value on the right-hand side of (14) can be rewritten in terms of the same coefficients. With the expansion of the one-particle reduced density matrix(20)nσr,r′=∑p,sφp∗rnpsσφsr′and the resolution of the identity(21)δr-r′=∑q,rφqrφr∗r′, all terms in (14) can thus be cast into an identical form and written as sums over products of four orbitals. Although the products of the orbitals exhibit linear degeneracies and do not form an orthogonal basis themselves, a sufficient condition for the equation of motion to be fulfilled is that the coefficients inside the sums on both sides of the equation agree, which eventually yields(22)∂∂tχpqrs+-t-t′=-iϵq-ϵpχpqrs+-t-t′+inps↑-nrq↓δt-t′.


## 2.4. Initial Value

From the Heaviside step function contained in the definition of the transverse dynamic spin susceptibility (1), it already follows that χ+-(r,r′;t-t′)=0 and hence(23)χpqrs+-t-t′=0fort<t′.The discontinuity at t=t′ is described by the second term on the right-hand side of (22), which establishes that the limit from the above equals(24)limt↘t′χpqrs+-t-t′=inps↑-nrq↓.The same result can also be directly derived from (1), which implies(25)limt↘t′χ+-r,r′;t-t′=iS^+r,t,S^-r′,t. It then follows by inserting the already established simplified expression for the commutator (12) in combination with the definition of the one-particle reduced density matrix (13) and its expansion in the chosen basis set, as described in the previous section.Starting from the initial value (24), the further forward evolution of the transverse dynamic spin susceptibility is continuous and governed by(26)∂∂tχpqrs+-t-t′=-iϵq-ϵpχpqrs+-t-t′fort>t′,which only retains the first term on the right-hand side of (22). Taken together, (23), (24), and (26) thus define the initial-value problem that lies at the center of this work. Due to the simple linear form of (26), it can in fact be integrated analytically, yielding the explicit solution(27)χpqrs+-t-t′=inps↑-nrq↓e-iϵq-ϵpt-t′Θt-t′.Due to the inclusion of the Heaviside step function, which ensures that the solution fulfills (23) and thus respects causality, this formula is valid for all values of the time argument. Its Fourier transformation to frequency space can also be carried out analytically and is given by(28)χpqrs+-ω+iη=∫-∞+∞χpqrs+-τeiω+iητdτ=nps↑-nrq↓ω-ϵq-ϵp+iηwith η>0. The location of the poles in the lower half of the complex frequency plane is a characteristic feature of retarded (or causal) correlation functions.


## 2.5. Discussion

The formal solution of the initial-value problem derived in the previous section takes the form of an expansion in products of orbitals from a complete orthonormal basis, given by (17), for the spatial variation of the transverse dynamic spin susceptibility and a superposition of harmonic oscillations, given by (27), for its time dependence. Despite the very suggestive appearance, especially of the Fourier transform (28), it is important to recall that the parameters ϵq do not constitute any actual energy levels of the physical system in question but are instead the eigenvalues of a different, auxiliary problem (15) and thus serve a purely mathematical purpose. For the same reason, the sums in (17) also run over the entire set of basis functions and are not restricted to transitions between occupied and unoccupied states. Even for systems with a finite number of electrons, the sums are thus infinite, necessitating truncations in practical applications. The description of the time dependence as a superposition of an infinite number of harmonic contributions is akin to a Fourier representation.As pointed out above, the basis set can be chosen freely with the sole requirement that it corresponds to the eigenfunctions of an auxiliary single-particle Hamiltonianh~(r) with a spin-independent local potential V~(r). In practical implementations, where truncations of infinite sums are unavoidable, this degree of freedom may be exploited in order to achieve good convergence even with a finite number of terms. The auxiliary Hamiltonian must then be chosen in such a way that the expansion (20) of the one-particle reduced density matrix of the actual physical system is numerically accurate with a manageable, finite number of coefficients. Consequently, it must reflect the characteristics of the specific physical system. Obvious choices are not admissible, however: The noninteracting Kohn-Sham system of spin-density-functional theory is ruled out because, for collinear magnetic systems, the effective potential contains an explicit spin-dependent exchange-correlation field that forces the breaking of the spin symmetry [15], which arises spontaneously even without magnetic fields in the real physical system due to the Coulomb interaction. Likewise, the so-called natural orbitals [16] that diagonalize the one-particle reduced density matrix and thereby enable its most efficient representation cannot be employed; although the natural orbitals constitute a complete orthonormal set, they do not, in general, correspond to the eigenfunctions of a single-particle Hamiltonian with known eigenvalues, and for a collinear magnetic system they are also spin dependent. While generic basis sets like plane waves are always an option, of course, a weighted average of the local effective potentials of the two spin channels, or the self-consistent effective potential obtained from constrained density-functional theory with a non-spin-polarized electron configuration, might be preferable in actual implementations, because fewer basis functions are required for an accurate representation. In any case, although the choice of an optimal auxiliary basis set remains a practical challenge, it ultimately affects only the speed of convergence, whereas the converged results themselves are independent of the basis set.A key feature of the constructed initial-value problem is that the characteristics of exchange and correlation are contained in the initial value in form of the one-particle reduced density matrix. This is plausible from a fundamental point of view, because the dynamic spin susceptibility equals the linear spin-density response function, which is a property of the stationary ground state. As the one-particle reduced density matrix is not obtained within the present scheme, it must be constructed independently. For practical purposes, it is important that it can in fact be derived self-consistently without explicit recourse to many-electron wave functions by means of reduced density-matrix functional theory [17], a nonlocal extension of conventional density-functional theory. Although the design of reliable quantitative energy functionals is still very much in progress, reduced density-matrix functional theory is already a practically useful scheme [18]. Alternatively, it can be derived from the single-particle Green function, which is obtainable from implementations of many-body perturbation theory, using suitable approximations [19].


## 3. Conclusions

In conclusion, this paper takes a new perspective on the calculation of the transverse dynamical spin susceptibility for collinear magnetic systems, which seeks to circumvent current bottlenecks that prevent systematic improvements in the accuracy of quantitative first-principles simulations. Based on the notion that the spin susceptibility is a property of the stationary state and that its inherent time dependence does not arise from external fields but from the internal phases of wave functions of the unperturbed system, it is here reformulated as the solution of an initial-value problem. An essential point in this process is the expansion in an auxiliary orthonormal basis, which allows all occurrent expectation values to be written either in terms of the coefficients of the transverse dynamic spin susceptibility itself or in terms of other observables that can be obtained independently. The time propagation is simple enough to allow a formal analytic solution, while the problem of accounting for exchange and correlation effects is shifted to the initial value. The challenges hence lie in the identification of efficient auxiliary basis sets that makeχpqrs+-(t-t′) an accurate representation of the true result with a manageable, finite number of basis functions and in constructing the initial value, specifically the spin-resolved one-particle reduced density matrix, with sufficient quantitative accuracy. Both of these problems are better understood at a fundamental level and appear more accessible than the design of nonadiabatic exchange-correlation functionals for time-dependent magnetic systems needed for further improvements of prevalent implementations based on the solution of Dyson-type integral equations, however.



---

*Source: 3732892-2018-01-03.xml*