Coherent two-dimensional THz magnetic resonance spectroscopies for molecular magnets: Analysis of Dzyaloshinskii–Moriya interaction

Electron paramagnetic resonance (EPR) and nuclear magnetic resonance (NMR) have a long history as typical examples of two-dimensional (2D) spectroscopy techniques that measure the varying time intervals of magnetic pulse trains applied to electron or nuclear spin systems. [1, 2] Although these techniques are powerful means for structural analysis of organic and inorganic materials, it is difficult to apply them to the investigation of spin dynamics because the time resolution of the magnetic pulses involved is limited to the order of microseconds.

Recently, coherent terahertz (THz) magnetic resonance (MR) spectroscopy was developed using magnetic pulses in the subpicosecond range generated by the magnetic-field component of THz light. [3, 4, 5] Just as ultrafast laser spectroscopy has made it possible to study the electronic excitation-state dynamics and intra- and intermolecular vibrational motions of complex molecular systems, [6] THz-MR spectroscopy has opened up the possibility of investigating strongly correlated spin dynamics in molecular magnets; this depends on the complex spin–spin interactions and the configurations of spins in the molecular environment, which cause dephasing and relaxation of the spin system.

The quasi-ferromagnetic (FM) and antiferromagnetic (AFM) precession modes in $\rm{YFeO_{3}}$ (YFO) that arise from the antisymmetric spin–orbit coupling [7] have been observed as free induction decay (FID) signals. [8, 9] Spin-wave excitation in AFM $\rm{NiO}$ has been detected with the aid of the Faraday effect. [10, 11] The skyrmion, which is a quasiparticle composed of vortex-like spin orientations, has been detected in several FM and AFM materials with thin-film structures. [12, 13, 14] THz electron spin resonance (ESR) and Hall conductivity spectroscopy have been employed to investigate the topological Hall effect and phase transitions. [15, 16, 3, 17]

Although THz-MR and THz-ESR spectroscopic approaches based on linear response theory have been successful for the classification of complex spin states in condensed phases, [18, 19, 20] it is unclear whether the spin states that are investigated are quantum-mechanically entangled, and whether the width of the peaks that are measured arise from the relaxation or dephasing processes, because these peaks are usually broad and overlap. Thus, an extension to coherent 2D THz-MR spectroscopy was developed. The observed 2D THz-MR spectrum for a YFO crystal made it possible to illustrate characteristic nonlinear spin responses, such as double-quantum (2Q) coherence and second-harmonic generation (SHG). [21]

As already illustrated in 2D optical spectroscopy, [22, 23] 2D THz-MR spectroscopy is not only useful for identifying complex spin interactions but also for monitoring the complex quantum spin dynamics under the influence of relaxation and dephasing arising from the environment at the femtosecond scale. Theoretical input regarding the complex profiles of spin–spin interactions is important for analyzing these 2D spectra under ultrafast nonlinear processes. In particular, this includes those in molecular magnets, which play a central role in spintronics and next-generation information technologies. [24, 25, 26]

In this paper, we provide a comprehensive theoretical framework for both 1D and 2D THz-MR measurements based on the response-function theory developed for nonlinear optical spectroscopy techniques. To illustrate our approach more closely, we employ a chiral spin chain described as a Heisenberg model with exchange coupling and Dzyaloshinskii–Moriya interaction (DMI) arising from the antisymmetric spin–orbit coupling. [27, 28, 29] Such anisotropy of the spin system leads to a series of unique magnetic and optical properties. [30] Unlike with conventional FM materials, in addition to the magnetic anisotropy, a non-centrosymmetric structure arises as the result of the chiral ordering of spins. As a result, a series of new phenomena arising from the nonlinear magnetic response, such as magneto-chiral dichroism and magnetization-induced SHG, were observed. [31, 32, 33, 34]

To investigate the properties of such materials as devices, it is necessary to investigate ultrafast spin dynamics under time-dependent external fields, in which quantum coherence and entanglement—not only among spins but also between spins and the environment—play a significant role. [2, 35] Thus, we include a harmonic heat bath of finite temperature in the spin system. The number of degrees of freedom of the heat bath is then reduced to obtain a time-irreversible equation of motion describing the effects of thermal fluctuations and dissipation. [36, 37] Because the motion of the spins is much faster than the thermal noise arising from the heat bath, the heat bath must be treated in a non-Markovian manner. Thus, the reduced equations of motion derived using Markovian approximations and other assumptions—such as the Bloch, Lindblad, and Redfield equations—are not suitable for the description of such ultrafast spin dynamics. We then employ the numerically “exact” hierarchical equations of motion (HEOM) approach, which can be used to treat non-Markovian and non-perturbative system–bath interactions at finite temperatures. [38, 39, 40, 41]

To demonstrate the applicability of the present theory, 1D and 2D THz-MR spectra were calculated for the chiral spin model with different DMI strengths describing the pitch and direction of chirality (clockwise or counterclockwise). We show that in 1D spectra, only the absolute strength of the DMI can be evaluated, whereas in 2D spectra, the sign of the DMI, which determines the direction of the chirality, can also be determined. While neutron-scatting techniques have been used to determine the structures of chiral materials, the present results indicate the possibility of determining the DMI through spin-dynamic processes. This finding should be valuable for the design of spintronic materials.

The rest of this paper is organized as follows. In Sec. II, we formulate linear (1D) and 2D THz-MR spectra based on the response-function theory. In Sec. III, we introduce the Heisenberg model with the DMI coupled to the harmonic heat bath. The HEOM are then presented. Numerical results and some discussion are presented in Sec. IV. Finally, Sec. V is devoted to our conclusions.

We consider a spin system coupled to a bath system driven by an external magnetic field $B(s)$. The total Hamiltonian is expressed as

where $\hat{H}_{tot}$ is the Hamiltonian of a composite system and $\hat{M}$ is the polarization operator for magnetic fields.

The observable in a magnetic measurement at time $t$ is expressed as $M(t)\equiv\langle\hat{M}(t)\rangle-\langle\hat{M}\rangle$, where $\hat{M}(t)$ is the Heisenberg representation of $\hat{M}$ for $\hat{H}^{\prime}(s)$. The thermal average for any operator $\hat{A}$ is defined as $\langle\hat{A}\rangle\equiv{\rm{tr}}\{\hat{A}\hat{\rho}_{tot}^{eq}\}$, with the equilibrium density operator expressed as $\hat{\rho}_{tot}^{eq}$. We can also express this as ${M}(t)={\rm{tr}}\{\hat{M}\mathcal{G}^{\prime}(t)\hat{\rho}_{tot}^{eq}\}-{\rm{tr}}\{\hat{M}\hat{\rho}_{tot}^{eq}\}$, where the Liouville operator is defined as

Here, the arrows indicate time-ordered exponentials.

Experiments such as 2D NMR and 2D EPS measurements [1, 2] are conducted using multiple magnetic pulses with finite time widths. In such experiments, the excitation by the external field is non-perturbative, and the desired spin dynamics are investigated by designing the profiles of pulse trains, as in the cases of spin-echo and correlation spectroscopy measurements. Theoretically, these signals are obtained by integrating the equations of motion for a spin system, such as the Bloch, [42] Redfield, [43] or stochastic Liouville equations, [44] or the HEOM, [40, 45, 46] under a sequence of magnetic pulses.

As in the case of coherent optical laser spectroscopies, the excitations of coherent THz-MR spectroscopy are assumed to be impulsive, which allows us to measure the signals in different orders of the field-system interactions separately. Thus, we can employ a response-function theory developed for ultrafast nonlinear laser spectroscopy. [6] Up to the third order, the signal is then expressed as

where the linear, second-order, and third-order response functions are defined as [40]

and

Here, $\mathcal{G}(s)$ is the Liouvillian without magnetic interaction, which is defined from $\mathcal{G}^{\prime}(s)$ with $B(s)=0$, and we have introduced the superoperator $\hat{M}^{\times}{{\hat{\rho}}}\equiv[\hat{M},{\hat{\rho}}]$. The right-hand side (rhs) of the second line in each of the above equations allows us to employ the equations of motion to calculate the response functions and give us an intuitive picture of the higher-order optical processes. [40, 41] For example, the rhs of the second line of Eq. (4) can be read from right to left as follows. The total system is initially in the equilibrium state $\hat{\rho}_{tot}^{eq}$. The initial state is then modified by the first magnetic pulses via the dipole operator as $(\hat{M}^{\times}\hat{\rho}_{tot}^{eq})=[\hat{M},\hat{\rho}_{tot}^{eq}]$ at $t=0$ and is propagated for time $t=s$ by $\mathcal{G}(s)$. The expectation value is then obtained by calculating the trace of $\hat{M}$.

Actual 2D experiments have been conducted using a pair of magnetic pulses $a$ and $b$ with inter-pulse delay $\tau$. [21, 47] Under the impulsive approximation, the magnetic field is expressed as

where $B_{a}$ and $B_{b}$ are the magnetic field strengths. The nonlinear element of the signal $M_{NL}(t,\tau)$ at $t+\tau$ is then evaluated as

where $M_{ab}(t+\tau)$, $M_{a}(t)$, and $M_{b}(t)$ are the total signal and the linear elements for pulses $a$ and $b$, respectively.

From Eqs. (3) and (7), the nonlinear element is evaluated as

Thus, the characteristic features of THz-MR spectroscopy are described by the nonlinear susceptibilities. Because each term on the rhs has different proportionality with respect to $B_{a}$ and $B_{b}$, they can be evaluated separately by changing their respective field strengths.

In nonlinear optical spectroscopies, $\chi_{2}$ is used to analyze 2D THz-Raman [48, 49] and 2D THz-IR-Raman signals [50, 51] that involve a 2Q transition process, while $\chi_{3}(s_{3},s_{2},s_{1})$ is used to analyze nonlinear 2D THz spectroscopies [52, 53] including 2D THz rotational spectroscopy. [54, 55] As illustrated in nonlinear optical spectroscopies, the $s_{1}$ and $s_{3}$ periods describe the time evolutions of coherent states, while $s_{2}$ describes the evolution of population states. [22, 23] While $\chi_{3}(\tau,0,t)$ is used for 2D spectroscopies with waiting time $s_{2}=0$, $\chi_{3}(0,\tau,t)$ is used for transient absorption spectra in which the time evolution of the excited-state population is measured.

Using $\chi_{1}(s)$ in Eq. (4), the linear absorption spectrum is defined as

where Im denotes the imaginary part. Following the experimental setup, we consider two kinds of 2D spectrum for different time configurations of $\tau$ and $t$ in Eq. (6), expressed as: (1) $\chi_{3}(\tau,0,t)$ and (2) $\chi_{3}(0,\tau,t)$. They are then expressed in the Fourier translation form as

where $\omega_{\tau}$ and $\omega_{t}$ represent the excitation and detection frequencies, respectively. [21]

To describe the novel phenomena of chiral magnets, we consider a linear chain system consisting of $L$ spins. The system Hamiltonian is defined as

where $\hat{\sigma}_{n}^{\alpha}$ $(\alpha=x,y,z)$ denotes the Pauli operator at the $n$-th site in the $\alpha$-th direction. Here, the first term with $J>0$ represents the AFM Heisenberg exchange coupling, and $\Delta$ is the anisotropic parameter. The second term in Eq. (13) is the DMI that is perpendicular to the $xz$ plane, with coupling strength $D_{y}$. The sign of $D_{y}$ determines the direction (i.e., clockwise or anticlockwise) or handedness (i.e., right- or left-handed circular configuration) of chirality, which is not easily determined by experimental measurements. [3, 29] The magnetization operator is defined as

A very important aspect of investigating ultrafast dynamics in magnetic materials is the inclusion of a quantum-mechanically consistent relaxation and dephasing mechanism. This can be achieved by including a harmonic heat bath in the spin system. The Hamiltonian of this heat bath is defined as

where $\hat{p}_{j}$, $\hat{x}_{j}$, $\omega_{j}$, and $m_{j}$ represent the momentum, position, frequency, and mass of the $j$-th oscillator, respectively. The system–bath interaction is defined by

where $\hat{V}$ represents the spin part of the system–bath interaction function, defined as

and $g_{j}$ is the coupling strength with the $j$-th oscillator. The total Hamiltonian is then given by

In open quantum dynamics theory, the time-irreversible process can be described using a reduced set of equations of motion. After reducing the number of degrees of freedom of the heat bath, the noise effects are characterized by the correlation function,

where $J(\omega)$ is the spectral density function (SDF), $\beta=1/k_{B}T$ is the inverse temperature, and $k_{B}$ is the Boltzmann constant. For the heat bath, we assume the Drude SDF, which is expressed as,

where $\zeta$ is the coupling strength and $\gamma$ is the inverse correlation time.

In non-perturbative and non-Markovian conditions, the time evolution can be described by the HEOM approach. [38, 39] We rewrite Eq. (19) in a linear-summation form,

where $c_{k}$ and $\nu_{k}$ are complex-valued coefficients. The HEOM can then be expressed as [40, 41]

where $\vec{n}=\{n_{1},n_{2},\cdots,n_{K}\}$ denotes the index vector, and $n_{k}$ are non-negative integers. Among all the $\hat{\rho}_{[\vec{n}]}(t)$, the $0$-th, with $\vec{0}=\{0,0,\cdots,0\}$, corresponds to the density operator of the reduced system, while all the others are introduced for ancillary purposes.

In the HEOM approach, the response functions Eqs. (4)–(6) are evaluated as the time evolution of the system under external excitation. The density matrix is replaced by the HEOM elements, and the Liouvillian $\mathcal{G}(t)$ is replaced using Eq. (22).

For example, we can evaluate $\chi_{3}(\tau,0,t)$ using the expression of the second line in Eq. (6) as follows. [40, 41] We first run the HEOM program for a sufficiently long period from a temporally initial condition at $t=-t_{i}$ (such as the factorized initial condition, $\hat{\rho}_{[\vec{0}]}(-t_{i})=\exp[-\beta\hat{H}_{S}]$ with all the other hierarchy elements set to zero) to time $t=0$ to reach a true thermal equilibrium, denoted by $\hat{\rho}_{[\vec{n}]}^{eq}(0)$. All the hierarchy elements have to be used to define a correlated initial condition. The system is next excited by the first magnetic interaction $\hat{M}$ at $t=0$ as $\hat{\rho}_{[\vec{n}]}^{\prime}(0)=[\hat{M},\hat{\rho}_{[\vec{n}]}^{eq}(0)]$. The perturbed elements $\hat{\rho}_{[\vec{n}]}^{\prime}$ then evolve in time by numerically integrating Eq. (22) up to $t$. At $t$, the system is excited by the second and third magnetic interactions as $\hat{\rho}_{[\vec{n}]}^{\prime\prime}(t)=[\hat{M},[\hat{M},\hat{\rho}_{[\vec{n}]}^{\prime}(t)]]$. Then, after $\hat{\rho}_{[\vec{n}]}^{\prime\prime}$ evolves in time with the initial condition $\hat{\rho}_{[\vec{n}]}^{\prime\prime}(t)$ to $t+\tau$, the response function is calculated from the expectation value of the magnetic moment as $\chi_{3}(\tau,0,t)=tr_{A}\{\hat{M}\hat{\rho}_{[\vec{0}]}^{\prime\prime}(t+\tau)\}$.

Notice that to take the system–bath coherence (or bath entanglement [41]) into account during the external perturbation, it is important to operate $\hat{M}$ to all of the hierarchy elements $\hat{\rho}_{[\vec{n}]}^{\prime\prime}(t)$. Although we only use $\hat{\rho}_{[\vec{0}]}(t)$ to calculate an expectation value, the other elements are essential to obtain an echo signal for non-Markovian noise in 2D spectroscopy. [40, 41]

In the following, we set the exchange-coupling strength as the base unit $J=1$ and calculate 1D (or linear absorption) and 2D spectra for $\Delta=0.1$ as a model for typical AFM metal oxides. [56, 57] The number of spin sites was $L=6$, with periodic boundary conditions $\hat{\sigma}_{n}^{\alpha}=\hat{\sigma}_{n+L}^{\alpha}$. Here, we consider cases without the DMI ($D_{y}=0$) and with the weak DMI ($|D_{y}|\leq 0.05$), which is appropriate for a small model system. Although the HEOM used in this study can investigate cases in which the system interacts strongly with the heat bath, here, we focus on the characterization of the chiral spin system and keep the coupling weak. Thus, the bath parameters were chosen as $\zeta=0.05$, $\gamma=1.0$, and $\beta\hbar=10$. The hierarchy parameters were chosen as $K=50$ and $\sum_{k}^{K}n_{k}\leq 10$, and the time step for numerical integration was set to $\delta t=0.01$.

A recently developed 2D THz-MR spectroscopic technique has created new possibilities for measuring complex molecular magnetic systems. In the present work, we have illustrated the key features of this technique based on nonlinear response-function theory and have described a method for simulating 2D THz-MR spectra through the use of the HEOM formalism. Using simulated 1D and 2D THz-MR spectra for a chiral magnetic material, we demonstrated that the 2D technique allows us to evaluate the pitch and direction of chirality determined from the strength and sign of the DMI, while only the absolute amplitude of the DMI can be determined from 1D measurements.

The reason the sign of the DMI was detected by 2D spectroscopy in this study is that there are eigenstates in which the phase of the wavefunction changes following the sign of the DMI, which is categorized as (III). Although there have been studies of eigenstates of spin systems with the DMI, [58] the existence of eigenstates that change phase with the sign of the DMI, as found in this study, has not been explored. It is important to know the causes of these eigenstates because they may lead to the appearance of novel phenomena.

To make a direct comparison between the results of our simulations and those obtained experimentally, however, we must increase the number of spins in accordance with those available in experimental systems. In investigating spin dynamics in the condensed phase, it is also important to consider the non-perturbative and non-Markovian system–bath interactions. Nevertheless, we believe that the present results elucidate the key features of 2D THz-MR spectroscopic methods with regard to probing the fundamental nature of a magnetic spin system.

For further investigations to monitor the ultrafast dynamical aspects of the spin system, such as the dynamics of spin waves, it is necessary to conduct a variety of advanced nonlinear spectroscopic approaches, such as pump-probe and transient absorption measurements to foster the development of this spectroscopic method. We leave such extensions to future studies, depending on progress in experimental and simulation techniques.