Bosonic Dirac Materials on a honeycomb antiferromagnetic Ising model

TL;DR

This paper demonstrates that a honeycomb antiferromagnetic Ising model can host Bosonic Dirac points in its magnon spectrum, making it a potential prototype for Bosonic Dirac materials controlled by external fields.

Contribution

It shows that the Ising model on a honeycomb lattice can host Bosonic Dirac points, expanding the understanding of bosonic Dirac materials and their realization in simple models.

Findings

Identification of Bosonic Dirac points in the magnon spectrum.

Massless magnons near Dirac points with linear dispersion.

External fields control the emergence of Dirac points.

Abstract

Motivated by the recent proposal of Bosonic Dirac materials (BDM), we revisited the Ising model on a honeycomb lattice in the presence of the longitudinal and transverse fields. We apply linear spin-wave theory to obtain the magnon dispersion and its degenerated points. These special degenerated points emerge on the excitation spectrum as a function of the external fields and can be identified as Bosonic Dirac Points (BDP). Since that, in the vicinity of these points the Magnons becomes massless with a linear energy spectrum as well as insensible in relation to weak impurity, exactly as it occurs with a Fermionic Dirac point. We also have calculated the quantum and thermal fluctuations over the ground state of the system using Effective Field Theory. Our results point out that this simple model can host Bosonic Dirac points and therefore is a suitable prototype to build a Bosonic Dirac material only controlled by external field.

Figures

Equations

$$\mathcal{H}=J\sum _{<i,j>}{S}_i^z{S}_j^z-h_x\sum _i{S}_i^x-h_z\sum _i{S}_i^z,$$

$$G(k,\omega )=g_0(k,\omega )+g_0(k,\omega )T_{imp}{g}_0(k,\omega ),$$

$${\tilde{S}}_i^{zp}$$

$${\tilde{S}}_i^{xp}$$

$${\mathcal{H}}^{LSWT}$$

$$E_c$$

$$\{\begin{array}{ll}& \tilde{A}=JzS\cos ({\varphi }_A)\cos ({\varphi }_B)+(h_x\sin ({\varphi }_A)+h_z\cos ({\varphi }_A))\\ & \tilde{B}=JzS\cos ({\varphi }_A)\cos ({\varphi }_B)+(h_x\sin ({\varphi }_B)-h_z\cos ({\varphi }_B))\\ & C_k=Jz\frac{S}{2}\sin ({\varphi }_A)\sin ({\varphi }_B)\gamma (k)\end{array},$$

$$H^{LSWT}=E_g+\frac{1}{2}\sum _k({\omega }_1(k){\alpha }_k^{†}{\alpha }_k+{\omega }_2(k){\beta }_k^{†}{\beta }_k)$$

$$E_g(T=0)=E_c+E_q$$

$${\omega }_{1,2}(k)$$

$${\bar{m}}_p^z$$

$$\omega (\mathbf{\Gamma })\sim D\sqrt{{\Delta }^{2\nu z}+k^{2z}},$$

$$\omega (\mathbf{K})\sim \tilde{A}\pm \frac{JzS\sin ({\varphi }_A{)}^2}{4}\delta k.$$

$$\mathcal{H}=E_q+\frac{1}{2}\sum _{k=\{k_x,k_y\}}{\psi }^{†}\bar{\mathcal{H}}(k_x,k_y)\psi ,$$

$$\mathcal{\bar{H}}(k_{x},k_{y})=\left(\begin{array}[]{cccc}\tilde{A}&C_{k}&0&C_{k}\\ C_{-k}&\tilde{B}&C_{-k}&0\\ 0&C_{k}&\tilde{A}&C_{k}\\ C_{-k}&0&C_{-k}&\tilde{B}\\ \end{array}\right).$$

$${\bar{\mathcal{H}}}_{I,J}(k_y)=\frac{1}{N_x}\sum _{k_x}{e}^{i{k}_x(m-m^{\prime })}{\bar{\mathcal{H}}}_{I,J}(k_x,k_y),$$

$${\psi }^{†}=$$

$$\tilde{A}$$

$$\tilde{B}$$

$$C_k$$

$$C_{-k}$$

$${\psi }_{zz}^{†}=\{b_{0,k_y}^{†},a_{1,k_y}^{†},b_{1,k_y}^{†}\cdots a_{m,k_y}^{†},b_{m,k_y}^{†},a_{m+1,k_y}^{†},b_{0,-k_y},a_{1,-k_y},b_{1,-k_y}\cdots a_{m,-k_y},b_{m,-k_y},a_{m+1,-k_y}\},$$

$$\bar{\mathcal{H}}(k_y{)}_{zz}$$

$$I_{-}=\left(\begin{array}[]{cc}I_{2L_{x}+2}&0_{2L_{x}+2}\\ 0_{2L_{x}+2}&-I_{2L_{x}+2}\end{array}\right),$$

$${\mathcal{H}}_{1p}=\left(J\stackrel{z}{\sum _{\delta }}{S}_{(1+\delta )p^{\prime }}^z-h_z\right)S_{1p}^z-h_x{S}_{1p}^x,$$

$$m_p=⟨\frac{h_z-a_{1p}}{\sqrt{(h_z-a_{1p}{)}^2+h_x^2}}\tanh \beta \sqrt{(h_z-a_{1p}{)}^2+h_x^2}⟩,$$

$$m_p=⟨\prod _{\delta \ne 0}^z({\alpha }_x+S_{(1+\delta )p^{\prime }}^z{\beta }_x)⟩{F(x)|}_{x=0},$$

$$F(x)=\frac{h_z-x}{\sqrt{(h_z-x{)}^2+h_x^2}}\tanh \beta \sqrt{(h_z-x{)}^2+h_x^2},$$

$$⟨S_{iA}^z{S}_{jB}^z\dots S_{lA}^z⟩\backsimeq m_A{m}_B\dots m_A,$$

$$m_p=\sum _{q=0}^z{A}_q(T_N,h_z,h_x)m_{p^{\prime }}^q,$$

$$A_q(T_N,H_z,h_x)=\frac{z!}{q!(z-q)!}{\alpha }_x^{z-q}{\beta }_x^q{F(x)|}_{x=0}.$$

$$m_A=X_0(T_N,h_z,h_x,m_0)+X_1(T_N,h_z,h_x,m_0)m_s,$$

Full text

I Introduction

The Ising model has been largely applied in statistical physics, quantum field theory, economy and biophysics Ising (1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 39, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26). In condensed matter physics, this model can be used to study magnetic phases and its phase transitions. In particular, to study quantum phase transitions in the Ising model is necessary to add an external magnetic transverse field. Generally, the transverse field induces a quantum phase transition from a ferromagnetic (F) or antiferromagnetic (AF) phase to a polarized one, In which the majority of the spins is along the direction of the field. In this case, the ground state before the transition (F/AF) can not be adiabatically connected with the ground state of the polarized phase.

Another intriguing result is the existence of Bosonic Dirac Materials (BDM). As pointed by Balatsky et al  Balatsky (5), the concept of fermionic Dirac materials can be extended to a bosonic version BDM. To perform this concept to a BDM is enough to show the following points; A linear characteristic of the energy dispersion in the vicinity of a point where the energy gap goes to zero; The robustness of these degenerated point against a weak disorder potential and a relative stability of these points against higher expansions of the spin wave theory.

In this work, we will demonstrate that a simple Ising model defined on the honeycomb lattice and subjected to the presence of transverse and longitudinal external field fulfils all these rules and therefore can be used as a prototype for BDM. Theoretically, the model can be written as

[TABLE]

such that, $J$ is the nearest neighbor exchange interaction, $S_{i}^{z}$ is the $z$ spin component in the site $i$, $h_{x}$ and $h_{z}$ are the transverse and longitudinal magnetic fields, respectively. The Eq. 1, was subjected to several methods and different lattices KinrossSachdev (2, 9, 10, 39).

Recent works point out that, the F and AF Heisenberg model on a honeycomb lattice exhibits degenerated points in the first Brillouin zone. At these points, the spectrum can be linearized and associated with BDP. Beside that, these BDP are robust against higher order Magnon-Magnon interaction, which are responsible to make a shift in the spectrum Balatsky (5).

The main difference between our results and the results from ref. Balatsky (5) is the presence of the Bosonic Dirac points only in the polarized phase. Therefore, is possible turn on or turn off these Bosonic Dirac Point by adjusting the external magnetic fields.

We find that by a fine tuning of the magnetic field, is possible to generate BDP located at the corners of the first Brillouin zone of the honeycomb lattice. From this, follows that is possible to transform a non-Bosonic Dirac phase into a Bosonic Dirac phase and vice verse. This is fascinating because at these points the magnons are massless, insensible to a weak disorder scattering potentials and can be propagated with very high speed $v_{B}$. A directly application of this result can be found in the generation of efficient spin transport Pires (27). It is worth to note, that some mechanism of the creation, localization and vanishing of these Dirac-like points, has attracted theoretical and experimental researchers Rechtsman (28, 29, 30, 31, 32, 33, 34, 35, 36, 37).

In addition, we determined the effects of the competition between the transverse and longitudinal fields, the influence of on-site potential disorder over the BDP, edge states behavior and temperature effects.

For this purpose, we used two different methods, linear spin-wave theory (LSWT) and effective-field theory 1 (EFT-1).

The LSWT has been used to build the ground state and the dispersion in the real and momentum space. In the real space, with open boundary conditions promoting a zig zag termination, we showed the edge states for a clean and disorder affected systems. On the other side, EFT-1 has been used to compute zero and finite temperature effects over the sublattice magnetization, as a generalization of the pervious works Creswick (13, 14, 15).

Experimentally, in order to test our results we suggest the scanning tunneling microscopy (STM) technique. Firstly, STM can be used to create a $2D$ honeycomb magnetic lattice Gomes (38) by depositing atoms or molecules in substrates. Secondly, the STM is able to change continuously the magnetic exchange interaction simply manipulating the position of the absorbantes (atoms, molecules etc..) which are deposited in a subtraction. Such that, for the present purpose an external field can be added to verify the controlling over the BDP.

The paper is organized as follows; In the section II, we will discuss about the relevant ground states; The excitation spectrum critical exponents, Bosonic Dirac Point will be discussed in the section III, while the impurity effects and edge states has been calculated in the sections IV and V, respectively. The magnetic properties and finite temperature effects will be showed in section VI. Finally, in the last section VII, we summarize our results.

II Phase Diagram and ground state

Using LSWT and EFT-1 (appendix A and D, respectively), we determined the phase diagram of the system in plane $h_{x}-h_{z}$, see Fig. 1. The critical lines separates a canted antiferromagnetic (CAF) (region below the critical line) phase from polarized phases (PP) (region above the critical line), see Fig. 2. The black and red lines were obtained by LSWT, while the blue line is a result from EFT-1.

We obtained the black line minimizing the classical energy $E_{g}=E_{c}$ with respect to angles $\phi_{A}$ and $\phi_{B}$ (LSWT-1). On the other side, adding the quantum correction $E_{q}$ in the the energy, we get the red line (LSWT-2). Note that, at the critical line and in the region above the critical line $\phi_{A}+\phi_{B}=\pi$ (Fig. 2 (c)).

Still about the angle spin configuration, see the quantum phase transition that occurs at $h_{x}=1.5$ for $h_{z}=0$ (red big solid circle in Fig. 1). Above $h_{x}\geq 1.5$ the spins are totally along the x-axis. Following, for $h_{x}=1.5$ and increasing the longitudinal field from zero to $h_{z}=0.5$ (red big solid circle) the spin configuration of the polarized phase is no more totally along the x-axis, getting an effective angle $\phi$ as in Fig. 2 (c). Along the critical line, this angle changes and finally for $h_{x}=0$ and $h_{z}=1.5$, $\phi=0$.

The results from LSWT-2 are qualitatively equal to LSWT-1, in the sense that it describes a continuous transition between CAF and P phases. However, one can observe a considerable reduction of the CAF region for LSWT-2 in comparison to LSWT-1, see Fig 1 (red line and triangle). This is a new result, in particular in the case of a Honeycomb lattice and corroborates with the idea that the spin wave theory tends to overestimate the ordered phase region in a first approach.

The critical line from EFT-1 (blue line and triangles) was also included Fig.1 and seems to be in great agreement with the critical line obtained from LSWT-2 (red squares). Here, the EFT-1 method takes into account a cluster of one spin interacting with their vicinity. We considered a staggered magnetization as a order parameter and the quantum critical point can be founded throughout the Eqs. 31- 32. It is worth to communicate that LSWT-1 provides results quantitative equal to a first approximation of Bogoliobov Mean Field Theory (MFA-1) for a cluster with one spin Neto2 (15).

III Linear Behaviour, Bosonic Dirac points and critical exponents

The linear behaviour can be identified around the special points called by $\Gamma=(0,0)$ and $K=(\frac{2\pi}{3a},\frac{2\pi}{3\sqrt{3}a})$ within the first Brillouin zone, see Fig. 3. Only the $K$ points are possible candidate to be “Bosonic Dirac points”, following the same nomenclature adopted by the Balatsky et al. Balatsky (5).

The low energy description of the magnon spectrum provides a linear dispersion around $\Gamma$ equal to $\omega(\mathbf{\Gamma})\sim D\sqrt{\Delta^{2\nu z}+k^{2z}}$. Where the gap function $\Delta$ goes to zero at the critical point as well as $\nu$ and $z$ are the correlation and dynamical critical exponents. Note that, $D$ and $\Delta$ are independent of $k$, they are given by $D=\sqrt{\frac{{\tilde{A}JzS{\sin(\phi_{A})}^{2}}}{4}}$ and $\Delta=\frac{\sin(\phi_{A})^{2}}{D^{2}}\left(h_{x}-h_{xc}\right)$ (since we choose to fix $h_{z}$), respectively. The critical value can be written as $h_{xc}=\frac{JzS-h_{z}\cos(\phi_{A})}{\sin(\phi_{A})}$ where the gap vanishes, as shown in more details in appendix B. These calculations leading to $\nu=\frac{1}{2}$ and $z=1$ for the critical exponents, the last ensure the linear behaviour.

For the critical exponent of the correlation length, we found $\nu=\frac{1}{2}$ as expected for a mean field approach. It is worth mentioning, for a 2D Ising model with a transverse field the value gets $\nu=0.63$ Elliott (39, 40).

Following, around the $K$ the dispersion gets a linear behavior as $\omega(\mathbf{K})\sim\tilde{A}\pm\frac{3}{4}JS\sin(\phi_{A})^{2}\delta k$. It is worth to note that in this case the gap function is $\Delta=(\omega_{1}(\mathbf{K})-\omega_{2}(\mathbf{K}))$, in particular, for $h_{z}=0$, $\Delta=(h_{xc}-h_{x})^{\nu z}$ Stinchcombe (40, 41) with $\nu z=1$.

The $z=1$ result for the Ising model in a presence of transverse, is known and works for any regular lattice since the dimension of the system be an integer SSachdev (42). In addition, this is also independent of the quantum critical point as demonstrated for the $\Gamma$ and $K$ points.

We are reporting in this work that theses linear behaviours appear in spectrum as a function of the external fields. The dominant set of critical exponents is dictated by the soft point at $\Gamma$. The number and location of these special points depends on criticality and the magnetic phase, see Fig. 4 and Fig. 1.

These results can be summarized as; In the CAF phase there is no Bosonic Dirac Points, sice the system is gapped. At the critical line, the energy dispersion becomes gapless and six Bosonic Dirac points emerge at corners of the Brillouin zone ($K$) and there is a linear soft mode at the center ($\Gamma$). These six Bosonic Dirac point at the corners of BZ persist in the polarized phase.

Our report is quite different from the case discussed in ref. Balatsky (5), since the existence of the Bosonic Dirac Points depends on the magnetic phase. For the present case, the controlling of these Dirac points resides in a simple adjusting of the magnetic fields instead of the exchange interaction $J$.

In Fig. 4, we explore the projection of the magnon excitation along of the path defined by $\Gamma-M-K-\Gamma$. This figure shows the spots of the Bosonic Dirac points in the BZ and how they are affected by the magnetic fields.

The color scheme follows the big solid circles in Fig. 1 for each pair $h_{x}-h_{z}$. The red line in the Fig. 4 a), is the quasiparticle excitations at the critical point $(h_{x}=1.5,h_{z}=0)$. For a fixed $h_{x}=1.5$ and varying $h_{z}$ from zero to 0.5, we get the point $(h_{x}=1.5,h_{z}=0.5)$ (polarized phase), where the dispersion is linear only around the $K$ points. The same behavior occur at the point $(h_{x}=1.5,h_{z}=1.0)$ (green curve). A similar behavior can be observed, but for a fixed $h_{x}=1.0$, see Fig. 4 b). Therefore, the Bosonic Dirac points seems to be a particular characteristic of the polarized phase.

IV Impurity effects and impurity resonance

Are these degenerate Bosonic Dirac points robust against impurities? To answer this question, we investigated the effects of an impurity scattering. The impurity simulates a local defect created by any removing/modification in the spin over the honeycomb lattice. For this reason, we considered a defect given by $H_{imp}=\sum_{i}V_{0}a^{\dagger}_{i}a_{i}$.

In order to present the physical conditions that the impurity scattering should obey to generate an impurity resonance, exactly at the Bosonic Dirac point $K$ and linear point $\Gamma$, we need to calculate the Green’s function of the system from the Dyson equation. The absence of an impurity resonance at these points implies robustness against the scattering potential $V_{0}$. So, this robustness need to be checked in order to make the association with “Bosonic Dirac point”. We write the impurity Green’s function using the T-matrix approach. The single magnon Green’s function for our problem is given by

[TABLE]

where, $g_{0}=(\omega+i\delta-\mathcal{H})^{-1}$, $T_{imp}=\frac{1}{2(V_{0}^{-1}+\bar{g}_{0})}\left(\sigma_{0}\otimes\sigma_{0}+\sigma_{0}\otimes\sigma_{z}\right)$, $\sigma_{i}$ are the usual Pauli matrix and $\bar{g}_{0}(\sigma_{0}\otimes\sigma_{0})=\frac{1}{N}\sum_{k}g_{0}(k,\omega)$.

To determine the impurity effects over the phase diagram in Fig. 1, we chose two points A and B; The A point is given by $(h_{x}=1.0,h_{z}=1.28)$, which corresponds to a critical point. The B point $(h_{x}=1.5,h_{z}=1.0)$ within the polarized region. Using the A point, we calculated the density of states for a clean case $\rho_{0}(\omega)=-\frac{1}{\pi}Im(g_{0}(\omega))$ (dashed blue lines of Fig. 5) and the impurity correction $\delta\rho(\omega)=-\frac{1}{\pi}Im(g_{0}(\omega)T_{imp}g_{0}(\omega))$( solid lines of Fig. 5). So, in the Fig. 5 an impurity resonance might emerges only in the case of very strong values of scattering potential $V_{0}$. This means that the correction in $\rho_{0}$ is precisely zero for low values of potential scattering $V_{0}$. Note yet, that the pick of $\delta\rho$ approaches from the resonance limit only when $V_{0}$ goes to infinity. This behaviour was observed along all critical line, ergo the critical line is robusts against this kind of defect.

Furthermore, the Fig. 6 shows the effects of an impurity scattering over de Bosonic Dirac point $K$ in the polarized phase. Again, low values of the impurity potential $V_{0}$ can not generate a resonance impurity. These results can be used to guarantee that these special degenerated points $\Gamma$ and $K$ are robust against a weak impurity potential in a similar way what as it occurs with their cousin, Fermionic Dirac points.

In the strong limit of scattering potential, the impurity resonance emerges only in the sublattice that does not hosts the scattering center. As discussed in ref. Balatsky (5), the scattering induces Friedel oscillations and the local magnon density $\rho(r,\omega)$ around the impurity describes waves emanating from the scattering center with the standard asymptotic decay proportional to $r^{-1}$.

V Real space and edge states

An interesting feature that arises in the finite honeycomb lattice is the presence of edges states. In this section we are devoted to understand the nature of these edges states and how a disorder potential can destroys it. The presence of the Bosonic Dirac Points at $K$ points in this finite lattice is also a relevant point to be observed.

We have written the Hamiltonian eq A in the real space with open boundary conditions along the $x$ direction, while in the $y$ direction we considered a periodic boundary condition (see appendix C). Thus, the lattice becomes a kind of a finite cylinder with zigzag terminations, in a very similar way to a nanoribbon in graphene context. The edge states are localized around the zigzag terminations of this cylinder, the highlighted red region, as depicted in Fig. 7 d).

Without disorder, the eigenvalues of the clean Hamiltonian in real space Eq. 19 are shown in the Fig. 7 a), b) and c) for $L_{x}=100$. Where $L_{x}$ counts the number of unit cell, containing a pair of sublattice A and B, along the finite $x$ direction. In Fig. 7, we have the bulk bosonic states (blue lines) and the zigzag edge states (red lines) of all magnetic phases. The edges states from CAF and polarized phase are shown in the parts a) and c), respectively. In the part b), we showed how these edges states behaves at the critical line. The location of the edge states are depicted in the part d). They are majority distributed around the terminations of the lattice as indicated by the red gradient color.

A coupling with a disorder potential can simulates an interaction with a substrate. The adding of this disorder potential allow us to investigate its effect over these Bosonic Dirac points and edge states within the polarized phase. Also, providing an another way to determine the robustness of the BDP against disorder inside a finite lattice study.

After this disorder calculations, we present the result of Fig. 8 for an ensemble with 30 realisations. As we can see, a weak on-site potential disorder for both bulk and edges, and only concentrated at the edges, respectively Fig. 8 a) and b), can not open a gap around the points K. So, even in this finite lattice study, these $K$ points are insensible in relation to this kind of disorder and so far they prove to be robust.

VI The temperature effects

Beyond the zero temperature results and in order to design a more complete study, we also investigate the behavior of the ground state in the presence of thermal fluctuations. Through the sublattice magnetization $m^{z}_{A,B}$ as a function of the reduced temperature $T=k_{B}T/J$ and the magnetic fields we developed nonzero phase diagrams.

Taking into account only the critical behavior, we present the critical reduced temperature as a function of the $h_{x}$ for several values of $h_{z}$ obtained by LSWT, and vise versa in Fig. 9 and Fig. 10, respectively. These results are able to be compared with Fig. 11 and Fig. 12, that comes from EFT-1.

Looking to the Fig. 9, we found that the critical temperature $T_{c}$ of the sublattice magnetization $m^{z}_{A}$ is identically to the melting temperature. So, above $T_{c}$ the system is paramagnetic without any sublattice magnetization. The EFT-1 results shown in the Fig. 11 are qualitative equal to results of the LSWT-1, therefore both methods are in excellent agreement.

On the other hand, analyzing the figures Fig. 10 and Fig. 12, which shows the behavior of the $T_{c}$ as a function of the $h_{z}$ for a fixed $h_{x}$, we can see a substantial difference between these two methods. Notice that, while the EFT-1 method found a monotonic behaviour for magnetization, the LSWT-1,2 suggests a more intriguing situation. The Fig. 10 indicates that the critical temperature of the two sublattice magnetizations possesses different behaviors, which depends on whether $h_{z}$ is less or greater than critical value $h_{xc}$, where the sublattice magnetizations go to zero. Below the critical point, the sublattice magnetization increases again and this occur because the polarized phase is always induced by the increasing of the $h_{z}$.

VII Conclusion

In this paper, we investigated the quantum phase transitions and the existence of Bosonic Dirac points in the Ising model on a honeycomb lattice. The quantum phase transitions are induced by external fields.

Using two methods, LSWT and EFT-1, we determined the phase diagram in plane $(h_{x},h_{z})$ at zero temperature. From LSWT, we determine the bosonic quasiparticle excitation. Depending of the magnetic phase, degenerated points can emerge within the spectrum at the $K$ points in the Brillouin zone. These points always occur along the critical line or within the polarized magnetic region. We demonstrated that around these special points the spectrum is linear and is robust against a weak impurity and disorder. The dynamical critical exponent was also investigated leading to $z=1$, that indicates a linear crossing around these points. These results often characterize Fermionic Dirac Materials and therefore can be used to classify our model as a Bosonic Dirac Material.

Following this, it is fascinate how the present simple model (inside the polarized phase) became a very strong candidate to describe a Bosonic Dirac material. As we see, this phase can be reached by a fine tune of the longitudinal and transverse fields. Which is very interesting, because in principle it makes possible to turn a non Bosonic Dirac material into a Bosonic Dirac material only adjusting the external parameters like magnetic fields. Beside that, the Bosonic Dirac materials should be robust against the magnon-magnon interaction as pointed by ref. Balatsky (5), from this perspective the existence of this Bosonic Dirac points depends only on the symmetry of the honeycomb lattice. Furthermore, differently from the graphene fermionic case, our bosonic quasiparticle excitation can also exhibit one additional linear behavior at the center of the BZ, a soft mode.

In the real space we found edge states for a lattice with zigzag termination, these edge states are a peculiar characteristic of the honeycomb lattice and can be observed in the AF and polarized phase. They proved to be robust against potential disorder and provides one more evidence for a Bosonic Dirac Material description.

To determine the temperature dependence and the possible quantum corrections, we calculating the magnetic properties using the LSWT-2 and EFT-1, which allow us to conclude that the methods are in great agreement.

It is worth to highlight yet the importance of the Bosonic Dirac points due its connection with production of spin wave excitations Pires (27), which can be used to breed efficient spin currents. Accordingly, these excitations could find application inside of Spintronics and optical lattices Offidani (43, 44, 45).

The real compound $Ba_{2}CoTeO_{6}$ can be described by an effective Ising model on a honeycomb lattice, and accordingly the present study, is a strong candidate to be a Bonsonic Dirac material Chanlert (51). Despite of the fact that this compound presents an approachable next-nearest-neighbor interaction $J_{2}$, the symmetry of honeycomb is preserved. Another promise compound was syntetized by Okabe (52). It is a single crystal of the verdazy radical which also presents a Neel order that can be described by a spin 1/2 antiferromagntic Heisenberg model in a honeycomb lattice.

We expect to motivate experimental efforts to study the nature of the Magnon Bosonic Dirac points and its relation with the quantum phase transitions in this simple model.

VIII Acknowledgment

We would like to thank you the prof. Mucio Continentino and prof. Octavio D. Rodriguez Salmon for all assistance as well as the Brazilian agencies CNPq, Capes and Fapeam for providing all supporting.

Appendix A The linear spin wave approach

The LSWT often presents some advantages over the usual mean field approximations, once the classical ground state (GS) and the origin of fluctuations are identified, this treatment provides the calculation of the dispersion relation that can be connected with neutron spectroscopy measurements and so on taking to some experimental background Pearson (46).

For the present Hamiltonian Eq. 1, in the absence of transverse and longitudinal fields, the classical ground state is antiferromagnetic. However, when the fields are turned on, the transverse field is responsible for the quantum fluctuation Dutta (47). Previous works has been applied LSWT in the presence of a transverse fields ColletaMila (48, 49). In this case, for low longitudinal fields, the transverse field deviates the spin into the $x$ direction. For this reason, the ground state, which should be considered is the canted Antiferromagnetic state (CAF). As result, these deviations produce a local rotation of the spin configurations over the $y$ axis, described by $\tilde{S}^{p}_{j}=S_{p}\mathcal{R}_{y}(S_{p}\phi_{p})S^{p}_{j}$, where $p$ accounts the sublattice. Namely, the sublattices A and B stands for spin up and down. The ground state configuration is obtained by the minimum value for energy in respect to the angles $\phi_{A}$ and $\phi_{B}$. In order to characterize the expected classical ground states as a function of the spin angle configurations, we defined the canted CAF phase when $\phi_{A}\neq\phi_{B}$, and polarized when $\phi_{A}+\phi_{B}=\pi$.

After performing the rotation above on the Hamiltonian Eq. 1, we can introduce the linear spin wave approach theory following Hostein-Primakoff transformation

[TABLE]

where, the $a^{\dagger}_{i,p}$($a_{i,p}$) is a bosonic creator (annihilator) operator of the vacuum state $|v\rangle$ and $S$ is the magnitude of the spin Mila (49). Substituting the spin wave mapping, Eq. 3, in the rotated Hamiltonian and performing a Fourier transformation, the quadratic Hamiltonian in a momentum space is finally obtained as

[TABLE]

The energy of the classical state $E_{c}$ and the coefficients are given by

[TABLE]

[TABLE]

with $\gamma(k)=\frac{1}{z}\sum_{\delta}e^{k\delta}$ being the structure factor, $\delta$ a vector that localizes all first neighbors from a site i and $z$ is the number of neighbors.

So, in order to find the quasi-particles excitations and spectrum dispersion of the Magnons, we used the Bogoliubov-Valatin transformation Mila (49) to obtain the following diagonal form of the Hamiltonian

[TABLE]

such that,

[TABLE]

is the ground state energy with the classical part $E_{c}$ and quantum correction equal to $E_{q}=-\frac{N}{2}(\tilde{A}+\tilde{B})+\frac{1}{2}\sum_{k}(\omega_{1,k}+\omega_{2,k})$.

The quasiparticles energies is given by

[TABLE]

where $\Psi_{k}=({{{\tilde{A}}}^{2}-{{\tilde{B}}}^{2}})^{2}+16{\tilde{A}}{\tilde{B}}{|C_{k}|^{2}}$.

In special, for the honeycomb lattice $z=3$ and $\gamma(k)=\frac{1}{3}(\,{{\rm e}^{\frac{i}{2}\left({\it kxa}+\,\sqrt{3}{\it kya}\right)}}+\,{{\rm e}^{\frac{i}{2}\left({\it kxa}-\,\sqrt{3}{\it kya}\right)}}+\,{{\rm e}^{-i{\it kxa}}})$, with all summation and integrals defined over the first Brillouin zone. From here on, the lattice constant assumes $a=1$.

The sublattice magnetization along the rotated quantized axis can be calculated by the following expression

[TABLE]

with $N_{uc}$ being the number of unity cells. In order to build the order parameter for each phase along the phase diagram in plane $h_{x}$-$h_{z}$, is convenient to define the magnetization along the old axis, i.e. $m^{x}_{p}=\bar{m}_{p}^{z}\sin(\phi_{p})$ and $m^{z}_{p}=\bar{m}_{p}^{z}\cos(\phi_{p})$.

Appendix B The low-energy description

An effective low-energy description of the Hamiltonian can be obtained by a series expansion of $\gamma(k)$ around the critical points $\mathbf{k}_{0}$. These points are located on the center $\mathbf{\Gamma}$ and corners $\mathbf{K}$ of the Brillouin zone. At the critical point $\phi_{B}=\pi-\phi_{A}$, that implies $\tilde{A}=\tilde{B}$ in Eq. 9. This constraint does not depend of BZ point choice.

Near of the $\mathbf{\Gamma}$ point, $\gamma(\mathbf{\Gamma})\sim(1-\frac{k^{2}}{4})$, for $k^{2}=k_{x}^{2}+k_{y}^{2}$. Substituting this expression in Eq. 9 (with $(+)$ within the first square), we get

[TABLE]

the coefficient is given by $D=\sqrt{\frac{{\tilde{A}JzS{\sin(\phi_{A})}^{2}}}{4}}$. While the gap function can be written as $\Delta=\frac{\sin(\phi_{A})^{2}}{D^{2}}\left(h_{x}-h_{xc}\right)$, since we choose to fix $h_{z}$ and looking for the critical value $h_{xc}=\frac{JzS-h_{z}\cos(\phi_{A})}{\sin(\phi_{A})}$ where the gap vanishes. The Eq. 11 was written as a function of the critical exponents, leading to $\nu=\frac{1}{2}$ and $z=1$ for the critical exponents of spatial and temporal correlation lengths. Note, $\omega(\mathbf{\Gamma})$ is linear “if and only if” $\Delta=0$, which is true for instance, at the critical point $h_{x}=\frac{3}{2}J$, $h_{z}=0$ and $\phi_{A}=\frac{\pi}{2}$, that leads $h_{xc}=JzS=\frac{3}{2}J$ and the gap function goes to zero. Is useful to realize that the dynamical critical exponent $z=1$ indicates a linear behaviour of the energy dispersion at the critical point and near of the soft point $\mathbf{\Gamma}$.

Now, performing the same but around of the $\mathbf{K}$ points, the series expansion for $\gamma(k)$ provides $\gamma(\mathbf{K}+\delta k)\sim(-\frac{\sqrt{3}}{4}\delta k_{x}-\frac{1}{4}\delta k_{y})+i(\frac{1}{4}\delta k_{x}-\frac{\sqrt{3}}{4}\delta k_{y})$, where $\delta k_{x}$ and $\delta k_{y}$ are low increments of $k$ in the vicinity of $K$. Substituting $\gamma(k)$ in Eq. 9, yields

[TABLE]

In this case, the gap goes to zero at finite value of the energy dispersion given by $\tilde{A}$. Therefore the gap function can be written as $\Delta=2(\omega(\mathbf{K})-A_{c})$, which renders $\nu=1$, $z=1$ and as a consequence an linear behaviour conclusion for $\omega(\mathbf{K})$. This result is one of the characteristics of a Bosonic Dirac Materials.

Appendix C Hamiltonian in an real space

In the momentum space $k=\{k_{x},k_{y}\}$, we have

[TABLE]

for a basis $\psi^{{\dagger}}=\{a^{{\dagger}}_{k},b^{{\dagger}}_{k},a_{-k},b_{-k}\}$ and

[TABLE]

After a Fourier transform to real space in $x$-axis, we need to define the $m$ index that runs over the $L_{x}$ pair of sublattice $A$ and $B$ along the now finite axis. Following this, each matrix element in 14 becomes a matrix of $L_{x}$-order, that can be written as

[TABLE]

where, $I$ and $J$ identify a matrix element from 14. Therefore, to complete bulk matrix in $x$-axis real space, $\mathcal{\bar{H}}(k_{y})$ must be of $4L_{x}$ order and the basis is increased as

[TABLE]

Accordingly 15, each matrix element in 14 becomes a matrix as

[TABLE]

for, $\tilde{C}=J\frac{S}{2}\sin(\phi_{A})\sin(\phi_{B})$ and $D_{k_{y}}=2\tilde{C}\cos(\frac{\sqrt{3}}{2}k_{y})$.

In order to reproduce the Zig-Zag edge configuration of the honeycomb lattice, we rewrite the $\mathcal{\bar{H}}(k_{y})$ in an appropriate basis 18, starting (left-edge) with a sublattice B and finishing (right-edge) with the sublattice A. The pair of operators $\{b^{{\dagger}}_{0,k_{y}},b_{0,-k_{y}}\}$ and $\{a^{{\dagger}}_{m,k_{y}},a_{m+1,-k_{y}}\}$, represent the left-edge and right-edge, respectively. The matrix $\mathcal{\bar{H}}(k_{y})$ on the Zig-Zag basis must assume the matrix form of 19 with $2(2L_{x}+2)$ order

[TABLE]

[TABLE]

where, $\sigma_{0}$ and $\sigma_{1}$ are the Pauli matrix. Now, the $m$-index must runs from zero to $L_{x}$, to take into account the pair formed by the left and right edges.

Since 19 represents a bosonic system, to calculate their eigenvalues the diagonalization process should be performs over the dynamic matrix $M_{D}=I_{-}\mathcal{\bar{H}}(k_{y})_{zz}$, where

[TABLE]

and $I_{{2L_{x}+2}}$ is an identity matrix of $(2L_{x}+2)$-order, same for the zero matrix $0_{2L_{x}+2}$.

Appendix D Effective-field theory

In order develop nonzero temperature phase diagrams, for comparison with the LSWT, we use the Effective Field Theory 1 (EFT-1) approach. This treatment consider clusters with a single central spin. The EFT-1 Hamiltonian for the present honeycomb lattice with coordination number $z$, is given by

[TABLE]

where, $p$ denotes the sublattices A and B. While, $p^{\prime}$ stands for an opposite sublattice respect to $p$.

From the Eq. 21 and using the approximate Callen-Suzuki relation barreto (50), we obtain the magnetizations $m_{p}=\langle S_{1p}^{z}\rangle$

[TABLE]

for $a_{1p}=J\overset{z}{\sum\limits_{\delta}}S_{(1+\delta)p^{\prime}}^{z}$.

Applying the identity $\exp(\alpha D_{x})F(x)=F(x+a)$ (with $D_{x}=\frac{\partial}{\partial x}$ is the differential operator) and the Van der Waerden identity for the spin up and down (i.e., $\exp(aS_{i}^{z})=\cosh(a)+S_{i}^{z}\sinh(a)$), the Eq. 22 leads

[TABLE]

with

[TABLE]

for $\alpha_{x}=\cosh(JD_{x})$ and $\beta_{x}=\sinh(JD_{x})$. The Eq. 23 is expressed in terms of multiple spin correlation functions. It is practically impossible to obtain all the spin correlation functions. For this reason, a decoupling of the right-hand side in Eq. 23) is needed, namely

[TABLE]

where $i\neq j\neq\ldots\neq l$ and $m_{p}=\left\langle S_{ip}^{z}\right\rangle$. Is worth to note, the approximation of Eq. 25 neglects correlations between different spins. However, relations as $\left\langle\left(S_{ip}^{z}\right)^{2}\right\rangle=1$ are exact, while in the usual MFA all the self- and multispin correlations are neglected. After this, the Eq. 23 can be written as

[TABLE]

with

[TABLE]

The coefficients $A_{q}(T_{N},h_{z},h_{x})$ are obtained by using the relation $\exp(\alpha D_{x})F(x)=F(x+a)$.

We can define the uniform and staggered magnetizations as $m=\frac{1}{2}(m_{A}+m_{B})$ and $m_{s}=\frac{1}{2}(m_{A}-m_{B})$, respectively. Accordingly, near the critical point we have $m_{s}\rightarrow 0$ and $m\rightarrow m_{0}$. So, the sublattice magnetization $m_{A}$ expanded up to linear order in $m_{s}$ (order parameter) is given by

[TABLE]

with

[TABLE]

and

[TABLE]

The present technique works only for second-order transitions. Therefore, to study the phase diagram we just analyze the Eq. 29) in the limit of $m_{\mathbf{s}}\rightarrow 0$. Thus, once the second-order line is located and using the fact that $m_{A}=m_{0}+m_{s}$ in Eq. 28), at the critical point we get

[TABLE]

[TABLE]

in which $m_{s}=0$, $\widetilde{h}_{x}=h_{x}/J$ and $\widetilde{h}_{z}=h_{z}/J$.

We can now obtain an analytical solution for the second-order transition, where $m_{s}$ is the order parameter, to describe the phase transition of the present model Eq. 1. The phases are identified following $m_{s}\neq 0$ (different sublattice magnetizations) for the antiferromagnetic phase and $m_{s}=0$ (equal sublattice magnetizations) for the polirized phase.