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

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 an ab 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.


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 groundstate 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 manyelectron 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 chargeneutral 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 spinwave 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][10][11][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 2 Advances in Mathematical Physics 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 exchangecorrelation 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 manybody 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 initialvalue 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 exchangecorrelation kernel. To simplify the notation, Hartree atomic units are used throughout.

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 where the angular brackets indicate the expectation value with respect to the ground-state many-electron wave function in second quantization and is the commutator of two quantum-mechanical operators. The spin-raising and spin-lowering operators are given bŷ and bŷ respectively, in terms of the time-dependent annihilation operator̂ (r, ) for an electron with spin ∈ {↑, ↓} and its adjoint, the creation operator̂ † (r, ), in the Heisenberg picture. The fermionic field operators with equal time arguments satisfy the anticommutation relations [3] {̂ (r, ) , as well as 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 Θ( − ὔ ) in (1) reflects causality and ensures that the spin density at time is only influenced by the magnetic field at times ὔ < , if the dynamic spin susceptibility is interpreted as the linear spin-density response function by means of Kubo's formula.

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 The derivative of the spin-raising operator (3) is in turn given by For a general interacting many-electron system governed by the Hamiltonian Advances in Mathematical Physics

3
where ℎ(r) = −∇ 2 /2 + (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] 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 righthand 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 In fact, only the kinetic energy in ℎ(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̂ + (r, ) 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 ℎ(r) were replaced by another Hamiltonianh(r) with a different spin-independent local potential̃ (r) instead of (r). This will be turned to an advantage below. It should be noted that although the time evolution of the operator̂ + (r, ) 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 systemdependent 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 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 and do not depend on the time variable for a stationary system in the ground state. If all terms are collected, one thus obtains the equation of motion 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 ὔ ; − ὔ ) or other, known quantities. Second, an initial value from which the time propagation can start must be established.

Expansion in an Auxiliary Basis.
In order to achieve the first of the two outlined tasks, a set of single-particle orbitals (r) is chosen, with the sole requirement that they are eigenfunctionsh (r) (r) = (r) of a Hamiltonianh(r) = −∇ 2 /2 +̃ (r) with some spinindependent local potential̃ (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 aŝ † (r, ) = ∑ * (r)̂ † ( ) , (16a) where the annihilation operator̂ ( ) and the creation operator̂ † ( ) now refer to the chosen basis. Inserting these expressions into the definition of the transverse dynamic spin susceptibility (1) yields the representation with the coefficients 4

Advances in Mathematical Physics
Crucially, if one exploits the relations then the hitherto elusive expectation value on the righthand side of (14) can be rewritten in terms of the same coefficients. With the expansion of the one-particle reduced density matrix and the resolution of the identity 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 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 ὔ ; − ὔ ) = 0 and hence The discontinuity at = ὔ is described by the second term on the right-hand side of (22), which establishes that the limit from the above equals The same result can also be directly derived from (1), which implies 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 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 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 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 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 Hamiltoniañ ℎ(r) with a spin-independent local potential̃ (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-densityfunctional theory is ruled out because, for collinear magnetic systems, the effective potential contains an explicit spindependent 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 manybody perturbation theory, using suitable approximations [19].

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 +− ( − ὔ ) an accurate representation of the true result with a manageable, finite number of basis functions and in constructing the initial value, specifically the spinresolved 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.

Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this article.