Dynamical and static spin structure factors of Heisenberg antiferromagnet on honeycomb lattice in the presence of Dzyaloshinskii-Moriya interaction

Quantum magnetism on geometrically two-dimensional (D) frustrated spin systems with $S=1/2$ have lately received massive attetions, due to their potential for realizing the quantum spin liquid, a magnetically disordered state which respects all the symmetries of the systems, even at absolute zero temperatureokuma ; zvygain ; fihlo ; ito ; balent . The spin model, recently attracted many interests, is the Heisenberg model with first and second antiferromagnetic exchange interaction in honeycomb lattice. In sufficiently low dimensional spin systems, the quantum mechanical zero point motion can forbid long range magnetic order and produce a quantum spin liquid state, a correlated state that breaks no symmetry and possesses topological properties, possibly sustaining fractionalized excitationsanderson1 ; fazekas ; liang ; sachdev1 ; sandvik . Although the triangular lattice was first theoretically proposed by Andersonanderson1 as an ideal benchmark to search for the quantum spin liquid, it was soon found that the $S=1/2$ antiferromagnetic Heisenberg model on a triangular lattice is magnetically ordered with a 120$\textdegree$ arrangement of the spins. Despite the intense activity, only a small number of triangular materials have been identified as possible candidates for quantum spin liquid behavior such as th layered organic materials. Hence, there is a need to find evidence for quantum spin liquid behavior in more compounds. Honeycomb lattice materials have attracted lots of attetion in recent years due to their interesting and poorly understood magnetic properties. Inorganic materials such as Na₂ Co₂TeO₆lefan , BaCo₂(AsO₂)₂ martin , are examples of honeycomb lattice spin half antiferromagnets. It is important then to understand theoretically the magnetic properties of interacting localized moments on the frustrated honeycomb lattice as has been previously done on triangular lattices. Althogh the numerical evidence for a quantum spin liquid in the half-filled Hubbard model on the honeycomb lattice has been questionedsorella , exact diagonalization studies on the $J_{1}-J_{2}$ Heisenberg model with $S=1/2$ have found evidence for short range spin gapped phases for $J_{2}=0.3-0.35$ suggesting the presence of a Resonance Valence Bond (RVB) state. Slow neutrons scatter from solids via magnetic dipole interaction in which the magnetic moment of the neutron interacts with the spin magnetic moment of electrons in the solid. The inelastic cross-section for magnetic neutron scattering from a magnetic system can be expressed in terms of spin density correlation functions. Specially, field induced effects on the dynamical spin correlation function in low dimensional quantum spin models have been attracting much interest from theoretical and experimental point of view in recent years affleck ; uimin ; dimitriev . AF Heisenberg chain in the presence of axial anisotropy at finite non zero values for magnetic field is a solvable model. Its ground state properties have been investigated within Bethe-Ansatzbethe ; faddeev . Static thermodynamic quantities, eg. the specific heat, the magnetic susceptibility and magnetization have been investigated by several methods including thermodynamic Bethe Ansatz, Quantum Monte Carlo and Density Matrix Renormalization Groupklumper ; johnston ; klumper1 .

Spin-orbit coupling induces both symmetric and antisymmetric anisotropic properties or Dzyaloshinskii-Moriya spin anisotropymoriya in the exchange coupling between nearest neighbour spins. A Dzyaloshiskii-Moriya (DM) interaction with a DM vector perpendicular to the layer produces an easy-plane spin anisotropic Hamiltonian for honeycomb lattice. DM interaction breaks $SU(2)$ symmetry of the model and reduces it to $U(1)$ symmetry around the axis perpendicular to the honeycomb plane. This DM interaction is believed to orient the spins in the 2D layer. For applied large magnetic fields ($B$) along perpendicular to the plane, the Zeeman term overcomes the antiferromagnetic spin coupling and the ground state is a field induced ferromagnetic state with gapped magnon excitations. Decreasing the magnetic field at the zero temperature, the magnon gap vanishes at the critical field ($B_{c}$) and a spiral transverse magnetic ordering develops. At finite temperature ($T$) and in the absence of exchange anisotropy, the frustration leads to an incommensurate spiral structure for transverse spin component below the critical line in the $B-T$ plane.

Hard core bosonic representation has been introduced to transform a spin Hamiltonian to bosonic one whereby the excitation spectrum is obtainedauerbach ; mahan . This mapping between the bosonic gas and original spin model is valid provided to add the hard core repulsion between particles in the bosonic Hamiltonianmahan . An anisotropic exchange interaction due to DM interaction adds a hopping term for bosonic particles to the main part of bosonic model Hamiltonian. Such term leads to different universal behavior in addition to change of critical points and thermodynamic properties.

In this paper, we intend to find the effects of antisymmetric exchange interaction (DM interaction) and magnetic field on dynamical spin structure factor of localized electrons on honeycomb lattice at finite temperature in the field induced gapped spin-polarized phase. Specially, we study the effects of Dzyaloshinskii-Moriya interaction strength and next nearest neighbor exchange constant on the spin structure factors of Heisenberg model on honeycomb structure for magnetic fields above critical field $B_{c}$. We have addressed the frequency behaviors of transverse and longitudinal structure factors of anisotropic antiferromagnetic Heisenberg model on honeycomb lattice for different values of magnetic field above threshold field. We have implemented the hard core boson transformation for spin operators which gives the excitation spectrum in terms of many body calculations for bosonic gasgorkov . We have used Brueckner approach gorkov to find the bosonic self-energy to get the spinon dispersion relation. In order to calculate structure factors of the spin model Hamiltonian including antisymmetric anisotropic Dzyaloshinskii-Moriya term, we have used one particle excitation spectrum of hard core bosonic model Hamiltonian. We obtain both transverse and longitudinal spin structure factors from one and two particle bosonic Green’s function. Also the effects of next nearest neighbor exchange constant and DM interaction strength on temperature dependence of transverse static spin susceptibility have been investigated in details.

The general form of anisotropic spin Hamiltonian on the monolayer honeycomb lattice which contains antisymmetric exchange interaction between spins in the presence of a magnetic field is written in the following form

The lattice structure of honeycomb lattice has been shown in Fig.(1). The symbol $\langle ij\rangle$ and $[ij]$ in Eq.(1) implies the nearest neighbor and next nearest neighbor sites in a honeycomb lattice, respectively. $J$ is the nearest neighbor coupling constant between spins, while $J^{\prime}$ denotes the next nearest neighbor coupling constant between spins. $z$ is thae axis perpendicular to the honeycomb plane. The antisymmetric spin exchange originating from the relativistic spin-orbit interaction is also known as Dzyaloshinskii-Moriya (DM) interaction. The third term in Eq.(1) with ${\bf D}=(0,0,D)$ describes the DM interaction between nearest neighbor sites. Also $g\simeq 2.2$ is the gyromagnetic constant and $\mu_{B}$ denotes the Bohr magneton. The gyromagnetic constant $g$ value is independent of the type of lattice and gets the same value for honeycomb lattices such as graphene. This constant is related to the structure of localized electron in the lattice. $B$ is the magnetic field strength. The magnetic field is applied as perpendicular to honeycomb plane. Anisotropy due to DM interaction and Zeeman term violate SU(2) symmetry of the isotropic Heisenberg model Hamiltonian.

In order to obtain the bosonic representation of the model hamiltonian two different bosonic operators are required. Therefore spin operators are transformed to bosonic ones as

$l$ denotes the unit cell index and $a,b$ label two sublattices. $a^{{\dagger}}_{l}(b^{{\dagger}}_{l})$ creates a boson in unit cell with index $l$ on sublattice $a(b)$. Exploiting above transformation, we have the following one particle bosonic hamiltonian

Based on Fig.(1), the symbol $\Delta=\vec{\Delta}_{1},\vec{\Delta}_{2},-\vec{\Delta}_{1},-\vec{\Delta}_{2}$ implies the characters of neighbor unit cells including nearest neighbor atomic sites. The translational vectors $\vec{\Delta}_{1},\vec{\Delta}_{2}$ connecting neighbor unit cells are given by

Also index $\Delta^{\prime}=\vec{\Delta}^{\prime}_{1},\vec{\Delta}^{\prime}_{2}$ implies the characters of neighbor unit cells including next nearest neighbor atomic sites. The translational vector $\Delta^{\prime}_{1},\vec{\Delta}^{\prime}_{2}$ are given by

where the length of until cell vectors is considered to be one. In terms of Fourier space transformation of Bosonic operators, the bilinear part of model Hamiltonian is given by

The coefficients in the above equation are

The wave vectors ${\bf k}$ are considered in the first Brillouin zone of the honeycomb lattice. Also the quartic part of the model Hamiltonian is obtained as follows

In order to reproduce the SU(2) spin algebra, bosonic particles must also obey the local hard-core constraint , i.e only one boson can occupy a single site of lattice. In terms of Fourier transformation of bosonic operators, we can write this part of the Hamiltonian as

where $\mathcal{U}$ implies the strength of the hard core repulsion. In fact $\mathcal{U}$ is the local intrasite integration between hard core bosons. The effect of hard core repulsion part ($\mathcal{U}\rightarrow\infty$) of the interacting Hamiltonian in Eq.(4) is dominant compared with the quartic term in Eq.(8). Under taking limit $\mathcal{U}\rightarrow\infty$ the excitation spectrum of bosonic gas describe the excitation spectrum of original spin model Hamiltonian.

Using a unitary transformation as

the bilinear part of the Hamiltonian is diagonalized as

The Bogoliubov coefficients $u,v$ are given by

According to bilinear part of bosonic model Hamiltonian in Eq.(6), Fourier transformation of the noninteracting Green’s function matrix elements are written in the following form

where $\omega_{n}=2n\pi k_{B}T$ denotes the bosonic Matsubara’s frequency. Also $T$ denotes the temperature of the system. In Eq.(13), symbol $\mathcal{T}$ means time-ordering operator. Since the hard core bosonic repulsion is on-site interaction between bosons, Wick’s theoremmahan concludes that only the diagonal elements of bosonic self-energy matrix gets non zero value. Thus we can write down a Dyson’s equation for the interacting Green’s function matrix (${\bf G}({\bf k},i\omega_{n})$) as

where $\Sigma^{Ret}_{aa}({\bf k},\omega)=\Sigma^{Ret}_{bb}({\bf k},\omega)$ is normal retarded self-energies diagonal element due to hard core repulsion between bosons. After using Dyson’s series mahan for the interacting Matsubara Green’s function matrix in Eq.(19) , the low energy limit of single particle retarded Green’s function matrix elements are

The renormalized excitation spectrum and renormalized single particle weight are given by

Since the Hamiltonian $\mathcal{H}_{int}$ in Eq.(4) is short ranged and $\mathcal{U}$ is large, the Brueckner’s approach (ladder diagram summation) rezania2008 ; fetter ; gorkov can be employed for calculating bosonic self-energies in the low density limit of bosonic gas and for low temperature. Firstly, the scattering amplitude (t-matrix) $\Gamma(p_{1},p_{2};p_{3},p_{4})$ of hard core bosons is introduced where $p_{i}\equiv({\bf p},ip_{m})_{i}$. ${\bf p}$ is momentum vector and $p_{m}\equiv 2m\pi k_{B}T$ implies the Matsubara’s frequency of the $i$th bosonic particle. The basic approximation made in the derivation of $\Gamma(K\equiv p_{1}+p_{2})$ is that we have considered the diagonal elements of bosonic Green’s function for making $\Gamma$. It is due to the strongly intrasite interaction between bosons. According to the Feynman’s rules in momentum space at finite temperature ($T$) and after taking limit $\mathcal{U}\longrightarrow\infty$, the scattering amplitude can be written as the following form

where $N$ is the number of unit cells and wave vector ${\bf Q}$ belongs to the first Brillouin zone of honeycomb lattice. $k_{B}$ is the Boltzmann constant. We can perform the summation over Matsubara frequencies $Q_{m}=2m\pi k_{B}T$ in the Eq.(22) according to Feynman rulesmahan and the final result for scattering amplitude is obtained by

where $n_{B}(x)=\frac{1}{e^{x/k_{B}T}-1}$ describes the Bose-Einstien distribution function. The hard core self-energy is obtained by using the scattering amplitude obtained in Eq.(22).

The hard core self-energy is obtained by integration over internal Matsubara’s frequency ($p_{m}$)

It should be noted that index $m$ is the integer number. We can simply obtain the retarded self-energy by analytic continuation ($i\omega_{n}\longrightarrow\omega+i0^{+}$) of Eq.(25). The contribution of $\mathcal{H}_{quartic}$ on the final results is very small beecause it is composed of quartic terms in the bosonic operators. It is therefore treated in mean-field approximation. This is equivalent to take only one-loop diagrams in to account (first order in $J$). On the mean field level we have $O_{1}O_{2}=\langle O_{1}\rangle O_{2}+\langle O_{2}\rangle O_{1}-\langle O_{1}\rangle\langle O_{2}\rangle$ where each $O_{1}$ and $O_{2}$ is a pair of boson operators. We can write for each pair of operators,

Thus, the effect of $H_{quartic}$ is to renormalize the coefficients of bilinear part of the Hamiltonian according to the following relations

The renormalized coefficients (Eq.(27)) will be considered to calculate the self-energy which are independent of energy (nonretarded in time representation).

The correlation function between spin components of localized electrons at different times can be expressed in terms of one and two particle bosonic Green’s functions. The frequency Fourier transformation of this correlation function produces the dynamical spin susceptibility. The frequency position of peaks in dynamical spin susceptibility are associated with collective excitation localized electrons described by Heisenberg model Hamiltonian. These excitations are related to the spin excitation spectrum of localized electrons on Honeycomb structure. In the view point of experimental interpretation, the dynamical spin susceptibility of the localized electrons of the system is proportional to inelastic cross-section for magnetic neutron scattering from a magnetic system that can be expressed in terms of spin density correlation functions of the system. In other words the differential inelastic cross section $d^{2}\sigma/d\Omega d\omega$ is proportional to imaginary part of spin susceptibilities. $\omega$ denotes the energy loss of neutron beam which is defined as the difference between incident and scattered neutron energies. The solid angle $\Omega$ implies the orientation of wave vector of scattered neutrons from the localized electrons of the sample. We can assume the wave vector of incident neutrons is along $z$ direction. The solid angle $\Omega$ depends on the polar angle between wave vector of scattered neutrons and the wave vector of the incident neutrons. A note is in order here. The solid angle $\Omega$ is different from the excitation spectrum of hard core bosonic particles , i.e. $\Omega_{\alpha}({\bf k})$, $\Omega_{\beta}({\bf k})$, in Eq.(21). The frequency position of peaks in $d^{2}\sigma/d\Omega d\omega$ determines the spin excitation spectrum of the magnetic systemdoniach . In order to study the general spin excitation spectrum of the localized electron of the systems, both transverse and longitudial dynamical spin-spin correlation functions have been calculated.

According to the linear-response theory, the components of spin susceptibility tensor are written by

that $\eta,\gamma=x,y,z$ impliy the spatial directions and $\omega_{n}=2n\pi k_{B}T$ is the bosonic frequency. ${\bf q}$ is the wave vector belonging to the first Brillouin one of honeycomb lattice. Since each unit cell in honeycomb lattice includes two different sublattices, the Fourier transformation of spin operators in Eq.(28) is exprssed in terms of both bosonic operators $a_{\bf q}$ and $b_{\bf q}$ as

Here, we define transverse dynamical spin suceptibility, i.e. $\chi_{+-}({\bf q},i\omega_{n})$, in terms of components of spin susceptibility tensor, $\chi_{\eta\gamma}({\bf q},i\omega_{n})$,

so that spin ladder operators are given as $S_{+}=S_{x}+iS_{y}$ and $S_{-}=S_{x}-iS_{y}$. Using the definitions of Fourier transformation of spin operators in Eq.(29), we can calculate the Matsubara representation of transverse dynamical spin suceptibility in terms of one particle bosonic Green’s function as

where the matrix elements of interacting Green’s function ($G$) have been expressed in Eq.(20). Also the Matsubara’ form of longitudinal dynamical spin susceptibility is written based on two particle Green’s function

After using Wick’s theoremmahan and performing Matsubara frequency summation rules, the following expression has been obtained for longitudinal dynamical spin susceptibility

The dynamical spin structure factor for both longitudinal and transverse spin directions are obtained based on retarded presentation of susceptibilities as

so that $\chi_{zz}^{Ret}$ and $\chi_{+-}^{Ret}$ are retarded dynamical spin structure factors for longitudinal and transverse components of spins, respectively. It is proportional to the contribution of localized spins in the neutron differential cross-section. For each ${\bf q}$ the dynamical structure factor has peaks at certain energies which represent collective excitations for bosonic gas which correspond to the spin excitation spectrum of the original model.

Static transverse spin structure factor ($s({\bf q})$) which is a measure of magnetic long range ordering for spin components along the plane, i.e. transverse direction, can be related to imaginary part of retarded dynamical spin susceptibility using following relation.

Substituting Eq.(20) into Eq.(31) and using Matsubara frequency summation rules Eq.(35), the final result for transverse static structure factor is given by

In the next section, the numerical results of dynamical spin structure and static spin structure have been presented for various magnetic field and DM interaction strength.

The numerical reuslts of spin structure factors of the spin 1/2 Heisenberg model on honeycomb lattice in the field induced spin-polarized phase have been presented in this section. The frequency behaviors of dynamical spin susceptibilities and temperature dependence of static spin structure factors are studied. The effects of both magnetic field and Dzyaloshinskii-Moriya term strength on the behaviors of structure factors are investigated. In the limit of $B/J\longrightarrow\infty$, the ground state of the original spin model Hamiltonian is a field induced spin-polarized state and a finite energy gap exists to the lowest excited state. The decrease in magnetic field lowers the energy gap which eventually vanishes at the critical magnetic field ($B_{c}$). The effects of hard core interaction on the excitation spectrum have been obtained by the Green’s function approach in the context of Brueckner’s formalism above threshold field $B_{c}$ where the density of bosonic gas is small.

The single particle excitation should be found from a self-consistent solution of Eqs.(21,23,25,27) with the substitutions $u_{\bf k}\longrightarrow\sqrt{Z_{\bf k}}U_{\bf k}$, $v_{\bf k}\longrightarrow\sqrt{Z_{\bf k}}V_{\bf k}$, $\omega_{\bf k}\longrightarrow\Omega_{\bf k}$ in the corresponding equations. The process is started with an initial guess for $Z_{\bf k},\Sigma_{aa}({\bf k},0)$ and by using Eq.(21) we find corrected excitation energy. This is repeated until convergence is reached. Using the final values for excitation spectrum we can calculate the dynamical and static spin structure factors by Eqs.(31,33,34,36). We discuss the numerical results for thermodynamic properties in the field induced spin polarized regime where energy spectrum of spin model hamiltonian includes a finite energy gap between ground state and first excited state. Therefore as long as excitation spectrum $\Omega_{{\bf Q}_{0}=(0,4\pi/3)}$ has non zero values, the system preserves its gapped spin polarized phase. In the following numerical results for both static and dynamical spin susceptibilities the wave vector has been fixed at generic wave vector ${\bf q}_{0}=(2\pi/\sqrt{3},0)$.

The effect of DM interaction strength, $D$, on critical magnetic field has been studied in Fig.(2). In this figure, energy gap ($\Delta$) versus magnetic field $g\mu_{B}B/J$ for different values of Dzyaloshinskii-Moriya interaction strength $D/J$ for $J^{\prime}/J=0.2$ by setting $k_{B}T/J=0.05$ has been plotted. Fig.(2) shows that the energy gap vanishes as the magnetic field approaches the critical value for $g\mu_{B}B_{c}/J$. For all values of $D/J$, the gap vanishes at the critical point $g\mu_{B}B_{c}/J$ where the transition from gapped spin liquid phase to the gapless one occurs. Moreover the critical field tends to higher value with increase of $D/J$ according to Fig.(2).

In Fig.(3), the dynamical transverse spin structure factor, i.e. $Im\chi_{+-}({\bf q}_{0},\omega)$, has been plotted as a function of frequency $\omega/J$ for different next nearest neighbor coupling exchange constants, namely $J^{\prime}/J=0.0,0.1,0.3$, at fixed magnetic field $g\mu_{B}B/J=6.0$ for $D/J=0.2$. For each value for $J^{\prime}/J$, there are two sharp peaks in imaginary part of transverse dynamical spin susceptibility at low and high frequencies. These frequency positions of sharp peaks describe the transverse magnetic excitation spectrum of the localized electrons on honeycomb structure. As shown in Fig.(3), the lower frequency peak have lower height in comparison with the higher frequency peak for each value of $J^{\prime}/J$. The height of each peak corresponds to the intensity of scattered neutrons from the sample which describes the inelastic scattering cross section of neutron beam. Moreover the frequency of each peak is equal to the difference between scattered neutron beam energy and incident neutron energy beam. This energy difference is absorbed by the localized electrons which causes to magnetic excitation of transvesre components of spins. According to Fig.(3), it is clearly observed that frequency position of peaks in the imaginary part of transverse spin susceptibility, $Im\chi_{+-}({\bf q}_{0},\omega)$, goes to lower values with increase of next nearest neighbor exchange constant $J^{\prime}/J$. It can be understood from this fact that the energy gap between ground state and first excited state decreases with $J^{\prime}$ and consequently the magnetic excitation appears at lower frequency. Moreover the Fig.(3) implies the intensity of scattered neutron beam from the sample is independent of $J^{\prime}/J$ and the heights of peaks in transverse spin structure factor are the same. The emergence of peaks in the dynamical spin susceptibilities of Heisenberg chain have been shown for magnetic fields below threshold one using quantum Monte Carlo method based on the stochastic series expansion and maximum entropy methodbrenig . At threshold magnetic field, the energy gap closes so that long range magnetic ordering develops for magnetic field below threshold one. Based on Fig.(2), the normalized threshold magnetic field is about $g\mu_{B}B_{c}/J\approx 3$. The estiamtion value for exchange coupling constant is around $J\approx 1meV$. Thus one can obtain the threshold magnetic field $B_{c}\approx 60$ Tesla. However the magnetic field for theoretical studies of Heisenberg model Hamiltonian on lattice structures has been considered as normalized value $g\mu_{B}B/J$langer ; zotos . In theoretical studies the qualitative behavior of physical parameters in terms of magnetic field is more important. Also the normalized magnetic field value reported in the theoretical works is in the region $0<g\mu_{B}B/J<5.0$. Also dynamical structure factor of two dimensional antiferromagnet in high fields has been investigated within spin-wave theory using self-consistent approach beyond the $1/S$ approximationfuhrman . The emergence of sharp peak in the frequency behavior of spin structure factor has been predicted in this work in analogy to our work.