Many magnetic materials do not conform to the (anti-)ferromagnetic paradigm where all electronic spins are aligned to a global magnetization axis. Unfortunately, most electronic structure methods cannot describe such materials with noncollinear electron spin on account of formally requiring spin alignment. To overcome this limitation, it is necessary to generalize electronic structure methods and allow each electron spin to rotate freely. Here, we report the development of an *ab initio* time-dependent non-relativistic two-component spinor (TDN2C), which is a generalization of the time-dependent Hartree-Fock equations. Propagating the TDN2C equations in the time domain allows for the first-principles description of spin dynamics. A numerical tool based on the Hirshfeld partitioning scheme is developed to analyze the time-dependent spin magnetization. In this work, we also introduce the coupling between electron spin and a homogenous magnetic field into the TDN2C framework to simulate the response of the electronic spin degrees of freedom to an external magnetic field. This is illustrated for several model systems, including the spin-frustrated Li_{3} molecule. Exact agreement is found between numerical and analytic results for Larmor precession of hydrogen and lithium atoms. The TDN2C method paves the way for the *ab initio* description of molecular spin transport and spintronics in the time domain.

## I. INTRODUCTION

Spin-based technology takes on an ever-increasing role in material design as scientists and engineers seek to utilize and manipulate spins to maximize the quantum yield of solar cells, to process quantum information, and to resolve complex molecular structure. Fundamental to these scientific and technological applications is a detailed understanding of spin-dependent many-electron dynamics.^{1–4} To date, most theoretical investigations on molecular spintronics rely on the use of effective model Hamiltonians which depend on experimental parameterization^{5–8} and focus on collinear spin systems in which all the spin magnetizations are aligned along a single anisotropy axis. While these methods have their merits, they are limited in their description of the time-dependent nature of the spin-dependent many-electron wave function. This is especially true for phenomena driven by strong spin non-collinearity, spin-orbit coupling, or by interaction with intense electromagnetic fields. A proper description of spin dynamics must come from a solution of the first-principles spin-dependent Hamiltonian in the time domain.

Most standard electronic structure methods, however, cannot treat systems involving non-collinear spins because they formally require electronic spins to be aligned (anti-)parallel with respect to each other. In both the unrestricted Hartree-Fock (UHF) and spin-density density functional theories (SDFT), the one-electron spin orbitals are required to be eigenfunctions of the spin operator

**m**(

**r**) is of the same direction at all points in space.

In order to treat non-collinear spins, one needs to retain the full vector form of the magnetization **m**(**r**) and allow each spin quantization axis to rotate. This is equivalent to writing the spin orbitals as a superposition of the spin-up and spin-down states. For Hartree-Fock, this leads to the generalized Hartree-Fock (GHF) method,^{9–13} which is similar in structure to the wave function used in two-component relativistic models. For density functional theory (DFT), a generalization of the functionals to account for non-collinear spins is required.^{14–22} The GHF and non-collinear DFT methods have been used to study the non-collinear spin states of dissociated molecules and geometrically frustrated systems.^{10–13,15,16,19–21,23–26} In general, the phenomenon of non-collinear magnetism arises from two different types of interactions: the first comes from nearest-neighbor magnetic exchange interaction and is usually the dominant contribution; the second stems from the anisotropic interaction due to spin-orbital coupling and is typically one order of magnitude weaker for light elements. The first kind of interaction, i.e., the magnetic exchange, comes from the Pauli principle which is embedded in the antisymmetrized electronic wave function, even though in the non-relativistic Hamiltonian the Coulomb repulsion does not explicitly depend on spin. Therefore, a non-collinear non-relativistic treatment is still necessary for a proper description of spin magnetizations, especially in cases of spin frustration. It should be noted that although the non-relativistic Hamiltonian commutes with the spin operators *S*^{2} and *S*_{z} and the exact wave function may possess all the spin symmetries, approximate wave functions, such as a single Slater determinant, need not have these symmetry constraints.^{9–11,13,27} In fact, adding symmetry constraints can only raise the energy of a variationally optimized solution. Thus, by removing some symmetry constraints, it is possible to obtain a lower energy variational solution. In such a situation, one is forced to choose between an approximate wave function that is variationally optimum and a higher energy approximate wave function that contains the proper symmetries of the Hamiltonian. This dilemma was first pointed out by Löwdin, and is known as Löwdin's symmetry dilemma.^{28}

From the point of view of electronic dynamics, the extension of the collinear DFT and HF to the time domain gives rise to the time-dependent density functional theory (TDDFT) and time-dependent Hartree-Fock (TDHF).^{29–39} In particular, the real-time formalism of TDDFT/TDHF (RT-TDDFT/RT-TDHF) has proven to be a promising method for revealing the physical underpinnings of dynamical multi-electron phenomena.^{40–44} However, these methods are unable to capture the spin dynamics driven by an external magnetic field, or an internal magnetic field generated by the motion of other electrons within the system. Implicit to the single Slater determinant based RT-TDDFT/RT-TDHF method are the constraints that spin state of each electron cannot change during the time evolution and that they are collinear with respect to each other. These constraints prevent conventional RT-TDDFT/RT-TDHF from fully accounting for spin-dependent many-electron dynamics.

In order to describe non-collinear spin dynamics, we extend the generalized Hartree-Fock to the time domain, giving rise to the real-time time-dependent non-relativistic two-component spinor (RT-TDN2C) method. We present in this work the mathematical formalism and numerical implementation with case studies of non-collinear spin dynamics in response to an external magnetic field. To the best of our knowledge, this work is the first attempt to address many-body spin dynamics within the full *ab initio* description.

## II. METHODOLOGY

In the following derivations, the wave function is represented in a set of atomic orbitals (AO) {χ} with μ, ν, … as indices. Primed notations are used for matrices in the AO basis and unprimed notations for matrices in the orthonormal basis. In particular, **F**′ and **P**′ (**F** and **P**) are Fock and density matrices in the AO (orthonormal) basis. **F** and **P** are related to their AO basis counterparts through an orthonormal transformation. σ and τ are used as spin variables (σ ∈ {α, β}; τ ∈ {α, β}). All equations are expressed in atomic units.

### A. Time-dependent non-relativistic two-component spinor formalism

In the TDHF theory, the *N*-electron wave function is represented by a single Slater determinant ψ(*t*) composed of *N* time-dependent spin orbitals, {ψ_{k}(**x**, *t*)}

where the variable **x** includes both the spatial coordinate **r** and the spin coordinate ** ω** for one electron.

Applying the Dirac-Frenkel time-dependent variational principle^{45,46} to this approximate wave function, we obtain the equations of motion for the single-particle spin orbitals,

where

where

In order to allow for non-collinearity, we express the time-dependent spin orbitals as a time-dependent two-component spinor,

where the spatial functions

_{μ}(

**r**)},

where the time-dependent expansion coefficients *C*’s may take complex values. With the definition of the two-component spinor, the time-dependent single-particle density operator can be formulated with the following spin-blocked structure:

In the AO basis, its matrix form can be written as

where each spin block **P**′^{στ}(*t*) is given by

The Fock matrix in the AO basis also takes a spin-blocked form analogous to the density matrix

where each spin block **F**′^{στ}(*t*) is given by

where **h**′^{στ} is the one-electron Hamiltonian matrix, whose time-dependence arises from the interaction of the spin system with external electromagnetic perturbations (to be discussed in Sec. II C). **J**′ and **K**′ are Coulomb and exchange matrices, respectively,

In this work, we use a set of atom-centered atomic orbital basis functions which form a non-orthogonal basis set. The density and Fock matrices are then transformed from the AO basis to an orthonormal basis using the following equations:

The transformation matrix **V** = **S**^{1/2} in the Löwdin orthonormalization method, or may be obtained by the Cholesky decomposition **S** = **V**^{T}**V**, where **S** is the overlap matrix of AO basis functions, *S*_{μν} = ⟨μ|ν⟩.

The final density-matrix based TDN2C Liouville-von Neumann equation in the orthonormal basis can be written as

Equation (17) is the working equation to describe the time evolution of the non-relativistic two-component spinor. It can be considered as a time-dependent extension of the generalized Hartree-Fock equation. As with most time-dependent electronic structure theories, Eq. (17) can be solved in a response formalism in the frequency domain or propagated explicitly in the time-domain. In Secs. II B–II D, we will present a time-propagation approach and perturbation terms arising from the interaction between a spin system and a static magnetic field, as well as a method designed to analyze spin magnetization in the time domain.

### B. Numerical propagation using modified midpoint and unitary transformation

To obtain a time-evolution dynamics of a spin system, the time-dependent equation for non-relativistic two-component spinors (Eq. (17)) is integrated with a modified midpoint and unitary transformation (MMUT) algorithm.^{40,42} In the MMUT method, the time-evolution operator is a unitary transformation matrix **U**(*t*_{n}) that is constructed from the eigenvectors **C**(*t*_{n}) and eigenvalues ** ε**(

*t*

_{n}) of the generalized Fock matrix at time

*t*

_{n},

where Δ*t* is the time step. The density matrix is then propagated from time *t*_{n−1} to *t*_{n+1} using the time-evolution operator **U**(*t*_{n}),

Note that the Fock and density matrices in Eqs. (18)–(20) are assumed to have a spin-blocked structure (Eqs. (15) and (16)). By computing the Fock matrix at the midpoint of the step, the MMUT method accounts for linear changes in the density matrix during the time step. The MMUT algorithm can be used with a large step size, while still maintaining very good control of numerical noise and integration errors.^{40}

### C. Two-component spinor in the presence of a static and uniform magnetic field

The formalism of the non-relativistic two-component spinor allows realistic simulations of non-collinear spin dynamics within the *ab initio* framework. In the non-relativistic framework, the interaction of an electron spin with external electromagnetic field is described by the so-called Schrödinger-Pauli Hamiltonian

where **A** and *U* are the vector potential and scalar potential of the electromagnetic field. **p** = −*i*** ∇** is the canonical momentum operator.

Embedding the Schrödinger-Pauli Hamiltonian into the time-dependent Hartree-Fock equation gives rise to the spin-dependent non-relativistic two-component spinor, i.e., the time-dependent Pauli-Hartree-Fock (PHF) equation

In this work, we consider a many-electron spin system interacting with an external static, uniform magnetic field **B**. Upon choice of the Coulomb gauge ** ∇** ·

**A**= 0 and

where

**l**=

**r**×

**p**is the orbital angular momentum. The last two terms represent the interaction of the orbital angular momentum with the magnetic field and the higher order term that is quadratic in the field strength. These two terms will affect the time evolution of the spatial distribution of the electrons and thus affect the total energy of the system. For weak magnetic fields (e.g., <10

^{−4}a.u.), the contributions of the last two terms are small and will be neglected here. Thus, we only consider the interaction of the electron spin with the magnetic field

** σ** are the Pauli matrices, and their action on the spin eigenfunctions gives

In the time-dependent non-relativistic two-component formalism, Eq. (27) enters the Fock matrix in Eq. (17). The effect of the interaction between the spin system and external homogeneous magnetic field will drive the spin dynamics propagated by Eqs. (18)–(20).

The effective spin Hamiltonian approach can be considered as a result of taking the non-relativistic limit of the Dirac equation for a spin in a homogeneous magnetic field. This is a simple step to test the time-dependent non-relativistic two-component formalism and algorithms presented here as the analytical solution is readily obtainable. The prospectus of including other perturbative terms will be discussed in Sec. IV.

### D. Spin magnetization analysis

The time evolution of the electron density *n*(**r**, *t*) and the spin magnetization density vector **m**(**r**, *t*) are defined through the reduced density matrix γ(*t*),

The full vector form of the spin magnetization density has been widely used in ground-state non-collinear DFT studies.^{19,20,22,23,47,48} In the AO basis, the electron density and magnetization density are expressed as

^{49}In this approach, one defines a pro-molecular or fuzzy magnetization density

where *A* gives the index of the atomic centers. The atomic contribution of the total magnetization density, i.e., the atomic magnetic moments, **m**_{A}(*t*) can be calculated through the weight function *w*_{A}(**r**),

The resulting time-dependent atomic magnetization moment **m**_{A}(*t*) provides a convenient way for analyzing the time-evolution of spin dynamics in a many-body system.

## III. BENCHMARK AND DISCUSSION

The Hirshfeld atomic moment is analyzed on-the-fly using the time-dependent magnetization densities. In this work, we choose three test cases to illustrate the spin dynamics within the framework of *ab initio* time-dependent non-relativistic two-component spinor: (1) one electron with one unpaired spin, (2) many electrons with one unpaired spin, and (3) many electrons with three non-collinear spins. Integration of the time-dependent non-relativistic two-component equation with the MMUT algorithm is implemented in the development version of GAUSSIAN program.^{50} The integration step size of 0.05 a.u. (∼1.21 as) is used for all test cases in this paper.

### A. Hydrogen atom in a uniform magnetic field

Evolution of a quantum mechanical spin-

_{0}= |

**B**|.

The hydrogen atom with a minimum STO-3G basis set provides an ideal test case for a simple two spin-state problem. In our simulation, the initial electron spin was aligned along *x* axis of the reference frame. A magnetic field was applied along *z* direction with a constant magnitude of 8.5 × 10^{−5} a.u., corresponding to ∼20 T. Figures 1(a) and 1(b) show the time evolution of the spin magnetization. As time evolves, the electron spin gradually changes its direction and exhibits a perfect precession about the applied magnetic field. Figure 1(c) displays the evolution of the *x*, *y*, and *z* components of the spin magnetization. As seen from Figure 1(c), time evolutions of spin magnetization (i.e., Larmor type precession) exhibit sinusoidal behavior, with an angular frequency ω = 8.5 × 10^{−5} a.u. This result is in excellent agreement with the analytic Larmor frequency ω_{0} = 8.5 × 10^{−5} a.u. In the absence of any spin dephasing mechanism, the Larmor type of precession of spin magnetization persists without any decay and the simulation also nicely demonstrates such a phenomenon.

### B. Lithium atom in a uniform magnetic field

In this section, we use lithium atom as an example to illustrate spin dynamics of a multi-electron system. The 3-21G basis set is used in this study to test the stability of the two-component spinor dynamics with a slightly larger basis set. The initial spin magnetization is aligned along the *y* axis of the reference frame, with the two 1*s* spins pointing in opposite directions with respect to each other, and the single 2*s* spin pointing in the *y* direction. The total spin magnetization is 1/2. A homogenous magnetic field is applied with a polar angle θ = 45° about *z* axis, and an azimuthal angle ϕ = 45° about *x* axis. The field strength was again chosen to be 8.5 × 10^{−5} a.u. (∼20 T). Figure 2(a) shows the time evolution of the total spin magnetization. As can be seen in Figure 2(a), the total spin magnetization precesses about the applied magnetic field with a constant angle and modulus, indicating the collinearity of the spins during time evolution. To have a better visualization of the precession, the total spin magnetization is projected on to the plane perpendicular to the applied magnetic field, referred to as the

*x*,

*y*, and

*z*components of the spin magnetization with respect to the reference frame. Sinusoidal oscillations of the spin magnetization were again observed. The computed angular frequency ω = 8.5 × 10

^{−5}a.u. is consistent with the analytic Larmor frequency.

### C. Non-collinear Li_{3} trimer in a uniform magnetic field

As a prototypical application of the time-dependent non-relativistic two-component spinor method, we investigate the spin dynamics of a two-dimensional lattice model with three lithium atoms forming an equilateral triangle. Such triangular systems are known to exhibit geometrical spin frustration: when the two spins at the first two lattice nodes are aligned anti-parallel, the spin at the third node cannot simultaneously minimize its interactions with the other two. These triangular systems may also possess a non-collinear antiferromagnetic Néel state, where the spin at each lattice node is oriented at angle of 120° with respect to each other. Experimentally, spin waves have been measured through neutron scattering experiments, as well as spin reorientation events in response to an applied external magnetic field. The spin dynamics of these systems have been studied with parameterized spin models such as the Heisenberg Hamiltonian.^{5–8} The *ab initio* time-dependent non-relativistic two-component spinor method offers a step towards the first-principles simulation of systems with non-collinear spins, and opportunities to obtain the insights into the complex nature of the spin dynamics.

The model system tested herein is an equilateral triangle consists of three Li atoms with the Li–Li bond length of 2.10 Å. The initial wave function at *t* = 0 is obtained by a self-consistent-field procedure using the generalized Hartree-Fock method and the 3-21G basis set. The GHF method properly identified the non-collinearity of the spin magnetization in the Li_{3} trimer, as shown in Figure 3(a). With this initial spin arrangement, at *t* = 0 a uniform magnetic field oriented perpendicularly to the trimer plane, with the same field strength of 8.5 × 10^{−5} a.u. (∼20 T) as used in previous simulations. Figure 3(b) illustrates the time evolution of the spin magnetization at each lattice point. Figure 3(b) shows that the magnetization of each individual Li atom precesses about the magnetic field, with the same angular frequency ω = 8.5 × 10^{−5} a.u. Due to the effect of the external magnetic field, the atomic magnetization vector on each node switches its direction in ∼1.0 ps, however, the overall non-collinearity is maintained throughout the simulation. This is not surprising since we are analyzing the atomic magnetizations in the absence of the spin-orbit and spin-spin coupling terms in our simulation.

## IV. CONCLUSIONS

In this paper, we have formulated an *ab initio* time-dependent non-relativistic two-component spinor theory, to study the non-collinear spin dynamics of many-electron systems in response to an external magnetic field. We employed a direct integration of the TDN2C Hartree-Fock equation using atom-centered basis functions and a unitary propagation approach with a modified midpoint algorithm. An analysis tool based on the Hirshfeld partitioning scheme has been developed to analyze the time-dependent spin magnetization.

For the simple one-electron system of the hydrogen atom and collinear multi-electron system of the lithium atom, the RT-TDN2C simulations yield the same results as the analytical Larmor precessions. As an important application of our methodology, we have simulated spin dynamics of a non-collinear Li_{3} trimer in response to an external magnetic field. The switching of the magnetization at each lattice node was observed during the dynamical simulation.

Presented in this work is a framework of *ab initio* time-dependent method that can be viewed as the bridge that connects the non-relativistic Schrödinger equation and the fully relativistic Dirac equation. The “spin-relaxed” approximation in the generalized equations is a step towards relativistic electronic structure theory, as the relativistic two-component methods can be formulated in the same way as the RT-TDN2C equation. This framework provides us with a unique opportunity to study spin dynamics from the first-principles perspective. As illustrated in this paper, the coupling between a many-spin system and a homogenous magnetic field can be introduced into this framework as a perturbation by taking the non-relativistic limit of the Dirac equation. Although including this coupling term in the time-dependent non-relativistic two-component spinor model does not offer much more than what can be learned from the parameterized models (e.g., Heisenberg Hamiltonian), the test cases presented in this work calibrate the robustness of the numerical implementation.

The non-relativistic spin-field coupling developed in this paper allows for electrons to respond to external magnetic fields by relaxing the spin constraints, but the TDN2C model still lacks the operators necessary to describe various spin angular momentum couplings, such as the spin-orbit and spin-spin couplings. These additional terms underlie many important experimental observations, such as spin-crossover and spontaneous magnetization. The incorporation of spin-orbit and spin-spin couplings in a non-relativistic time-dependent formalism requires additional theoretical scrutiny and will be pursued in near future.

The time-dependent non-relativistic two-component spinor formalism introduced herein lacks electron correlation, a problem inherited from the HF formalism. In a forthcoming paper, we will examine the correlated non-relativistic two-component method using non-collinear DFT functionals.

## ACKNOWLEDGMENTS

Financial support from the U.S. National Science Foundation (NSF) (CHE-1265945 supporting F.D., Graduate Research Fellowship DGE-1256082 to J.J.G.), the (U.S.) Department of Energy (DOE) (DE-SC0006863 supporting X.L.), and additional support from Gaussian, Inc., and the University of Washington Student Technology Fund are gratefully acknowledged.