Ultrafast reorientation of the N\'eel vector in antiferromagnetic Dirac semimetals

Antiferromagnets exhibit distinctive characteristics such as ultrafast dynamics and robustness against perturbative fields, thereby attracting considerable interest in fundamental physics and technological applications. Recently, it was revealed that the N\'eel vector can be switched by a current-induced staggered (N\'eel) spin-orbit torque in antiferromagnets with the parity-time symmetry, and furthermore, a nonsymmorphic symmetry enables the control of Dirac fermions. However, the real-time dynamics of the magnetic and electronic structures remain largely unexplored. Here, we propose a theory of the ultrafast dynamics in antiferromagnetic Dirac semimetals and show that the N\'eel vector is rotated in the picosecond timescale by the terahertz-pulse-induced N\'eel spin-orbit torque and other torques originating from magnetic anisotropies. This reorientation accompanies the modulation of the mass of Dirac fermions and can be observed in real time by the magneto-optical effects. Our results provide a theoretical basis for emerging ultrafast antiferromagnetic spintronics combined with the topological aspects of materials.

Several years ago, another mechanism that modulates AFM orders via electric currents was proposed and demonstrated [39,40]. In antiferromagnets with a combined spaceinversion and time-reversal symmetry, electric currents induce a staggered spin density, and thus, a staggered torque termed the Néel spin-orbit torque (NSOT). The NSOT can efficiently switch the staggered magnetization, i.e., the Néel vector, thereby yielding significant progress in AFM spintronics [41][42][43][44]. Recently, it was found that a nonsymmorphic crystalline symmetry in addition to the parity-time symmetry preserves the crossing of doubly degenerate bands; this implies that Dirac fermions are controlled by the direction of the Néel vector [45]. The discovery of this close relationship between the magnetic and electronic structures opened a research field called topological AFM spintronics [46,47]. Some experiments have shown that not only electric current pulses but optical pulses can control the AFM order [48,49]; these studies performed transport measurements such as the anisotropic magnetoresistance and the planar Hall effect to observe the magnetic structure [44]. Meanwhile, the topological elec- * Deceased. tronic structure has received less attention, and therefore, the real-time dynamics of the Néel vector and the Dirac fermions induced by the NSOT are still unclear.
From a theoretical viewpoint, two approaches have been established to study the real-time dynamics of magnets. One is based on microscopic models such as the Hubbard model and Kondo lattice model [18], which allows considering electron correlation and quantum effects; however, the cluster size is severely limited and approximations are needed to deal with the symmetry-broken ordered states. The other approach is to use the Landau-Lifshitz-Gilbert (LLG) equation of the classical vectors of magnetic moments, which has been widely adopted in spintronics [50][51][52][53][54][55]. This equation enables the simulation of large-scale magnetic structures, while the electronic degree of freedom is integrated out. Therefore, a unified framework is required to reveal the relationship between the magnetic and electronic structures.
In this work, we propose a theory based on a minimal twodimensional model of AFM Dirac semimetals to investigate the real-time magnetic and electronic dynamics. We show that the Néel vector is rotated by optically-induced NSOT and other torques originating from magnetic anisotropies. This reorientation accompanies the modulation of the mass of Dirac fermions. Furthermore, the magneto-optical effect is found to be promising for the time-resolved measurement of the Néel vector.

Theoretical model
We consider a minimal tight-binding model of AFM Dirac semimetals projected onto the two-dimensional square lattice. The Hamiltonian is divided into three parts as H = H ele + H exc + H mag . The first term H ele = −2ℎ 1 cos 2 cos 2 − ℎ 2 (cos + cos ) consists of the nearest-neighbor (ℎ 1 ) and next-nearestneighbor (ℎ 2 ) hoppings, and the spin-orbit coupling ( ) of itinerant electrons with momentum k. Here, the spin and sublattice degrees of freedom of the itinerant electrons are described by the Pauli matrices σ and τ , respectively. The localized magnetic moments on each sublattice denoted by S A and S B couple to the itinerant electrons through the exchange interaction Hereafter, we define the uniform and staggered magnetizations as m = (S A + S B )/2 and n = (S A − S B )/2, respectively. The vector n is termed the Néel vector. The exchange Hamiltonian is an extension of the model for S A = −S B = n proposed by Šmejkal et al. [45]. While both the parity (P) and timereversal (T ) symmetries are broken because of the presence of the localized magnetic moments, the parity-time (PT ) symmetry is preserved when m = 0. This guarantees the Kramers degeneracy of the energy bands in the whole Brillouin zone. In addition to the PT symmetry, this model is invariant under a nonsymmorphic glide symmetry operation only when the Néel vector is directed along the 100 directions, which causes the crossings of the doubly degenerate energy bands, i.e., the gapless Dirac points, to be protected at the Brillouin zone boundary [45], as shown in Fig. 1a. The energy gap at the Dirac points depends on the value of exc and the direction of n( ). The spin-orbit coupling and the AFM ordering of the localized moments lead to the spin-momentum locking of the itinerant electrons. The expectation value of the staggered spin moments of the itinerant electrons in the lower bands is displayed in Fig. 1b, which indicates that the net staggered moment takes a negative value when n = (1, 0, 0) because of the AFM exchange interaction.
To stabilize the initial state in which the Néel vector is aligned with the [100] direction, we introduce magnetic anisotropy energy where and represent coefficients of the easy-plane and biaxial anisotropies, respectively. Here, other magnetic interactions such as the exchange interaction between the localized moments, which can be implemented in a straightforward manner, are omitted because the AFM state is stabilized by the indirect exchange mediated by the itinerant electrons at half filling.
The time-dependent electric field F ( ) parallel to the twodimensional plane is incorporated in the model via the Peierls substitution. The time evolutions of the itinerant electrons and the localized magnetic moments are described by the von Neumann and LLG equations, respectively; these equations are solved simultaneously by the fourth-order Runge-Kutta method. The parameter values and units are presented in Methods.

Application of a constant electric field
First, we discuss real-time dynamics induced by a constant electric field that is switched on at time = 0 and applied in the [110] direction, i.e., F ( ) = 0 ( ) (cos /4, sin /4, 0) with being the step function. Figure 2a shows the real-time profiles of the in-plane -component of the Néel vector and the out-of-plane -component of the uniform magnetization. Figure 2b displays the excited electron density defined by the number density of electrons that occupy single-particle energy levels above the chemical potential of the initial state. The energy bandgap shown in Fig. 2c is defined by the difference between the single-particle energies of the lowest level above and the highest level below ; this difference is equal to the direct gap when there is the particle-hole symmetry with ℎ 2 = 0. As time evolves, abruptly decreases towards zero at time = 1 , 2 , or 3 , which is indicated by arrows in Fig. 2a and is hereafter termed a reorientation time r . Meanwhile, the component of the uniform magnetization appears with values of the order of ∼ 0.03; this means that the magnetic structure is slightly canted as illustrated in Fig. 2e. As the field magnitude increases from 0 = 0.005 to 0.009, the reorientation time r becomes short. The excited electron density shown in Fig. 2b increases as time evolves. Figure 2b shows some plateaus in which the excited electron density is constant in time and the energy bandgap opens simultaneously (Fig. 2c). Once the excited electron density becomes greater than approximately 0.01 (indicated by the shade in Fig. 2b), the Néel vector starts to rotate. Note that the bandgap opens through the Néel-vector dynamics induced by the constant field, which is in contrast to the gap opening of the Floquet-Bloch bands under the high-frequency driving of electrons [56]. These phenomena induced by the constant electric field occur in a timescale of the order of ∼ 10000ℏ/ℎ 1 , which corresponds to 6.6 ps for ℎ 1 = 1 eV. Therefore, the Néel vector and the energy bandgap can be controlled by the external field in such a short timescale.
The mechanism of the reorientation of the Néel vector is understood as follows. The external electric field induces the nonequilibrium staggered spin density of the itinerant electrons through the spin-orbit coupling, and it generates torques in the [001] direction (see Fig. 2d). Then, the two sublattice magnetic moments are slightly canted towards the [001] direction, which means > 0. The easy-plane anisotropy, whose energy per sublattice moment is estimated to be ( ) 2 , induces another effective field (0, 0, −2 ) that rotates the localized moments in the anticlockwise direction as depicted in Fig. 2e. Here, we have shown the NSOT-driven reorientation of the Néel vector on the basis of the microscopic model in which the real-time dynamics of the itinerant electrons and the localized moments are explicitly considered. Figure 2f shows the polar plots of the in-plane angle of the Néel vector as a function of time multiplied by the field magnitude for different values of 0 . The Néel vector rotates towards the [010] direction for sufficiently large field magnitudes 0 0.004. For the weak magnitude 0 0.004, the trajectories of ( 0 ) are on a universal curve independent of 0 ; this implies that the timescale of the dynamics is determined by the vector potential A ∼ F rather than the electric field F . This is because the NSOT is induced by the electric current due to the diamagnetic response. The reorientation is regarded as a deviation from the weak-field universal trajectory. Even before the reorientation, the Néel vector slightly deviates from the axis because of the NSOT; once the excited electron density exceeds approximately 0.01, the NSOT-induced staggered torque shown in Fig. 2d overcomes the other torques originating from the magnetic anisotropies. These behaviors are summarized in Fig. 2g, where the reorientation time is plotted as a function of 0 for different values of the biaxial anisotropy and the Gilbert damping constant . Overall, r is inversely proportional to 0 , which means that the vector potential ∼ 0 r governs the reorientation dynamics; r becomes shorter as or decreases. It is also found that the product r is a function of 0 / for ≥ 0.01 and = 1. The reorientation time shows the step-like feature indicated by the dotted lines, which reflects a discontinuous behavior of ( 0 ) with respect to 0 as seen in the inset of Fig. 2f. Although the anisotropy energies and the Gilbert damping constant are highly material dependent, this reorientation can occur in the timescale of subpicoseconds with a sufficiently large 0 ( 0.01) irrespective of and .

Irradiation of a monocycle pulse
We show the real-time dynamics induced by a monocycle pulse of which the frequency is in the terahertz range but is off-resonant (see Methods for parameters). The electric field is parallel to the [110] direction. Figure 3a shows the inplane components of n and the out-of-plane component of m for 0 = 0.02. The Néel vector rotates towards the [010] direction during pulse irradiation, and it stops the rotation after irradiation. This reorientation accompanies the modulation of the bandgap as shown in Fig. 3b. When the bandgap is opened, the excited electron density does not increase as in the case of the constant field, which means that the pulse is off-resonant.
The dashed lines in Fig. 3a shows and for the same field amplitude but different polarization, i.e., F 0 [110], whose signs are opposite to those in the F 0 [110] case (solid lines). This polarization dependence is also seen in the case of the constant field (not shown), and it implies that the direction of the Néel vector is determined by the polarization and inverted by the mirror operation on F 0 against a plane perpendicular to n( = 0).
The pulse amplitude dependence of the Néel vector dynamics is summarized in Fig. 3c. When the amplitude is weak ( 0 < 0.01), the Néel vector cannot climb over the potential barrier of the magnetic anisotropy. However, the in-plane angle exceeds 45 • for a sufficiently large 0 (≥ 0.01) for which the excited electron density becomes greater than approximately 0.01, and the Néel vector remains in the [010] direction after irradiation for even larger 0 (≥ 0.02). Figures 3d and 3e show the in-plane angular velocity of the Néel vector and the component of the uniform magnetization, both of which are divided by 0 . For weak 0 < 0.01, the time profiles of / ( 0 ) and / 0 are on each universal curve, whereas the deviations from the curves are found for 0 ≥ 0.01 and they lead to the reorientation towards the [010] direction. These time profiles are quite similar because / is approximated by 2 as long as m n .

Optical conductivity and magneto-optical effect
Finally, we show the optical longitudinal and transverse (Hall) conductivities for different Néel vector angles, and we discuss the resulting magneto-optical Voigt effect, assuming that the system is in the collinear AFM ground states. is equivalent to under the rotation ↦ → + /2, and is equal to . The sign of the transverse conductivity is unchanged by ↦ → − or exc ↦ → − exc (not shown), which implies the conductivity is an even function of and exc . Since the present model holds the D 4h symmetry in the limit of exc n → 0, the conductivity linear in the electric field can be expanded with respect to a time-reversal-breaking field B ∼ exc n as (B) = (0) + (1) + (2) , where (1) vanishes when B ⊥ [001] and is a symmetric tensor. The off-diagonal components of (B) can be written as (B) = (B) ∝ sin 2 for n = (cos , sin , 0). This is consistent with the numerical result of , although the present values of (= 0.8) and exc (= 0.6) are beyond the perturbative regime. From these properties of and , we can identify the direction in which the Néel vector rotates by measuring the longitudinal optical conductivity and magneto-optical signal. Figures 4c  and 4d show the magneto-optical rotation angle and ellipticity calculated from and (see Methods). The rotation angle exhibits large values up to 1 • . In real materials, the values of and exc may be smaller than the values adopted here, which causes a reduction in the rotation angle down to a few millidegrees. Even so, the magneto-optical rotation can be measured using modern optical apparatus.
The transverse conductivity appears when ≠ 0 and exc ≠ 0 even though the magnetic structure is collinear AFM. Recently, the transverse response in collinear antiferromagnets is attracting considerable attention, and it has been discussed for the case in which the conductivity tensor is antisymmetric [57,58]; our calculations demonstrate the transverse optical response attributed to the symmetric tensor in collinear antiferromagnets.

DISCUSSIONS
We investigated the real-time dynamics of the magnetic and electronic structures induced by the constant and pulse electric fields, and we revealed the microscopic mechanism of the NSOT-driven reorientation of the Néel vector. Our method can be applied to any tight-binding model coupled to localized moments. The experimental verification of the present results is necessary for a better understanding of the NSOT-induced dynamics and for further progress in ultrafast AFM spintronics combined with the topological aspects of the materials. Candidate materials include orthorhombic CuMnAs [45,[59][60][61] and MnPd 2 [62], which have been studied in the field of AFM spintronics and are expected to host Dirac fermions. Note that the reorientation of the Néel vector demonstrated in this study is not directly related to what has been observed in experiments because of the discrepancy in the timescale. To make more quantitative discussions on, e.g., the timescale of the reorientation and the required magnitude of the field, we should consider the crystal structure and dimensionality of materials, obtain realistic estimates of model parameters as well as the Gilbert damping constant from the ab initio calculations and experiments, and take other dissipation processes caused by the electron-electron interaction into account; these are left for future work.
Recent intense terahertz pulse sources, whose oscillation period is 1 ps and peak amplitude is more than 1 MV cm −1 [63], are promising to drive the reorientation of the Néel vector in the picosecond timescale. In our calculations, a typical value of 0 is of order 0.01, which corresponds to 0.26 MV cm −1 (see Methods). There are many earlier studies on the electrical modulation and detection of the Néel vector [44]; their time resolution was limited to 1 ns. However, it may be possible to observe the Néel vector dynamics in real time via optical measurements because the magnetic structure is closely related to the electronic structure in the present model. One approach is the time-and angle-resolved photoemission spectroscopy; this can directly access the energy band structure, which reflects the direction of the Néel vector as the positions and energy gaps of the Dirac points or nodal lines. Another approach is to use magneto-optical effects. This system exhibits the magneto-optical rotation with angles up to 1 • . The time-resolved magneto-optical measurements can be carried out even in the subpicosecond timescale, which provides a sufficient time resolution to observe the real-time dynamics discussed in this study.

Real-time dynamics
We calculate the real-time dynamics of the itinerant electrons and the localized magnetic moments as follows [64,65]. The quantum state of the electrons is exactly described by a one-body density matrixˆk because the electronic Hamiltonian H ele contains no many-body interaction. Once the localized magnetic moments are given, the Hamiltonians in Eqs. (1) and (2) can be diagonalized for each k as (H ele + H exc )|k = k |k , where |k is obtained by the unitary transform as |k = =A,B =↑,↓ |k * , . In the initial state in which the external field is absent, the density matrix is given bỹ k = |k k k |, where k = ( − k ) is the step function with being the chemical potential chosen such that the number of electrons is ele . The time evolution of the density matrix is governed by the von Neumann equation with the initial conditionˆk ( = 0) = −1 k˜k k . The expectation value is defined by · = −1 k Tr[ ·ˆk], where represents the number of k points in the Brillouin zone. The localized moments on each sublattice denoted by S A and S B are treated as classical unit vectors, and therefore, the Hamiltonian H ( ) depends on the classical spins and the external field. The classical spins follow the LLG equation for = A, B, where denotes the Gilbert damping constant and b ≡ − H / S represents the effective field given by Here, σ A/B = (1 ± )σ/2 denotes the sublattice spin density. The time-dependent electric field F ( ) is incorporated in the model via the Peierls substitution , and denote the vector potential, elementary charge, and lattice constant, respectively. Throughout this paper, ℎ 1 , , , and ℏ are set to unity. The physical quantities of time, electric-field magnitude, and conductivity are expressed in units of ℏ/ℎ 1 = 0.66 fs, ℎ 1 /( ) = 26.3 MV cm −1 , and 2 /(ℏ ) = 6.4 × 10 3 S cm −1 , respectively, for ℎ 1 = 1 eV and = 0.38 nm. In most calculations, we use ℎ 2 = 0.08, = 0.8, exc = 0.6, = 0.1, = 0.01, and = 1. The easy-plane anisotropy coefficient determines the in-plane angular velocity of the Néel vector / . For the above parameters, the difference in the total energies of the ground states is H | n=(0,0,1) − H | n=(1,0,0) = 0.13ℎ 1 . This may be considerably larger than the anisotropy energy of real materials ( 1 meV) [60,61]. However, we confirmed that the reorientation of the Néel vector occurs within 10 ps even for smaller values: = = 0.001, = 0.08, exc = 0.06, and 0 = 0.001. Equations (5) and (6) are simultaneously solved using the fourth-order Runge-Kutta method with a time step = 0.01. We perform the numerical calculation with = 512 × 512point mesh in the Brillouin zone, and we fix the number of the electrons to ele = 2 (half-filled). The external field w )] for the monocycle pulse. We adopt c = 2000 and w = 300.

Optical conductivity and magneto-optical effect
We calculate the optical conductivity from where J k = k |J k |k denotes the matrix element of the electric current operator J k = − H (k − A)/ A| A=0 and represents a broadening factor. The dielectric tensor is given by with 0 = 8.854 pF m −1 being the electric constant. Since the optical conductivity tensor is symmetric ( = ) and anisotropic ( ≠ ) in the present model, we cannot use the well-known formula of the complex Kerr rotation angle K = /[(1− ) √ ]. According to Maxwell's equations in material media, the electric-field amplitude F 0 satisfies whereN represents a complex refractive index vector. From Eq. (11) we have the relation (ˆ2 − ) (ˆ2 − ) = 2 , assuming that the propagation vector of incident light is perpendicular to the two-dimensional -plane. The symmetric property of the conductivity and dielectric tensors leads to the eigenmodes of Eq. (11) being linearly polarized, contrary to the antisymmetric case in which the eigenmodes are circularly polarized. The boundary condition is that the tangential components of the electric field F and magnetic fieldN × F / 0 ( 0 is the speed of light) are continuous at the surface, from which we obtain two eigenmodes and complex refractive indices. The incident linearly polarized light is decomposed into a linear combination of the eigenmodes. The amplitude vector of the reflected light is proportional to ( pp , sp , 0), where pp and sp are the amplitude reflection coefficients of the incident p-polarized light to the reflected p-and s-polarized light. The magneto-optical rotation angle and ellipticity are determined by the fit of the electric-field trajectory to an ellipse rotated around the axis by angle , and the relation = sgn[arg( sp / pp )] × | sp |/| pp |, respectively.

ACKNOWLEDGMENTS
One of the authors, Sumio Ishihara, passed away on 7 November 2020 during the preparation of the manuscript. We thank H. Matsueda and Y. Masaki for their valuable comments and critical reading of the manuscript. This work was supported by JSPS KAKENHI Grant Nos. JP19K23419, JP20K14394, JP17H02916, JP18H05208, and JP20H00121. The numerical calculations were performed using the facilities of the Supercomputer Center, the Institute for Solid State Physics, the University of Tokyo.

AUTHOR CONTRIBUTIONS
A.O. conceived the idea, performed the calculations, and wrote the manuscript. A.O. and S.I. discussed and interpreted the results.