\documentclass[amsmath,amssymb,showpacs,twocolumn,dvips,epsfig,aps,prb]{revtex4-2}
\usepackage{graphicx}
\usepackage{dcolumn}
\usepackage{bm}
%\usepackage{hyperref}
%
\begin{document}
%\title{Spin transport theory of non-Hermitian quantum spin systems}
\title{Non-Hermitian linear-response theory for spin diffusion in quantum systems}
\author{Leonardo S. Lima}
\affiliation{Department of Physics, Federal Technological Education Center of Minas Gerais, 30510-000, Belo Horizonte, MG, Brazil}%, e-mail: lslima@cefetmg.br, Tel. +5531998996706.}

\date{\today}
\begin{abstract}
Spin transport theory of non-Hermitian quantum systems, where the non-Hermiticity  led to a series of exotic phenomena when the system experiences dissipation to an environment is proposed. The
goal is a better understanding of the spin diffusion in the one-dimensional non-Hermitian XXZ model which allows the transport coefficients be computed analytically. Moreover, we analyze the electric transport of the  one-dimensional non-Hermitian Hubbard model which is a very important model of electrons strongly correlated, where we have investigated the effect of non-Hermitian parameters like the imaginary hopping  on AC and DC conductivities of the system. We analyze the large $U$ limit of this model, where such non-Hermiticity contributes with a minus sign in the virtual exchange of quasiparticles and the behavior of the ground state energy and low-lying excitations is reversed.
\end{abstract}

\keywords{transport, non-Hermitian, XXZ model, Hubbard model.}
\pacs{75.10.Pq, 75.40.Mg}
\maketitle
%
\section{Introduction}{\label{QMC}}

It is well known that the non-Hermitian quantum mechanics can provide a description  of dissipative systems in a minimalist fashion\cite{N,XZZhang,kevin,Lei}. It has gained much attention recently due to the rapid progresses in experimental implementations of non-Hermiticity as photonic experiments \cite{CPoli,JMZeuner,LXiao,XZhan,SWei,MPra,HZhou,MABan} and ultra-cold atoms\cite{JLi}, where has been revealed that the non-Hermiticity changes the properties of a large number of well-known quantum phenomena well explained in the Hermitian physics, ranging from quantum and topological phase transitions, to quantum magnetism\cite{CM1,CM2,MN,YA,KK,MS,TE,HS,FK,SY,ZG,KK2,KK3,XZ,TE}.
The purpose of the present paper is to study the spin transport theory and the electric transport theory for the  one-dimensional non-Hermitian XXZ model and the one-dimensional non-Hermitian Hubbard model respectively, using the non-Hermitian linear-response theory. The effect of splitting magnon bands, induced by the non-Hermitian parameters in properties as quantum correlation has been analyzed by us in the Refs.\cite{lslima_Physe,lslima20212,lslima20214}.
 The paper is organized as follows.  In section \ref{4}, we discuss about the non-Hermitian spin transport theory of quantum spin systems. In section \ref{5}, we  analyze the case of the  half-integer spin and integer spin one-dimensional non-Hermitian XXZ model and the very important case of strongly correlated electrons systems modelled by the one-dimensional non-Hermitian Hubbard model, where we present the analytical and numerical results for the transport coefficients.  In the last section \ref{7}, we present our conclusions and final remarks.

\section{Non-Hermitian linear-response theory for spin currents}{\label{4}}

In the linear response theory for non-Hermitian systems with Hamiltonian $\mathcal{H}=\mathcal{H}_0+\mathcal{H}_{dis}$, where the non-Hermitian dissipation term corresponding to a It\^o's stochastic differential equation\\ $dx(t)=\gamma a(x(t),t) dt+b(x(t),t)\circ dW(t)$, with the Wiener increment $dW(t)=\xi(t)dt$ and white noise term $\xi(t)$
\begin{equation}{\label{nonher}}
\mathcal{H}_{dis}=\sum_{j}\left(-i\gamma\mathcal{B}^{\dag}_j\mathcal{B}+\mathcal{B}_j^{\dag}\xi_j+\xi_j^{\dag}\mathcal{B}_j\right)
\end{equation}
where $\Gamma=2\gamma$, $\gamma\in\mathbb{R}$ and $\xi(t)$ is a stochastic noise that satisfies $\langle\xi_i(t)\xi^{\dag}_j(t')\rangle=\Gamma\delta_{ij}\delta(t-t')$, $\mathcal{B}_j$ is a Hermitian operator. The linear response to Eq.~(\ref{nonher}) is given by\cite{Lei}
\begin{equation}
\delta\mathcal{H}(t)=-\Gamma\sum_{j}\int_{0}^{t}\langle\{\mathcal{W}(t),\mathcal{B}^{\dag}_j(t')\mathcal{B}_j(t')\}-2\mathcal{B}^{\dag}_j(t')\mathcal{W}(t)\mathcal{B}_j(t')\rangle dt',
\end{equation}
where the physical observable $\mathcal{W}(t)$ is given by
\begin{eqnarray}
\mathcal{W}=\langle\hbox{Tr}\left(\rho\mathcal{W}_{\mathcal{H}}(t)\right)\rangle_{\hbox{noise}},\\
\rho=\frac{e^{-\mathcal{H}_0/T}}{\hbox{Tr}\left\{e^{-\mathcal{H}_0/T}\right\}}
\end{eqnarray}
 and the operator $\mathcal{W}(t)$ in the Heisenberg picture is given by $\mathcal{W}(t)=e^{i\mathcal{H}^{\dag}t}\mathcal{W}e^{-i\mathcal{H}^{\dag}t}$.
Thus, the response of the system to the frequency-dependent gradient of the external magnetic field $\mathbf{H}$ generates a spin current given by $\mathcal{J}=\sigma\nabla\mathbf{H}$, where the response linear to the external field in $x$ direction is
\begin{equation}
\langle\mathcal{J}_x(l,t)\rangle=\sum_{j}\int_{-\infty}^{\infty}dt'\chi_{jS}(l,j,t-t')h^z(j,t'),
\end{equation}
being the response function defined as
\begin{equation}
\chi_{jS}(l,j,t-t')\equiv i\Theta(t-t')\langle 0|[\mathcal{J}_x(l,t),S^z(j,t')]|0\rangle,
\end{equation}
where $\Theta$ is the Heaviside step function.
On the other hand, the non-Hermitian response function is given by\cite{kevin}
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  \chi_{jS}^{NH}(l,j,t-t')=-\frac{1}{\hbar}\Theta(t-t')\bigg[\langle 0|\{\mathcal{J}_x(l,t),S^z(j,t')\}|0\rangle\nonumber\\-2\langle 0|\mathcal{J}_x(l,t)|0\rangle\langle 0|S^z(j,t')|0\rangle\bigg],\nonumber\\
\end{eqnarray}
where $\{\cdot\cdot\cdot\}$ is the unequal-time anti-commutator to establish the link between the response function and the correlation function. We have the non-Hermitian dynamic susceptibility as the Fourier transform
\begin{equation}
\chi_{jS}^{NH}(\tau,\omega)=\int_{-2\tau}^{2\tau}d\Delta t\chi_{jS}^{NH}\left(\tau+\frac{\Delta t}{2},\tau-\frac{\Delta t}{2}\right)e^{i\omega\Delta t},
\end{equation}
where $\tau=it$. $\chi_{jS}$ is the non-Hermitian response function
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  \chi_{jS}^{NH}(t,t')=\frac{i}{\hbar} \Theta(t-t')\langle 0|\{\mathcal{J}_x(t),S^z(t')\}|0\rangle.
\end{eqnarray}

The wave-vector-dependent susceptibility is given by
\begin{equation}\label{su}
\chi_{jS}^{NH}({k},\omega)\equiv\frac{i}{N}\int_{0}^{\infty}dte^{i(\omega+i0^{+})t} \langle 0|\{\mathcal{J}_x({k},t),S^z(-{k},t')\}|0\rangle.
\end{equation}
From continuity equation for the spin current:\\ $\dot{S}^z({k},t)+i{k}\cdot\mathcal{J}_x({k},t)=0$, $\chi_{jS}^{NH}$ can be transformed as follows:
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  \chi_{jS}^{NH}({k},\omega)=\frac{i}{N}\frac{1}{i(\omega+i0^{+})}\bigg[ik_x\int_{0}^{\infty}dte^{i(\omega+i0^{+})t}\nonumber\\
  \times\langle 0|\{\mathcal{J}_x({k},t),\mathcal{J}_x(-{k},0)\}|0\rangle%\nonumber\\
  -\langle 0|\{\mathcal{J}_x({k},0),S^z(-{k},0)\}|0\rangle\bigg].\nonumber\\
\end{eqnarray}
Using the representation of the spin current operator in terms of spin operators
\begin{eqnarray}\label{current}
\mathcal{J}_x(j)=\frac{iJ}{2}\sum_x(S_{j}^{+}S_{j+x}^{-}-S_{j}^{-}S_{j+x}^{+}),
\end{eqnarray}
where $j+x$ is the nearest-neighbor site of the site $j$ in the positive $x$ direction, we can transform the second term as
\begin{eqnarray}
\langle 0|\{\mathcal{J}_x({k},0),S^z(-{k},0)\}|0\rangle=\frac{i}{2}\sum_{l,x}\mathcal{J}_{l,l+x}\left(e^{ik_x}-1\right)\nonumber\\
\times\langle S_{j}^{+}S_{j+x}^{-}+S_{j}^{-}S_{j+x}^{+}\rangle.\nonumber\\
\end{eqnarray}
In the long-wavelength $k_x\rightarrow 0$ limit the susceptibility $\chi_{jS}^{NH}({k},\omega)$ is thus proportional to $ik_x$ and we can write
\begin{equation}
\langle\mathcal{J}_{x}({k},\omega)\rangle=-\frac{\langle-\mathfrak{K}_x\rangle-\mathfrak{F}({k},\omega)}{i(\omega+i0^{+})}ik_xh^z({k},\omega),
\end{equation}
where $\langle-\mathfrak{K}_x\rangle$ is the kinetic energy % being given by
\begin{eqnarray}
\langle \mathfrak{K}_x \rangle=-\frac{J}{N}\sum_{j}\left(S_j^{+}S_{j+x}^{-}+S_j^{-}S_{j+x}^{+}\right),
\end{eqnarray}
and $\mathfrak{F}$ is the Green's function defined in  $T=0$ by\cite{Mahan}
\begin{equation}\label{Green}
\mathfrak{F}({{k}},\omega)=\frac{i}{ N}\int_{0}^{\infty}dte^{i\omega t}\langle 0| [\mathcal{J}_x({{k}},t),\mathcal{J}_x(-{{k}},0)]|0\rangle.
\end{equation}
%being ${\mathbf{T}}$, the time ordering operator.

The regular part of the conductivity $\sigma$ (continuum conductivity) in the context of Hermitian quantum mechanics is given by\cite{Kubo,Mahan,Sentef,Pires1,lslima2013,Limadz}:  $\hbox{Re}\left[\sigma(\omega)\right]=D_S(T)\delta(\omega)+\sigma^{reg}(\omega)$, where
\begin{eqnarray}\label{sigmareg}
\langle\mathcal{J}_{\alpha}({k},\omega)\rangle=\sum_{\beta}\sigma_{\alpha\beta}({k},\omega)ik_{\alpha} h_{\beta}({k},\omega),\nonumber\\
\sigma_{\alpha\beta}({k},\omega)=\hbox{Re}[\sigma_{\alpha\beta}({k},\omega)]+i\hbox{Im}[\sigma_{\alpha\beta}({k},\omega)]\nonumber\\
\sigma^{reg}(\omega)=\frac{\hbox{Im} \{\mathfrak{F}({{k}}=0,\omega)\}}{\omega}
\end{eqnarray}
and $\alpha,\beta=x,y,z$. $D_S(T)$ is the spin Drude's weight. The effective $T$ that best relates the susceptibilities via fluctuation dissipation relation for a fixed waiting time $t_w$ is given as\cite{kevin}
\begin{eqnarray}{\label{temper}}
\hspace{-0.0cm}T=\arg\min_{\Theta}\int d\omega\left[-\chi'^{NH}(t_w,\omega)\tanh\left(\frac{\hbar\omega}{2k_B\Theta}\right)-\chi"(t_w,\omega)\right]^2,\nonumber\\
\end{eqnarray}
where $\chi^{NH}=\chi'^{NH}+i\chi"^{NH}$, $\chi=\chi'+i\chi"$. % For $\delta=0$ we have the Hermitian model and $\delta\ne0$ the model is non-Hermitian.  %We obtain a small difference in the behavior of the curves for the two models (Hermitian and non-Hermitian) due to transformation of the non-Hermitian Hamiltonian  in  Hermitian,  Eq.~(\ref{model}). Moreover, for $T$ non-zero, $D_S(T)$ rises with $T$ however, this description is only qualitative due to approach used.


\begin{figure}
    \centering
%\includegraphics[width=6.0cm]{spintransp_nonherm_XXZ} \\
\includegraphics[width=6.0cm]{spintransp_nonherm_XXZ_corr}\\
\caption{$\sigma^{reg}(\omega)$ at $T=0.0J$ for different values of non-Hermitian coupling $\delta$ for the spin-1/2 non-Hermitian XXZ model. For $\delta=0$ we have that the model is Hermitian and $\delta\ne0$ the model is non-Hermitian. We find that conductivity tends to zero at DC limit.}\label{fig_33}
\end{figure}

\section{Results and discussion}{\label{5}}

As example of application of the method, we consider the non-Hermitian XXZ model given by the Hamiltonian
\begin{equation}
\mathcal{H}=-\sum_{i,j\ne i}\left(S_i^{+}S_{j}^{-}+S_i^{-}S_{j}^{+}\right)+\Delta\sum_{i,j\ne i} S_{i}^{z}S_{j}^z+\sum_{j}\mathbf{H}\cdot S_j,
\end{equation}
where the non-Hermiticity originates from complex magnetic field $\mathbf{H}=\left(1,-i\delta,0\right)$, which can be understood as the spin-dependent losses and are within the reach of ultracold atom experiments.
The degeneracy of the system breaks down when the external complex field presents the form $\mathfrak{T}\mathbf{H}\cdot S_j\mathfrak{T}^{-1}=-\mathbf{H}^{*}\cdot S_j$. The external field also spoils the commutation relation $\left[\sum_{j}S_j^{z},\mathcal{H}\right]=0$.

We can recast the Hamiltonian above as
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
 \tilde{\mathcal{H}}=-\sum_{i,j\ne i}\left(\tau_i^{+}\tau_{j}^{-}+\tau_i^{-}\tau_{j}^{+}+\Delta \tau_{i}^{z}\tau_{j}^z\right)+\sum_{i}\sqrt{1-\delta^2}\tau^x_j,\nonumber\\
\end{eqnarray}
which is the XXZ model in a transverse field. The transformation hold if and only if $|\delta|\ne 1$.
Considering a local complex field $\mathbf{H}\cdot S_N$, the Hamiltonian shares two eigenstates $|\psi\rangle=\pm\sqrt{1-\delta}|\uparrow\rangle+\sqrt{1+\delta}|\downarrow\rangle$, which satisfies $\mathcal{H}|\psi\rangle=(-N\Delta/4\pm\sqrt{1-\delta^2}|\psi\rangle$, where $|\psi\rangle$ coalesce at $|\delta|=1$ where the ground state is dramatically changed from degeneracy to coalescence by the local complex field.  Under a homogeneous global complex field, the Hamiltonian of the free spins in the absence of the interaction describes the $\mathcal{PT}$-symmetric hypercube graph of $N$ dimension, and the system can be projected onto several invariant subspaces denoted by $S$. $\delta=1$ us the $n=2S+1$ of $n$-eigenstate coalescence in each subspace.


\textit{Half-integer spin S}:
Performing the Jordan-Wigner transformation
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  \tau_j^x =\frac{1}{2}-\bar{f}_jf_j, \nonumber\\
  \tau_j^y =\frac{i}{2}\sum_{j<l}(1-2\bar{f}_jf_j)(\bar{f}_j-f_j), \nonumber\\
  \tau_j^{z}=-\frac{1}{2}\sum_{j<l}(1-2\bar{f}_jf_j)(\bar{f}_j+f_j)
\end{eqnarray}
to replace the quasi-spin operators by the new non-Hermitian operators  $f_j$ and $\bar{f}_j$, where $\bar{f}_j=\mathfrak{A}_jc_j^{\dag}\mathfrak{A}_j^{-1}$, $f_j=\mathfrak{A}_jc_j\mathfrak{A}_j^{-1}$, $c_j^{\dag}$, $c_j$ are the creation and annihilation operators of spinless fermions. The new operators satisfy the fermionic anticommutation relation $\{\bar{f}_j,f_j'\}=\delta_{j,j'}$. The parity of the number of fermions is a conservative quantity such that the Hamiltonian can be expressed as
$\tilde{\mathcal{H}}=\tilde{\mathcal{H}}_{+}\mathbf{I}=\tilde{\mathcal{H}}_{-}\mathbf{I}$, where $\tilde{\mathcal{H}}_{+}=\tilde{\mathcal{H}}_{-}=-2(\bar{f}_{N}\bar{f}_1+\bar{f}_{N}f_1+\bar{f}_1f_{N}+f_1f_{N})$ and the Hamiltonian can be rewritten as
\begin{eqnarray}
\tilde{\mathcal{H}}=\frac{1}{4}\sum_{j=1}^{N}[2g\sqrt{1-\delta^2}(1-2\bar{f}_jf_j)\nonumber\\
+\Delta(\bar{f}_j\bar{f}_{j+1}+\bar{f}_jf_{j+1}+\bar{f}_{j+1}f_{j}+f_{j+1}f_{j})\nonumber\\
-(\bar{f}_j\bar{f}_{j+1}-\bar{f}_jf_{j+1}-\bar{f}_{j+1}f_{j}+f_{j+1}f_{j})\nonumber\\
+(1-2\bar{f}_{j+1}f_{j+1}-2\bar{f}_{j}f_j+4\bar{f}_{j+1}f_{j+1}\bar{f}_{j}f_j)].
\end{eqnarray}
%
\begin{figure}
    \centering
%\includegraphics[width=6cm]{nonhermitian_Drude_XXZ} \\
\includegraphics[width=6cm]{nonhermitian_Drude_XXZ_corr}\\
\caption{Behavior of the Drude's weight  $D_S(T)$ as a function of $T$ for the spin-1/2 non-Hermitian XXZ model. For $\delta=0$ we have that the model is Hermitian and $\delta\ne0$, the model is non-Hermitian. The Drude's weight is finite at $T=0$ indicating thus an ideal spin conductor at this limit.}\label{fig_46}
\end{figure}
Taking the Fourier transform
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  f_j=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}f_{\mathbf{k}}e^{i{\mathbf{k}\cdot\mathbf{r}_j}}, \hspace{0.5cm}\bar{f}_j=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}f_{\mathbf{k}}e^{-i{\mathbf{k}\cdot\mathbf{r}_j}},
\end{eqnarray}
\begin{eqnarray}
\tilde{\mathcal{H}}=\frac{1}{4}\sum_{j=1}^{N}[2g\sqrt{1-\delta^2}(1-2\bar{f}_{\mathbf{k}}f_{\mathbf{k}})\nonumber\\
+\Delta(e^{-i{\mathbf{k}}}\bar{f}_{\mathbf{k}}\bar{f}_{-{\mathbf{k}}}+2\cos({\mathbf{k}})\bar{f}_{\mathbf{k}}f_{{\mathbf{k}}}+e^{i{\mathbf{k}}}f_{{\mathbf{k}}}f_{-{\mathbf{k}}})\nonumber\\
-(e^{-i{\mathbf{k}}}\bar{f}_{\mathbf{k}}\bar{f}_{-{\mathbf{k}}}-2\cos({\mathbf{k}})\bar{f}_{\mathbf{k}}f_{{\mathbf{k}}}+e^{i{\mathbf{k}}}f_{{\mathbf{k}}}f_{-{\mathbf{k}}})\nonumber\\
+(1-4\bar{f}_{{\mathbf{k}}}f_{{\mathbf{k}}}+4\bar{f}_{{\mathbf{k}}}f_{{\mathbf{k}}}\bar{f}_{{\mathbf{k}}}f_{\mathbf{k}})].
\end{eqnarray}
where $g=1$, $|\mathbf{k}|=k=2\pi(n+1/2)/N$, $n=0,1,2,...,N-1$, the Hamiltonian can be written as
\begin{equation}\label{ham}
  \tilde{\mathcal{H}}_{+}=\bar{\mathcal{H}}_{-}=\sum_{0<k<\pi}\bar{\chi}_{\mathbf{k}}\bar{\mathcal{H}}_{+}(\mathbf{k})\chi_{\mathbf{k}},
\end{equation}
with $\bar{\chi}_{\mathbf{k}}=(\bar{f}_{\mathbf{k}},f_{-\mathbf{k}})$, $\chi_{\mathbf{k}}=(f_{\mathbf{k}},\bar{f}_{-\mathbf{k}})^{T}$, and
\begin{eqnarray}
\hspace{-0.00cm} \bar{\mathcal{H}}_{+}({\mathbf{k}})=\left(
  \begin{array}{cc}
    (\Delta+1)\cos(k) -\lambda-1 & i(\Delta-1)\sin(k) \\
    -i(\Delta-1)\sin(k) & -((\Delta+1)\cos(k) -\lambda-1)  \\
  \end{array}
\right),\nonumber\\
\end{eqnarray}
where $\lambda=2g\sqrt{1-\delta^2}$ and we make $g=1$ to simplify the notation. In following, making the non-Hermitian Bogoliubov transformation
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  \bar{\psi}_{\mathbf{k}}=\cos\frac{\varphi_{\mathbf{k}}}{2}\bar{f}_{\mathbf{k}}+i\sin\frac{\varphi_{\mathbf{k}}}{2}f_{-\mathbf{k}}, \\
  \psi_{\mathbf{k}}=\cos\frac{\varphi_{\mathbf{k}}}{2}f_{\mathbf{k}}-i\sin\frac{\varphi_{\mathbf{k}}}{2}\bar{f}_{-\mathbf{k}},
\end{eqnarray}
where $[\bar{\psi}_{\mathbf{k}},\psi_{\mathbf{k'}}]=\delta_{\mathbf{k},\mathbf{k'}}$ and\\ $\varphi_{\mathbf{k}}=\tan^{-1}[(\Delta+1)\sin({k})/(\lambda+1-(\Delta+1)\cos({k}))]$, the Hamiltonian is recast in the diagonal form\\ $\tilde{\mathcal{H}}_{+}=\sum_{k}\left(\bar{\psi}_{\mathbf{k}}{\psi}_{\mathbf{k}}-1/2\right)$, with the dispersion relation of quasi-particles given by
\begin{equation}
\omega_{\mathbf{k}}=\sqrt{\left[(\Delta+1)\cos({k})-\lambda-1\right]^2-(\Delta-1)^2\sin^2({k})},
\end{equation}
If $|\delta|< 1$, the single-particle spectrum is real and $|\delta|> 1$, the system presents a complex single-particle spectrum regardless of $\mathbf{k}$.
%

The regular part of the conductivity $\sigma$ (continuum conductivity) in the context of Hermitian quantum mechanics is given by\cite{Kubo,Mahan}:  $\hbox{Re}\left[\sigma(\omega)\right]=D_S(T)\delta(\omega)+\sigma^{reg}(\omega)$, where
 $D_S(T)$ is the spin Drude's weight, being given by
\begin{equation}\label{drude}
\hspace{-0.0cm}D_S(T)\approx-\frac{1}{N}\sum_{{k}}\frac{\cos(k_{x})}{\omega_{\mathbf{k}}}[1+n(\omega_{\mathbf{k}})],
\end{equation}
where $n(\omega_{\mathbf{k}})=1/(e^{\beta\omega_{\mathbf{k}}}- 1)$ is the boson occupation number and $\beta=1/T$.

In Fig.~\ref{fig_33}, we present the behavior of $\sigma^{reg}(\omega)$ for different values of non-Hermitian coupling $\delta$. We obtain the AC conductivity tending to zero at $\omega\rightarrow 0$  however, as we have $\sigma(0)=D_S\delta(\omega)$ and since that we obtain $D_S$ finite, we must have a divergence in the DC current. However, the scattering among quasi-particles must introduce a spreading in the conductivity where in a real system the conductivity must to stay finite. The large peaks obtained for the conductivity are due to the behavior of the dispersion relation at range $(0.0<\omega/J<0.3)$, generating so, resonance effects on conductivity. In Figs.~\ref{fig_34} and \ref{fig_35}, we analyze the AC conductivity for the non-Hermitian XXZ model. In this case, we get the continuum conductivity tending to zero at DC limit, $\omega\rightarrow 0$.

The behavior of the Drude's weight  $D_S(T)$ as a function of $T$ is displayed in Fig.~\ref{fig_46} for the one-dimensional spin-1/2 non-Hermitian XXZ model.  Since $D_S(T)$ is finite for all $T$ values, we get that the system is an ideal spin conductor. For $\delta=0$ we have that the model is Hermitian and $\delta\ne0$ the model is non-Hermitian. Moreover, we get a significant difference among the behavior of the curves for the two models (Hermitian and non-Hermitian), due to transformation of the non-Hermitian Hamiltonian in Hermitian. For non-zero temperature, $D_S(T)$ rises with the temperature due to spinons thermally excited for $T$ non-zero.

\begin{figure}
    \centering
%\includegraphics[width=6.0cm]{spintransp_nonherm_XXZ} \\
\includegraphics[width=6.0cm]{spintransp_nonherm_XXZ_spin1}\\
\caption{$\sigma^{reg}(\omega)$ at $T=0.0J$ for different values of non-Hermitian coupling $\delta$ for the spin-1 non-Hermitian XXZ model. For $\delta=0$ we have the Hermitian model and $\delta\ne0$, the model is non-Hermitian. We find that conductivity tends to zero at DC limit.}\label{fig_43}
\end{figure}
\begin{figure}
    \centering
%\includegraphics[width=6cm]{nonhermitian_Drude_XXZ} \\
\includegraphics[width=6cm]{nonhermitian_Drude_XXZspin1_corr}\\
\caption{Behavior of the Drude's weight  $D_S(T)$ as a function of $T$ for the spin-1 non-Hermitian XXZ model. For $\delta=0$ we have the Hermitian model and $\delta\ne0$, the model is non-Hermitian. The Drude weight is finite at $T=0$ indicating thus an ideal spin conductor at this limit.}\label{fig_56}
\end{figure}

\textit{Integer spin S}:
 Considering the ferromagnet case, we perform the Holstein-Primakoff transformation linearized
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
\hspace{-1cm}  S^{y}_j\approx\frac{\sqrt{2S}}{2}(a^{\dag}_j+a_j),\hspace{0.3cm} S^{z}_j\approx\frac{\sqrt{2S}}{2i}(a^{\dag}_j-a_j),\hspace{0.3cm}S^{x}_j=S-a^{\dag}_ja_j\nonumber\\
\end{eqnarray}
and in following, the Fourier transformation
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  a_{\mathbf{k}}=\frac{1}{\sqrt{N}}\sum_{j=1}^{N} a_je^{i\mathbf{k}\cdot r_j},\hspace{0.5cm}a^{\dag}_{\mathbf{k}}=\frac{1}{\sqrt{N}}\sum_{j=1}^{N} a^{\dag}_je^{-i\mathbf{k}\cdot r_j}.\nonumber\\
\end{eqnarray}
We get
\begin{eqnarray}
%\hspace{-1.7cm} {\mathcal{H}}=\left(
%  \begin{array}{cc}
%    (\Delta+1)\cos(k)-\frac{S}{2}(\Delta-1)+S(S+\lambda) & i(\Delta+1)\sin(k) \\
%    -i(\Delta+1)\sin(k) & -((\Delta+1)\cos(k)-\frac{S}{2}(\Delta-1)+S(S+\lambda)  \\
%  \end{array}
%\right),\nonumber\\
\hspace{-0.25cm}{\mathcal{H}}=\left[(\Delta+1)\cos(k)-\frac{S}{2}(\Delta-1)+S(S+\lambda)\right]\sigma_3
-(\Delta+1)\sin(k)\sigma_2,\nonumber\\
%{\mathcal{H}}=\left[(\Delta+1)\cos(k)-\frac{S}{2}(\Delta-1)+S(S+\lambda)\right]\left(
%                                                                                 \begin{array}{cc}
%                                                                                   1 & 0 \\
%                                                                                   0 & -1 \\
%                                                                                 \end{array}
%                                                                               \right)\nonumber\\
%+(\Delta+1)\sin(k)\left(
%                    \begin{array}{cc}
%                      0 & i \\
%                      -i & 0 \\
%                    \end{array}
%                  \right),\nonumber\\
\end{eqnarray}
where $\sigma_1$, $\sigma_2$ and $\sigma_3$ are the (spin-1/2) Pauli's Matrices. We get the diagonal Hamiltonian
\begin{equation}
\mathcal{H}=\sum_{{k}}\varepsilon_{\mathbf{k}}a^{\dag}_{\mathbf{k}}a_{\mathbf{k}}
\end{equation}
with the dispersion relation of quasi-particles
\begin{equation}\label{quasi}
\hspace{-0.5cm}  \varepsilon_{\mathbf{k}}=\sqrt{\left((\Delta+1)\cos(k)-\frac{S}{2}(\Delta-1)+S(S+\lambda)\right)^2-(\Delta+1)^2\sin^2(k)}.
\end{equation}

In Figs.~\ref{fig_43}, we analyze the AC conductivity for the non-Hermitian spin-1 XXZ model. In this case, we get the continuum conductivity tending to zero at DC limit, $\omega\rightarrow 0$. We display the behavior of $\sigma^{reg}(\omega)$ for different values of non-Hermitian coupling $\delta$. As $\sigma(0)=D_S\delta(\omega)$ and since that we obtain $D_S$ finite for all $T$ values, we must have a divergence in the spin current. However, the scattering among magnons must introduce a spreading in the conductivity where in a real system the conductivity must to stay finite. The large peaks obtained for the conductivity are due to the behavior of the dispersion relation at range $0.0<\omega/J<0.5$ generating so, resonance effects on conductivity into this range of $\omega$.

The behavior of Drude's weight  $D_S(T)$ as a function of $T$ is displayed in Fig.~\ref{fig_56} for the one-dimensional spin-1 non-Hermitian XXZ model.  Since $D_S(T)$ stays finite  for all $T$ values, we get an ideal spin conductor behavior for the system for all $T$ values.  Moreover, we get a significant difference among the behavior of the curves for both models (Hermitian and non-Hermitian), due to transformation of the non-Hermitian Hamiltonian in Hermitian. For non-zero temperatures values, $D_S(T)$ rises with the temperature due to magnons thermally excited at this range of $T$. However this description for $T$ non-zero is only qualitative due to limitations of the spin wave approach used.

\textit{Non-Hermitian Hubbard model}: The model is given by the Hamiltonian
\begin{equation}{\label{model}}
\mathcal{H}=i\sum_{i,j}\sum_{\sigma=\uparrow,\downarrow}t_{ij}\left(c_{i,\sigma}^{\dag}c_{j,\sigma}+c_{j,\sigma}^{\dag}c_{i,\sigma}\right)+U\sum_{j}n_{j,\uparrow}n_{j,\downarrow},
\end{equation}
where $c_{i,\sigma}$ $(c_{j,\sigma}^{\dag})$ are the usual annihilation (creation) operators of fermions with spin $\sigma\in\{\uparrow,\downarrow\}$ at site $n_{i,\sigma}=c_{i,\sigma}^{\dag}c_{i,\sigma}$ is the number operator for a particle of spin $\sigma$ on site $l$. In addition, $U$, $t_{ij}\in\mathbb{R}$ are the interaction and kinetic energies respectively.
 The number of sites and particles are $\mathcal{N}$ and $\mathcal{M}$, respectively. The Hamiltonian above can describe a 1D homogeneous ring system in which $it_{ij}=it\delta_{i,j+1}$. Each subspace is labelled by ${k}=2n\pi/N$ which is the momentum indexing the subspace. Each subspace labeled by ${k}$ can be decomposed into four subspaces with $(S, S^z)=(0,0), (1,0)$ and $(1,\pm 1)$ in term of the spin symmetry. In each invariant subspace we have the two-particle solution, where the two-particles state is calculated in Ref.~\cite{XZZhang}. Based on the symmetry of the system, we have the paring mechanism through the exact solution with the two-particle subspace. The bound pair emerges in the $(0,0)$ subspace with energy of the bound pair given by $\xi_{{\mathbf{k}}}=\pm\sqrt{U^2-16t^2\cos^2({k}/2)}$. The details about the bound pair are given in Ref.~\cite{XZZhang}. For the case of $U$ negative, the lowest energy of a bound pair is $\omega_{\pi}=-U$ locating on the subspace with ${k}=\pi$. Furthermore, the bound pair state is $|\varphi^{b}_{{\mathbf{k}}}\rangle=\sum_rg_{{\mathbf{k}}}|\phi_r^{-}({{k}})\rangle$ with $g_{{\mathbf{k}}}^{-}(j)=1/\sqrt{2}$ if $j=0$ and $e^{-\gamma j}$, if $j\ne 0$, where $\gamma=\ln\left[\left(-U\pm \sqrt{U^2+4\lambda_{{\mathbf{k}}}^2}\right)/4it\cos({k}/2)\right]$. In the absence of on-site interaction $(U=0)$, only the scattering eigenstate with imaginary eigenenergy is present and the system does not accommodate a bound pair state. A nonzero interaction $(U\ne 0)$ leads thus, the emergence of a bound pair. When $|U|>|4t|$, the system possesses the full real bound pair spectrum however, a small $U$ results in the appearance of the imaginary bound pair energy.
\begin{figure}
    \centering
\includegraphics[width=6.0cm]{non_Hermitian_Hub_cond.eps} \\
\caption{$\sigma^{reg}(\omega)$ at $T=0.0$ for different values of non-Hermitian couplings $U$, for the one-dimensional non-Hermitian Hubbard model. We find $\sigma^{reg}(\omega\rightarrow 0)$ diverging at DC limit indicating thus, a super-current behavior.}\label{fig_34}
\end{figure}
%

The continuum part of the spin conductivity $\sigma^{reg}(\omega)$, is defined in terms of the Green's function $\mathfrak{F}(k,\omega)$.
We obtain the spin current operator in terms of the operators Fourier transformed $c_{\mathbf{k}}^{\dag}$ and $c_{\mathbf{k}}$ given by
%
\begin{eqnarray}\label{Jcurrent}
\hspace{-1.0cm}\mathcal{J}_x(t)=\sum_{k}\mathcal{J}_x(k,t)= -\sum_{{k}}\frac{\cos(k_{x})}{\xi_{{\mathbf{k}}}}c_{{\mathbf{k}}}^{\dag}{\Lambda}_{{\mathbf{k}}}(t)c_{{\mathbf{k}}}.
\end{eqnarray}
The spin current response function $\mathfrak{F}({{k}},\omega)$ at non-zero $T$ is given by\cite{Mahan}
\begin{equation}\label{Green}
\mathfrak{F}({{k}},\omega)=\frac{i}{ N}\int_{0}^{\infty}dte^{i\omega t}\langle 0| [\mathcal{J}_x({{k}},t),\mathcal{J}_x(-{{k}},0)]|0\rangle,
\end{equation}
\noindent
where
$\mathfrak{F}({{k}}=0,\omega\rightarrow 0)$ is the susceptibility or retarded Green's function\cite{Mahan}.
The retarded Green's function or dynamical correlation function is obtained after performing an analytical  calculation, where we obtain the result
%\begin{eqnarray}
%\mathfrak{F}({{k}},\omega)=\sum_{{k},{k}'}\frac{\sin(k_{x}')}{\xi_{{k}'}}\frac{\sin(k_{x})}{\xi_{{k}}}{\Lambda}_{{k}}(\omega),\nonumber\\
%\end{eqnarray}
%being
\begin{eqnarray}
{\Lambda}_{{\mathbf{k}}}(\omega)=\frac{1}{\pi^2}\int_{0}^{2\pi}d\omega_1\mathfrak{F}_0({k},\omega_1)\tilde{\mathfrak{F}}_0({k},\omega-\omega_1),\nonumber\\
\end{eqnarray}
where $\mathfrak{F}_0$, $\tilde{\mathfrak{F}}_0$ are the bare propagators.
\begin{eqnarray}
  \mathfrak{F}_0({k},\omega)= \frac{1}{\omega-\xi_{{\mathbf{k}}}+i0^{+}},\hspace{0.5cm}\tilde{\mathfrak{F}}_0({k},\omega)= -\frac{1}{\omega+\xi_{{\mathbf{k}}}-i0^{+}}.\nonumber\\
\end{eqnarray}
${\Lambda}_{{\mathbf{k}}}(\omega)$ is the Fourier transform of  ${\Lambda}_{{\mathbf{k}}}(t)$, which is the dynamical correlation function.
%\begin{eqnarray}
%\hspace{-0.5cm}{\Lambda}_{{k}}(t)-i\langle 0|\bigg[\psi_{{k}}(t){\bf{\sigma}_{\alpha\beta}}\psi^{\dag}_{{k}}(t)\psi_{{k}}(0){\bf{\sigma}_{\alpha\beta}}\psi^{\dag}_{{k}}(0)\nonumber\\
%\hspace{-0.5cm}+\psi_{{k}}(t){\bf{\sigma}_{\alpha\beta}}\psi^{\dag}_{{k}}(t)\psi_{{k}}(0){\bf{\sigma}_{\alpha\beta}}\psi^{\dag}_{{k}}(0)+\psi_{{k}}(t){\bf{\sigma}_{\alpha\beta}}\psi^{\dag}_{{k}}(t)\psi_{{k}}(0){\bf{\sigma}_{\alpha\beta}}\psi^{\dag}_{{k}}(0)\nonumber\\
%\hspace{-0.5cm}+\psi_{{k}}(t){\bf{\sigma}_{\alpha\beta}}\psi^{\dag}_{{k}}(t)\psi_{{k}}(0){\bf{\sigma}_{\alpha\beta}}\psi^{\dag}_{{k}}(0)\bigg]|0\rangle.\nonumber\\
%\end{eqnarray}
Thus, we obtain the regular part of the longitudinal spin conductivity $\sigma^{reg}(\omega)$ as being  given by
\begin{eqnarray}\label{sigma}
\hspace{-0.0cm}\sigma^{reg}(\omega)=\sum_{{k}}\frac{\cos^2(k_x)}{\xi_{{\mathbf{k}}}} n(\xi_{{\mathbf{k}}})\left(1-e^{\xi_{{\mathbf{k}}}/T}\right)\delta(\omega-2\xi_{{\mathbf{k}}}),\nonumber\\
\hspace{-1.0cm}\sigma^{reg}(\omega)=\sum_{{k}}\frac{\cos^2(k_x)}{\xi_{{\mathbf{k}}}}\tanh\left(\frac{\xi_{{\mathbf{k}}}}{2T}\right)\delta(\omega-2\xi_{{\mathbf{k}}}).\nonumber\\
%\hspace{-0.8cm}\sigma^{reg}(\omega)=\sum_{{k}}\sin^2(k_{x}) \frac{\left[1+2n(\xi_{{k}})\right]}{\xi_{{k}}^3}\delta(\omega-\xi_{{k}}).
\end{eqnarray}
The Drude's weight is given by
\begin{eqnarray}\label{drude}
\hspace{-0.0cm}D_S(T)\approx-\frac{1}{N}\sum_{k}\sin(k_x)\tanh\left(\frac{\xi_{{\mathbf{k}}}}{2T}\right),
%\hspace{-0.0cm}D_S(T)\approx-\frac{1}{N}\sum_{k}\frac{\cos(k_{x})}{\xi_{{k}}}\left(n(\xi_{{k}})+1\right),
\end{eqnarray}
where $n(\xi_{{\mathbf{k}}})=1/(e^{\beta\xi_{{\mathbf{k}}}}+ 1)$ is the fermions occupation number.
In Fig.~\ref{fig_35} and \ref{fig_34}, we present the behavior of Drude's weight and $\sigma^{reg}(\omega)$ for different values of non-Hermitian coupling $\delta$. We obtain the AC conductivity tending to zero at $\omega\rightarrow 0$  however, as we have $\sigma(0)=D_S\delta(\omega)$ and since that we obtain a $D_S$ finite, we must have a divergence for the DC current. However, the scattering among particles must introduce a spreading in the conductivity where in a real system the conductivity must to stay finite. In Figs.~\ref{fig_34} and \ref{fig_35}, we analyze the conductivity for the non-Hermitian model Eq.~(\ref{model}). In this case, we get a divergence in the continuum conductivity at DC limit, $\omega\rightarrow 0$. The behavior obtained for the AC conductivity is due to the form of the Eq.~(\ref{sigma}), which is a very complicated expression of ${k}$, involving thus, many processes that depends on ${k}$. Furthermore, as we obtain a finite Drude's weight for all values of $T$, we have a Dirac's delta peak for the conductivity at $\omega=0$ and consequently, we obtain that the transport is ideal  in this point ($\omega=0$) for all values of $T$. For  values nonzero of $\omega$ ($\omega\neq 0$), we obtain a decreasing in the conductivity for higher values of $T$ and $\omega$, although this behavior is only qualitative  due to approach used.
\begin{figure}
    \centering
\includegraphics[width=6.0cm]{D_nonhermit_Hub.eps} \\
\caption{Drude's weight  $D_S(T)$ as a function of $T$ for different values of non-Hermitian couplings $U$, for the one-dimensional non-Hermitian Hubbard model.  The Drude's weight is finite for all $T/J$ indicating thus, an ideal conductor for all $T$ values.}\label{fig_35}
\end{figure}

In all cases analyzed, the influence of dispersionless flat modes on  longitudinal spin conductivity is only to give rise to a Dirac's delta-like peak at  frequency $\omega=\xi_{{\mathbf{k}}}$, where $\xi_{{\mathbf{k}}}$ is a plane mode in each case. Furthermore, the presence of large peaks in the AC spin conductivity and a finite Drude's weight  $D_S(T)$, indicate a supercurrent behavior for the system although, for one has a superconductor behavior  is  necessary that the system  exhibits the  Meissner effect as well\cite{scientific}.

\section{Outlook}{\label{7}}
In summary, we present the spin transport theory of the  spin-1/2 and spin-1 one-dimensional non-Hermitian  XXZ model and the electrical transport theory of strongly correlated electrons systems modelled by the one-dimensional non-Hermitian Hubbard model. We get the variation of the transport coefficients with the non-Hermitian coupling parameters of the system.
  In a general way, in quantum spin systems either real fields or complex fields generate a splitting of degenerate ground states, where the spins are aligned along the direction of the external magnetic field. However, the eigenvalues and eigenvectors of the system with real spectrum will not suffer many changes with the external magnetic field and the initial state display a oscillating behavior and periodic among all possible spin configurations.



\appendix{\bf Acknowledgment}
 This work was partially supported by National Council for Scientific and Technological Development (CNPq) Brazil.
%
\appendix{\bf Conflict of Interest}\\
The author declares none conflict of interest.\\
%
\appendix{\bf Author Contribution Statement }\\
Leonardo S. Lima is the sole author of this manuscript.\\
%
\appendix{\bf Data Availability Statement}\\
All data generated or analysed during this study are included in this paper.\\
%
\bibliography{Bibliography}
\begin{thebibliography}{10}
%
\bibitem{N} N. Moiseyev, non-Hermitian Quantum Mechanics, Cambridge University Press, Cambridge, UK (2011).
\bibitem{XZZhang} X. Z. Zhang, Z. Song, $\eta$-pairing ground states in the non-Hermitian Hubbard modelPhys, Rev. B 103, 235153 (2021).
\bibitem{kevin} Geier, Kevin T., Hauke, Philipp, From Non-Hermitian Linear Response to Dynamical Correlations and Fluctuation-Dissipation Relations in Quantum Many-Body Systems, PRX Quantum 3, 030308 (2022).
\bibitem{Lei} L. Pan, X. Chen, Y. Chen, H. Zhai, Non-Hermitian linear response theory, Nat. Phys. 16, 767 (2020).
\bibitem{CPoli} C. Poli, M. Bellec, U. Kuhl, F. Mortessagne e H. Schomerus, Selective enhancement of topologically induced interface states in a dielectric resonator chain, Nat. Comum. 6 , 6710 (2015).
\bibitem{JMZeuner} J.M. Zeuner, M.C. Rechtsman, Y. Plotnik, Y. Lumer, S. Nolte, M.S. Rudner, M. Segev e A. Szameit, Observation of a Topological Transition in the Bulk of a Non-Hermitian System, Phys. Rev. Lett. 115 , 040402 (2015).
\bibitem{LXiao} L. Xiao, X. Zhan, Z. H. Bian, K. K. Wang, X. Zhang, X. P. Wang, J. Li, K. Mochizuki, D. Kim e N. Kawakami et al. L. Xiao, X. Zhan, Z. H. Bian, K. K. Wang, X. Zhang, X. P. Wang, J. Li, K. Mochizuki, D. Kim e N. Kawakami, W. Yi, H. Obuse, B. C. Sanders, P. Xue, Nat. Phys.  13 , 1117 (2017), Nat. Phys.  13 , 1117 (2017).
\bibitem{XZhan} X. Zhan, L. Xiao, Z. Bian, K. Wang, X. Qiu, B. C. Sanders, W. Yi e P. Xue, Detecting Topological Invariants in Nonunitary Discrete-Time Quantum Walks, Phys. Rev. Lett. 119 , 130501 (2017).
\bibitem{SWei} S. Weimann, M. Kremer, Y. Plotnik, Y. Lumer, S. Nolte, K. G. Makris, M. Segev, M. Rechtsman e A. Szameit, Topologically protected bound states in photonic parity-time-symmetric crystals, Nat. Mater. 16 , 433 (2017).
\bibitem{MPra} M. Parto, S. Wittek, H. Hodaei, G. Harari, M. A. Bandres, J. Ren, M. C. Rechtsman, M. Segev, D. N. Christodoulides e M. Khajavikhan, Edge-Mode Lasing in 1D Topological Active Arrays, Phys. Rev. Lett. 120, 113901 (2018).
\bibitem{HZhou} H. Zhou, C. Peng, Y. Yoon, C. W. Hsu, K. A. Nelson, L. Fu, J. D. Joannopoulos, M. Solja\v{c}\'{c} and B. Zhen, Observation of bulk Fermi arc and polarization half charge from paired exceptional points, Science 359 , 1009 (2018).
\bibitem{MABan} M. A. Bandres, S. Wittek, G. Harari, M. Parto, J. Ren, M. Segev, D. N. Christodoulides e M. Khajavikhan, Topological insulator laser: Experiments, Science 359 , eaar4005 (2018).
\bibitem{JLi} J. Li, A. K. Harter, J. Liu, L. de Melo, Y. N. Joglekar e L. Luo, Observation of parity-time symmetry breaking transitions in a dissipative Floquet system of ultracold atoms, Nat. Comum. 10 , 855 (2019).
\bibitem{CM1} C. M. Bender and S. Boettcher, Real Spectra in Non-Hermitian Hamiltonians Having PT Symmetry, Phys. Rev. Lett. 80, 5243 (1998).
\bibitem{CM2} C. M. Bender, Making sense of non-Hermitian Hamiltonians, Rep. Prog. Phys. 70, 947 (2007).
\bibitem{MN} M. Nakagawa, N. Kawakami, and M. Ueda, Non-Hermitian Kondo Effect in Ultracold Alkaline-Earth Atoms, Phys. Rev. Lett. 121, 203001 (2018).
\bibitem{YA} Y. Ashida, S. Furukawa, and M. Ueda, Parity-time-symmetric quantum critical phenomena, Nat. Commun. 8, 15791 (2017).
\bibitem{KK} K. Kawabata, Y. Ashida, and M. Ueda, Information Retrieval and Criticality in Parity-Time-Symmetric Systems, Phys. Rev. Lett. 119, 190401 (2017).
\bibitem{MS} M. S. Rudner, L. S. Levitov, Topological Transition in a Non-Hermitian Quantum Walk, Phys. Rev. Lett. 102, 065703 (2009).
%\bibitem{TE} T. E. Lee, Phys. Rev. Lett. 116, 133903 (2016).
\bibitem{HS} H. Shen, B. Zhen,  L. Fu, Topological Band Theory for Non-Hermitian Hamiltonians, Phys. Rev. Lett. 120, 146402 (2018).
\bibitem{FK} F. K. Kunst, E. Edvardsson, J. C. Budich, and E. J. Bergholtz, Biorthogonal Bulk-Boundary Correspondence in Non-Hermitian Systems, Phys. Rev. Lett. 121, 026808 (2018).
\bibitem{SY} S. Yao, Z. Wang, Edge States and Topological Invariants of Non-Hermitian Systems, Phys. Rev. Lett. 121, 086803 (2018).
\bibitem{ZG} Z. Gong, Y. Ashida, K. Kawabata, K. Takasan, S. Higashikawa, M. Ueda, Topological Phases of Non-Hermitian Systems, Phys. Rev. X 8, 031079 (2018).
\bibitem{KK2} K. Kawabata, T. Bessho, and M. Sato, Classification of Exceptional Points and Non-Hermitian Topological Semimetals, Phys. Rev. Lett. 123, 066405 (2019).
\bibitem{KK3} K. Kawabata, S. Higashikawa, Z. Gong, Y. Ashida, and M. Ueda, Topological unification of time-reversal and particle-hole symmetries in non-Hermitian physics, Nat. Commun. 10, 297 (2019).
\bibitem{XZ} X. Zhang and J. Gong, Non-Hermitian Floquet topological phases: Exceptional points, coalescent edge modes, and the skin effect, Phys. Rev. B 101, 045415 (2020).
\bibitem{TE} T. E. Lee and C.-K. Chan, Heralded Magnetism in Non-Hermitian Atomic Systems, Phys. Rev. X 4, 041001 (2014).
%
%
\bibitem{lslima_Physe} L. S. Lima, Spin Nernst effect and quantum entanglement in two-dimensional antiferromagnets on checkerboard lattice, Physica E 128, 114580 (2021).
\bibitem{lslima20212} L. S. Lima, Quantum correlation and entanglement in the Heisenberg model with biquadratic interaction on square lattice, Eur. Phys. J. D 75, 28 (2021).% https://doi.org/10.1140/epjd/s10053-021-00044-4.
%\bibitem{lslima20213} L. S. Lima, Eur. Phys. J. Plus 136, (2021) 789.
\bibitem{lslima20214} L. S. Lima, Quantum Phase Transition and Quantum Correlation in the Two-dimensional Honeycomb-bilayer Lattice Antiferromagnet, J. Low Temp. Phys. 205, 112 (2021).
\bibitem{Kubo}
    Kubo, R., Toda, M., \& Hashitsume, N., Statistical Physics II, Springer-Verlag, New York, (1985).
\bibitem{Mahan}
    Mahan, G. D., Many Particles Physics, Plenum, New York, (1990).
\bibitem{Pires1} Pires,  A. S. T., Lima, L. S., Spin transport in antiferromagnets in one and two dimensions calculated using the Kubo formula, Phys. Rev. B 79, (2009) 064401 (2009).
\bibitem{lslima2013} Lima, Leonardo. S., Spin transport of the quantum integer spin S one-dimensional Heisenberg antiferromagnet coupled to phonons, Eur. Phys. J. B 86,  99 (2013).
%\bibitem{lslima20201} Leonardo. S. Lima, J. Magn. Magn. Mater. 505 (2020) 166751.
\bibitem{Sentef}  Sentef, M., Kollar, M., Kampf, A. P., Spin transport in Heisenberg antiferromagnets in two and three dimensions, Phys. Rev. B 75, 214403 (2007).
\bibitem{Limadz} Lima, L. S., Low-temperature spin transport in the $S=1$ one- and two-dimensional antiferromagnets with Dzyaloshinskii–Moriya interaction, Phys. Status Solidi B, 249, 1613 (2012).
%
\bibitem{scientific} Lima, L. S., Antiferromagnetic and ferromagnetic spintronics and the role of in-chain and inter-chain interaction on spin transport in the Heisenberg ferromagnet, Scientific Reports 11, 20442 (2021).
\bibitem{Nija} Nita, M., Ostahie, B., Aldea, A., Spectral and transport properties of the two-dimensional Lieb lattice, Phys. Rev. B 87, 125428 (2013).
\bibitem{Rui} Mao, R., Dai, Yan-Wei, Cho, Sam Young, Zhou, Huan-Qiang, Quantum coherence and spin nematic to nematic quantum phase transitions in biquadratic spin-1 and spin-2 XY chains with rhombic single-ion anisotropy, Phys. Rev. B 103,  014446 (2021).
\bibitem{Cao} Cao, X., Chen, K., He, D., Magnon Hall effect on the Lieb lattice, J. Phys. Condens. Matter 27, 166003 (2015).
\bibitem{Micheli} Micheli, A., Brennen, G. K., Zoller, P., A toolbox for lattice-spin models with polar molecules, Nat. Phys. 2, 341 (2006).
\bibitem{JT} Chalker, J. T., Holdsworth, P. C. W., Shender, E. F.,  Hidden order in a frustrated system: Properties of the Heisenberg Kagomé antiferromagnet, Phys. Rev. Lett. 68, 855 (1992).
\bibitem{DA} Huse, D. A., Rutenberg, A. D., Classical antiferromagnets on the Kagom\'e lattice Phys. Rev. B 45, 7536(R) (1992).
\bibitem{Meghadeepa} Adhikary, M., Ralko, A., Kumar, B., Quantum paramagnetism and magnetization plateaus in a kagome-honeycomb Heisenberg antiferromagnet, Phys. Rev. B 104,  094416 (2021).
\end{thebibliography}


\end{document}































\begin{eqnarray}\label{sigmareg}
\langle\mathcal{J}_{\alpha}(\mathbf{k},\omega)\rangle=\sum_{\beta}\sigma_{\alpha\beta}(\mathbf{k},\omega)ik_{\alpha} h_{\beta}(\mathbf{k},\omega),\nonumber\\
\sigma_{\alpha\beta}(\mathbf{k},\omega)=\hbox{Re}[\sigma_{\alpha\beta}(\mathbf{k},\omega)]+i\hbox{Im}[\sigma_{\alpha\beta}(\mathbf{k},\omega)]\nonumber\\
%\hbox{Re}\left[\sigma_{\alpha\beta}(\omega)\right]=D_S(T)\delta(\omega)+\sigma^{reg}(\omega)\nonumber\\
\sigma^{reg}(\omega)=\frac{\hbox{Im} \{\mathfrak{G}({\mathbf{k}}=0,\omega)\}}{\omega}
\end{eqnarray}
and $\alpha,\beta=x,y,z$.
%\bibitem{Qi} Qi-Hui Chen, F. J. Huang, Young-Ping Fu, Phys. Rev. B 105, (2022) 224401.
%\bibitem{colpa} J. H. P. Colpa, Physica A 93, (1978) 327.
\bibitem{Ballentine} L. E. Ballentine, Am. J. Phys. 55, (1986) 785.
\bibitem{Divincenzo} D. P. DiVincenzo, Science 270, (1995) 255.
\bibitem{Fabio} F. Benatti, M. Fannes, R. Floreanini, D. Petritis, {\it{Quantum Information, Computation and Cryptography An Introductory Survey of Theory}},  Springer Heidelberg, Germany (2010).
\bibitem{Dan} D. C. Marinescu, G. M. Marinescu, {\it{Approaching quantum computing}}, Pearson Prentice Hall, New Jersey, USA (2004).
\bibitem{Michael} M. A. Nielsen, I. L. Chuang,  {\it{Quantum Computing and Quantum Information}}, Cambridge University Press, Cambridge, UK (2000).
\bibitem{Latorre3} J. I. Latorre, A. Riera, J. Phys. A.:Math. Theor. 42, (2009) 404002.
%
\bibitem{A} A. R. Its, B-Q Jin, V. E. Korepin, J. Phys. A: Math. Gen. 38, (2005) 2975.
\bibitem{Dagmar} D. Bruss, F. Leuchs,  {\it{Lectures on Quantum Information}}, WILEY-VCH Verlag, Weinheim, Germany (2007).
\bibitem{Latorre} J. I. Latorre, E. Rico, G. Vidal, Quant. Inf. Comput. 4, (2004) 48.
\bibitem{Calabrense} P. Calabrese, John Cardy J. Stat. Mech. (2004) P06002.
%
\bibitem{Davide} D. Bianchini, O. A. Castro-Alvaredo, B. Doyon, E. Levi, F. Ravanini, J. Phys. A: Math. Theor. 48 04FT01 (2014).
\bibitem{Vidal} G. Vidal, J. L. Latorre, E. I. Rico, A. Kitaev, Phys. Rev. Lett. 90, (2003) 227902.
\bibitem{Calabrense2} P. Calabrese, Physica A 504, (2018) 31.
\bibitem{Laflorence} N. Laflorencie, Physics Reports 646, (2016) 1.
%
%
\bibitem{K} K. Zyczkowski, P. Horodecki, A. Sanpera, M. Lewenstein, Phys. Rev. A 58, (1998) 883.
\bibitem{G} G. Vidal, R. F. Werner, Phys. Rev. A 65, (2002) 032314.
%
\bibitem{Plenio} M. B. Plenio, Phys. Rev. Lett. 95, (2005) 090503.
\bibitem{C} C. Castelnovo, Phys. Rev. A 88, (2013) 042319.
\bibitem{Y} Y. A. Lee, G. Vidal, Phys. Rev. A 88, (2013) 042318.
%\bibitem{Hill} S. A. Hill, W. K. Wootters, Phys. Rev. Lett. 78, (1997) 5022.
%\bibitem{W} W. K. Wootters, Phys. Rev. Lett. 80, (1998) 2245.


\textit{Concurrence}: Bipartite density matrices $\rho$ can be decomposed into infinitely many convex combinations of pure bipartite states
\begin{equation}{\label{conc}}
E[\rho]=\inf\left\{\sum_{j}\lambda_{j}\mathcal{S}(\rho_{1}^{j}):\rho=\sum_{j}\lambda_{j}|\psi^{j}\rangle\langle\psi^{j}|\right\},
\end{equation}
which is the smallest convex combination of the von Neumann entropies $\mathcal{S}(\rho_{1}^{j})$ of the reduced density matrices $\hbox{Tr}_2(|\psi^{j}\rangle\langle\psi^{j}|)$ of the pure states $|\psi^{j}\rangle\langle\psi^{j}|$ in terms of which $\rho=\sum_{j}|\psi^{j}\rangle\langle\psi^{j}|$, $\lambda_{j}\leq 0$, $\sum_{j}\lambda_{j}=1$, that is the entanglement of formation.
For a pair of qubits, the variational quantity Eq.~(\ref{conc}) can be substituted by the concurrence quantity\cite{Hill,W}
\begin{equation}
C=\max\{0,\lambda_1-\lambda_2-\lambda_3-\lambda_4\},
\end{equation}
where $\lambda_1\leq\lambda_2\leq\lambda_3\leq\lambda_4$  are the square roots of the eigenvalues of the operator $R=\rho(\sigma_y\otimes\sigma_y)\rho^{*}(\sigma_y\otimes\sigma_y)$ and $\lambda_1$ is the greatest square root of the four eigenvalues of $R$. This formula uses the "spin flip" transformation that is applicable to the states of an arbitrary number of qubits.\cite{Hill,W} For a state of two qubits $\rho$,  a nonzero concurrence means that the qubits 1 and 2 are correlated or entangled.  The concurrence $C=0$ corresponds to an unentangled state and $C=1$ corresponds to a maximally entangled state. For $\mathcal{N}$ qubits, we apply the above transformation to each individual qubit. \\
%{\it{Proposition}}: The concurrence of the nearest-neighboring qubits in the Heisenberg model is given by $C=1/2\max\left[0,-U/JN-1\right]$ for the antiferromagnet (AFM), $C=1/2\max\left[0,U/3JN-1\right]$ for the ferromagnet (FM).\cite{Vidal}
The relation between  concurrence and internal energy $U$ of the system is given by $C=1/2\max\left[0,-U/J\mathcal{N}-1\right]$ for the antiferromagnetic system (AFM) and $C=1/2\max\left[0,U/3J\mathcal{N}-1\right]$ for the ferromagnetic system (FM).\cite{Vidal}. Thus, the entanglement is uniquely determined by the partition function of the system $|\langle\sigma_z\otimes\sigma_z\rangle|=\mathcal{Z}^{-1}|\sum_ie^{-\beta E_i}\langle i|\sigma_z\otimes\sigma_z|i\rangle|\leq\mathcal{Z}^{-1}|\sum_ie^{-\beta E_i}\|\sigma_z\otimes\sigma_z|i\rangle\|\leq\mathcal{Z}^{-1}\sum_{i}e^{-\beta E_i}=1$, where we have $-1\leq U/3J\mathcal{N}\leq 1$. For the AFM case, the increase of the temperature or internal energy will generate a decreasing of the concurrence, being for a value as $-\mathcal{N}J$, the concurrence will become zero. The temperature where the concurrence vanishes is called the threshold temperature.



































Quantum entanglement is the quantum mechanical property that Schr\"{o}dinger singled out many year ago as "the characteristic trait of quantum mechanics" and that has been many analyzed in connection with Bell's inequality\cite{Ballentine,Divincenzo,Fabio,Dan,Michael}.  A pure pair of quantum systems is called entangled whether it is unfactorable, in convey a mixed state is entangled if it can not be represented as a mixture of factorable pure states. It is a well known fact that the quantum information theory can be used together condensed matter physics in characterizing of quantum phase transitions (QPT) that are characterized by the ground-state energy  of  quantum many-particle systems. Thus, the quantifying of quantum correlations in these many-body systems enhances  condensed matter physics and quantum information theory, being the measure of quantum correlation or entanglement in a system  given by the von Neumann entropy.\cite{Latorre3}.
\begin{figure}
    \centering
%\includegraphics[width=6.0cm]{von_neuman} \\
\caption{von Neumann entropy $S(\rho)$  as a function of $T$ for the model Eq.~(\ref{model}), for different strength of interaction $t=1$, $U=-0.8$, $U=-2.0$ and $U=-4.5$. The parent Hermitian system can be obtained by assuming $it\rightarrow t$. For the Hermitian system, the bound pair with lowest energy lies in the $\mathbf{k}=0$ subspace while the ground state of two-particle non-Hermitian setting locates on the subspace indexed by $\mathbf{k}=\pi$. }\label{fig_37}
\end{figure}

\textit{Von Neumann entropy}: The von Neumann entropy (VN) of a density matrix $\rho$ is defined by the formula\cite{Fabio,Dan,Michael} $\mathcal{S}\equiv-\mathrm{Tr}\left(\rho\log_2\rho\right)=-\sum_{j}r_j\log_2 r_j$, where $\rho=\sum_{j}r_j|r_j\rangle\langle r_j|$, $r_j\geq 0$, $\sum_{j}r_j=1$, $\langle r_i|r_j\rangle=\delta_{ij}$. We can define the quantum version of the entropy by the relative entropy $\mu$ to $\nu$ defined by $\mathcal{S}\left(\mu||\nu\right)\equiv\mathrm{Tr}\left(\mu\log_2\mu\right)-\mathrm{Tr}\left(\nu\log_2\nu\right)$, where the quantum relative entropy is non-negative $\mathcal{S}\left(\mu||\nu\right)\geq 0$, with equality if and only if $\mu=\nu$. The VN entropy is a quantifier of  entanglement between two different partitions of a system nominated as $A$ and $B$.  The ground state $|\Psi\rangle_{AB}$ belongs to a Hilbert space composed $\mathcal{H}=\mathcal{H}_A\otimes\mathcal{H}_B$. Thus, from the Schmidt's decomposition procedure, we can write $|\Psi\rangle_{AB}=\sum_{i}^{\mathcal{N}}r_i|\varphi_i\rangle_A\otimes|\phi_i\rangle_B$, where $\alpha_i$ are the Schmidt's coefficients and $\mathcal{N}\leq\min(\dim\mathcal{H}_A,\mathcal{H}_B)$.
In following, the whole system is considered as  a binary system,  being the block of $\mathcal{N}$ spins as sub-system $A$ and the rest of the chain as sub-system $B$ \cite{A,Latorre,Dagmar}.
Thus, the VN entropy between the two partitions is defined by $\mathcal{S}_A=\mathcal{S}_B\equiv-\sum_{i}r_i^2\log_2r_i^2$,
where $\mathcal{S}_{A}$ is the entropy of the subsystem $A$ and
$\mathcal{S}_{A}=\mathcal{S}(\rho_A)\equiv -\mathrm{Tr}_{A}(\rho_{A}\log_2\rho_{A})$.

In general, the Gibbs distribution $\rho$ is given by $\rho\propto e^{-\beta \mathcal{H}}$, where the statistical ensemble describing the system for long time is expected to be the canonical ensemble, being the density matrix of the canonical ensemble given by\cite{Calabrense,Davide,Vidal,Calabrense2}\\ $\rho=\frac{e^{-\sum_{_\mathbf{k}}\omega_\mathbf{k}\hat{n}_{\mathbf{k}}}}{\mathcal{Z}}$,
where $\mathcal{Z}$ is the partition function of the system. Thus, we have the following expression for the thermal entanglement for the non-symmetrical or non-Hermitian (non-self adjunct) model of  noninteracting fermions
%$\rho=\frac{e^{-\sum_{\mathbf{k}}\omega_\mathbf{k}\hat{n}_{\mathbf{k}}}}{\mathcal{Z}}$, being $\mathcal{Z}$ the partition function given by\\
$\mathcal{Z}=\mathrm{Tr}\left(e^{-\sum_{\mathbf{k}}\omega_\mathbf{k}\hat{n}_{\mathbf{k}}}\right)=\prod_{\mathbf{k}}(1+e^{-\beta\omega_{\mathbf{k}}})$, where we consider  $\beta=1$.
By using the form of $\mathcal{Z}$, we obtain  the thermodynamical potential as $\phi=-\frac{1}{\beta}\log_2\mathcal{Z}$ and hence
\begin{equation}{\label{ent}}
\mathcal{S}=\sum_{\mathbf{k}}\log_2(1+e^{-\beta\omega_{\mathbf{k}}})-\beta\sum_{\mathbf{k}}\frac{\omega_{\mathbf{k}}e^{-\beta\omega_{\mathbf{k}}}}{(1+e^{-\beta\omega_{\mathbf{k}}})}.
\end{equation}
Thus, the density entropy $\mathfrak{s}(\mathbf{k})$  that each mode contributes to the thermodynamic entropy in the thermodynamic limit is given by  $\mathcal{S}/\mathcal{N}^2\equiv \int_{-\pi}^{\pi}\int_{-\pi}^{\pi}\mathfrak{s}(\mathbf{k})d^2k$. Moreover, in the infinite time limit, the thermodynamic entropy and the VN entropy have the same density, representing the contribution that each mode has to the entropy of the density matrix. Consequently, in the large time limit, the results of the VN entropy  must be equal to  classical thermodynamic entropy since the results of the approach used are accurate in this limit $\mathcal{N}\rightarrow\infty$, which the identification between the two kinds of entropies can be made\cite{Calabrense2}.
%
\begin{figure}
    \centering
%\includegraphics[width=5.5cm]{von_neuman22} \\
\caption{$S(\rho)$  as a function of the complex hopping $U$ for the model Eq.~(\ref{model}).}\label{fig_38}
\end{figure}

In Figs.~\ref{fig_37}, we get the entropy of the density matrix as a function of $T$ given by Eq.~(\ref{ent}). The results are obtained for different strength of interaction $t=1$, $U=-0.8t$, $U=-2.0t$ and $U=-4.5t$. In a general way, the entanglement  entropy is a good entanglement measure if the whole system is in a pure state. However, non-hermitian Hamiltonians do describe open systems, coupled to the environment and therefore, are in  a mixed state. So the entanglement entropy is not a good measure in the investigated situation. Of course, the von Neumann entropy of the density matrix can be defined and studied in any case, it is often a useful quantity, but is not entanglement entropy in this case. Therefore, in the situation described here we have the entropy of the density matrix. The parent Hermitian system can be obtained by assuming $it\rightarrow t$. For the Hermitian system, the bound pair with lowest energy lies in the $\mathbf{k}=0$ subspace while the ground state of the two-particle non-Hermitian setting locating on the subspace indexed by $\mathbf{k}=\pi$. The presence of the imaginary hopping not only makes scattering energy bands imaginary but also reverses the whole bound band. Such non-Hermiticity alters significantly the paring mechanism and hence favors superconductivity.  We get a divergence of  $S(\rho)$ at $T\rightarrow 0$ due to large rising of the quantum fluctuations near to $T=0$ and hence, a lost of quantum information at this limit. The behavior at range of high $T$ is only qualitative due to limitations of the approach used.  The behavior of the quantum correlation given by the von Neumann entropy is determined by the behavior of the energy bands that depend on  imaginary hopping $U$ and which generates a large effect on quantum entanglement. In Fig.~\ref{fig_38}, we obtain $S(\rho)$  as a function of the strength of $U$ for the model Eq.~(\ref{model}). The behavior is monotonically decreasing  for the range of small values of complex hopping up to $U\approx3.5t$ when the behavior change suavely and $S(\rho)$ start to rise. This behavior is a consequence of the behavior of imaginary dispersion relation of quasiparticles at this range where occurs a larger lost of quantum information for large values of $U$.


\textit{Entanglement negativity}: The entanglement negativity is the linear and partial transpose whose the trace norm is convex and monotone function however, not additive. Besides, it present a large deficiency i.e. a failure in satisfying the discriminant property, either that the entanglement $E(\rho)=0$ if and only if $\rho$ is separable\cite{K}.  The entanglement negativity\cite{K,G,Laflorence} is given for a mixed state $\rho_{GE}$ by
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  {N}(\rho) = \frac{\|\rho_A^{T}\|_1-1}{2},
\end{eqnarray}
where $\rho_A^{T}$ is the partial transpose of $\rho_{GE}$ with respect to the subsystem $A$, and $\|\cdot\cdot\cdot\|_1$ is the trace norm. Furthermore, the logarithmic negativity
\begin{equation}
E_{N}(\rho)=\log_2\|\rho_A^{T}\|_1,
\end{equation}
is used much often as a measure of thermal entanglement for disjoint intervals\cite{Plenio}. Consequently, the negativity has been proven to be useful to detect topological order,\cite{C,Y} where one makes $\rho_A=\rho_{GE}$, being the result for the entanglement negativity obtained analytically as
\begin{eqnarray}{\label{neg}}
E_{{N}}(\rho)=-\log_2\bigg\|\frac{e^{-\beta\sum_{_\mathbf{k}}\omega_\mathbf{k}\hat{n}_{\mathbf{k}}}}{\mathcal{Z}}\bigg\|\nonumber\\
=\beta\sum_{\mathbf{k}}\omega_{\mathbf{k}}\hat{n}_{\mathbf{k}}+\log_2\prod_{\mathbf{k}}(1+e^{-\beta\omega_{\mathbf{k}}}).%\nonumber\\
%=\sum_{\mathbf{k}}\frac{\beta\omega_\mathbf{k}}{(1+e^{-\beta\omega_{\mathbf{k}}})}+\sum_{\mathbf{k}}\log_2(1+e^{-\beta\omega_{\mathbf{k}}}).
\end{eqnarray}
%
\begin{figure}
    \centering
%\includegraphics[width=5.5cm]{Hubbard} \\
\caption{Entanglement negativity $E_{{N}}(\rho)$,  as a function of $T$ for the model Eq.~(\ref{model}), for different strength imaginary hopping $t=it$ such as: $t=1$, $U=-0.8t$, $U=-2.0t$ and $U=-4.5t$. The negativity entanglement suffers a rising with $T$ for all values of complex hopping $U$.  }\label{fig_35}
\end{figure}

In Fig.~\ref{fig_35}, we present the behavior of the entanglement negativity given by Eq.~(\ref{neg}) for the model Eq.~(\ref{model}) as a function of $T$. The aim here is to verify the effect of splitting of the magnon modes introduced by the imaginary hopping ($it=t$) on quantum correlation. As we can see, the quantum correlation rises with the increase of the magnon splitting. The same behavior occurs with the increase of the chemical potential $\mu$ where the magnon splitting becomes higher. The cross of the curves at range $T\approx0.7J$ for low values of $U$, reflects in an equality of the imaginary dispersion relation  at this range.
%In Fig.~\ref{fig_36}, we present the behavior of the entanglement negativity for the model Eq.~(\ref{model}) as a function of interaction potential $U$. We obtain a very fast oscillation of the negativity with potential $\mathcal{U}$.





\bibitem{Brataas} Arne Brataas, Bart van Wees, Olivier Klein, Gr\'egoire de Loubens, Michel Viret, Physics Reports 885, (2020) 1.
\bibitem{JJ} J. Holanda, O. Alves Santos, J. B. S. Mendes, S. M. Rezende, J. Phys. Condens. Matter 33, (2021) 435803.
\bibitem{EE} E. Erlandsen, A. Sudbo, Phys. Rev. B 105, (2022) 184434.
\bibitem{Anderson} P. W. Anderson, Science 235, (1987) 1196.
\bibitem{Raul} R. Soni, N. Kaushal, C. Sen, F. A. Reboredo, A. Moreo, E. Dagotto, New J. Phys. 24, (2022) 073014.
\bibitem{white1} Steven R. White, Phys. Rev. Lett. 69, (1992) 2863.
\bibitem{white2} Steven R. White, Phys. Rev. B 48, (1993) 10345.


\textit{Metal-insulting antiferromagnetic bilayer model}:
Antiferromagnetic insulators have had interest as alternatives to ferromagnetic insulators as active components in spintronics\cite{Brataas,JJ,EE}. A model of interest consists in a bilayer structure consisting of an antiferromagnetic insulator on top of a normal metal. A voltage bias is applied to the normal metal in order to produce an electron current along the $x$ axis. The electrons interact with the spins leading to an induced magnon spin current\cite{EE}. %We consider the system illustrated in Fig.~\ref{fig_38} considering the layers to be two-dimensional and apply square lattice models. We start out from a tight-binding description of electrons hopping between lattice sites in a normal metal. For the antiferromagnetic insulator, we consider localized spins with easy-axis anisotropy, interacting with each other through a nearest-neighbor exchange interaction and a next-nearest-neighbor interaction. For a sufficiently small and isotropic Fermi surface in the normal metal, the Hamiltonian describing the electrons takes the form $\mathcal{H}_N=\sum_{\mathbf{k}\sigma}\varepsilon_{\mathbf{k}\sigma}c_{\mathbf{k}\sigma}^{\dag}c_{\mathbf{k}\sigma}$, where $\varepsilon_{k\sigma}=tk^2a^2-\mu-\sigma H_e$, where $c_{\mathbf{k},\sigma}^{\dag}$ is the operator for an electron with momentum $\mathbf{k}$ and spin $\sigma=\uparrow,\downarrow$. $t$ is the electron hopping amplitude, $a$ is the lattice constant, $\mu$ is the chemical potential, and $H_e$ is a spin-splitting field.

 The Hamiltonian describing the magnons is given by
%\begin{figure}
%    \centering
%\includegraphics[width=10.0cm]{fig2} \\
%\caption{Bilayer structure consisting of an antiferromagnetic insulator on top of a normal metal of the model Eq.~(\ref{elect}). A voltage is applied to the normal metal in order to produce an electron current along the $x$ %axis. }\label{fig_38}
%\end{figure}
\begin{eqnarray}{\label{elect}}
% \nonumber % Remove numbering (before each equation)
  \mathcal{H} =\sum_{\mathbf{k}}(\xi_{\mathbf{k}}+h_{\nu})\alpha^{\dag}_{\mathbf{k}}\alpha_{\mathbf{k}}+(\xi_{\mathbf{k}}-h_{\nu})\beta^{\dag}_{\mathbf{k}}\beta_{\mathbf{k}},
\end{eqnarray}
where we consider the lattice spacing $a=1$. The magnon spectrum is given by $\xi_{\mathbf{k}}=\sqrt{\Delta+\mu^2\mathbf{k}^2}$, where $\Delta$ is the gap in the spectrum. $h_{\nu}$ is a splitting of the magnon modes through an external field. $\alpha^{\dag}_{\mathbf{k}}$ are the creation operators for spin down and $\beta^ {\dag}_{\mathbf{k}}$ is the creation operator for a spin up. The model is represented in Fig.~\ref{fig_38}.

\begin{figure}
    \centering
\includegraphics[width=6.0cm]{metal-insulating} \\
\caption{$E_{{N}}(\rho)$ for the model Eq.~(\ref{elect}), as a function of $h$ (right-side), for a value of coupling $\mu$ as  $\mu=0.9$ and for a value of $h=1.0$, held fixed (left-side).}\label{fig_37}
\end{figure}




\textit{Two-dimensional Heisenberg model with bilinear-biquadratic-bicubic terms}: %{\label{771}}
%The interest in spin Heisenberg models with spin $S>1/2$ started many years ago where valence bond states serves as toy model related to high-$T_c$ superconductivity\cite{Anderson}. The AKLT model extended the notion of valence bond states to spins higher than $1/2$. The AKLT model for $J_2/J_1=1/3$, with $J_1$ and $J_2$ both positive, the ground state is exactly solvable and is defined as
%\begin{equation}
%\mathcal{H}_{AKLT}=\sum_{\langle i,j\rangle}\mathbf{S}_i\cdot \mathbf{S}_j+\frac{1}{3}\sum_{\langle i,j\rangle}\left(\mathbf{S}_i\cdot\mathbf{S}_j\right)^2.
%\end{equation}
The higher-order Heisenberg model for any spin-$S$ can be written generically as
\begin{equation}
\mathcal{H}=\sum_{\langle i,j\rangle}\sum_{\nu=1}^{2S}J_{\nu}\left(\mathbf{S}_i\cdot\mathbf{S}_j\right)^{\nu}.
\end{equation}
Thus, for the bilinear-biquadratic-bicubic model with $S=3/2$ the higher-order Heisenberg Hamiltonian reads:
\begin{equation}{\label{Bicubic}}
\mathcal{H}_{3/2}=\sum_{\langle i,j\rangle}\left[J_1\mathbf{S}_i\cdot \mathbf{S}_j+J_2\left(\mathbf{S}_i\cdot\mathbf{S}_j\right)^2+J_3\left(\mathbf{S}_i\cdot\mathbf{S}_j\right)^3\right].
\end{equation}
%For $S=3/2$, the diagonalization of this Hamiltonian for a two-site system gives four energy levels\cite{Raul}
%\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
%  \mathfrak{E}_1 &=& -\frac{15}{64}\left(16J_1-60J_2+225J_3\right),\hspace{0.25cm} \hbox{Singlet}\hspace{0.1cm} \\
%  \mathfrak{E}_2 &=& -\frac{11}{64}\left(16J_1-44J_2+121J_3\right),\hspace{0.25cm} \hbox{Triplet}\hspace{0.1cm} \\
%  \mathfrak{E}_3 &=& -\frac{3}{64}\left(16J_1-12J_2+9J_3\right),\hspace{0.25cm} \hbox{Quintuplet}\hspace{0.1cm}\\
%  \mathfrak{E}_4 &=& \frac{9}{64}\left(16J_1+36J_2+81J_3\right),\hspace{0.25cm} \hbox{Singlet}\hspace{0.1cm}.
%\end{eqnarray}
%\section{Results}{\label{77}}


\textit{Analysis by DMRG}: Density matrix renormalization group (DMRG) is a well known numerical technique suited to treat the one-dimensional spin-1/2 Heisenberg model \cite{white1,white2}. However, any finite two-dimensional lattice can be mapped in a one-dimensional lattice where the sites of the lattice are numbered and therefore long range interactions are introduced. Since in mean field theories, the dynamics of the operators $\vec{S}(\vec{r})$ is omitted, the variational principle  assumes an expectation value of the operator $\langle\vec{S}(\vec{r})\rangle$ where we neglect the fluctuations and drop higher-order terms.


\begin{figure}
    \centering
\includegraphics[width=6cm]{Bicubi_Conc2} \\
\caption{Plot of $C$ vs $J_3$ ($J_3<0$) for the model Eq.~(\ref{Bicubic}) as a function of $J_{3}$  obtained by DMRG. We obtain a small change in the behavior for different values of $J_2$. The calculations were performed for a lattice size $\mathcal{N}=192$ sites.}\label{fig_55}
\end{figure}
\noindent

In Fig.~\ref{fig_55}, we present the concurrence $C$ as a function of  $J_{3}$ ($J_3<0$), for different values of $J_2$ ($J_2<0$) coupling using DMRG. The sign of biquadratic and bicubic strength ($J_2$ and $J_3$) depends on relation between energies $\mathfrak{E}_1$, $\mathfrak{E}_2$, $\mathfrak{E}_3$ and $\mathfrak{E}_4$. We have  %$E_3-3E_2+2E_1<0$, we have  $J_2><0$ and $E_3-3E_2+2E_1>0$, $J_2>0$ of
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  \frac{J_2}{J_1} &=\frac{4}{3} \frac{29\mathcal{E}_4-85\mathcal{E}_3+81\mathcal{E}_2-25\mathcal{E}_1}{(81\mathcal{E}_4+115\mathcal{E}_3-351\mathcal{E}_2+155\mathcal{E}_1)} \\
 \frac{J_3}{J_1} &=\frac{4}{3} \frac{\mathcal{E}_4-5\mathcal{E}_3+9\mathcal{E}_2-5\mathcal{E}_1}{(81\mathcal{E}_4+115\mathcal{E}_3-351\mathcal{E}_2+155\mathcal{E}_1)},
\end{eqnarray}
where $\mathcal{E}_{\alpha}=\mathfrak{E}_{\alpha}+\mathfrak{E}_{off}$, being $\mathfrak{E}_{off}$ an offset energy\cite{Raul}. The calculations were performed for a lattice size $\mathcal{N}=192$. We obtain a very small variation of the results with different lattice sizes $\mathcal{N}=64,128,192$. Furthermore, we obtain a strong influence of the strength of the bicubic term  $J_{3}$ on concurrence with the extinction of $C$ near to $J_3\approx1.1$.  Thus, the analysis by DMRG in Fig.~\ref{fig_55} seems to confirm the decreasing of the concurrence $C$ with the coupling $J_{3}$.















\section{XXZ model on triangular lattice}{\label{modelo}}

The XXZ model on triangular lattice with magnetic field along the $z$ axis is given by the spin Hamiltonian
\begin{equation}\label{model}
  \mathcal{H}=\sum_{\langle i\alpha,j\beta\rangle}\left[J\mathbf{S}_{i\alpha}\cdot\mathbf{S}_{j\beta}+\Delta S_{i\alpha}^{z}S_{j\beta}^{z}\right]-\sum_{i\alpha}H_{\alpha}S_{i\alpha}^{z},
\end{equation}
where $\Delta>0$ indicates an anisotropic exchange coupling. The notation $\alpha,\beta=A,B,C$, denote the three sublattices in consideration. The last term is an external applied magnetic field. For $0<H<3$, the system is in $Y$ phase. For $3<H<H_{cl}$, the system is in the up-up-down phase (UUD). The phase boundary between UUD phase and V phase is given by $H_{cl}=\frac{3}{2}\left[1+2\Delta+\sqrt{1+12\Delta+4\Delta^2}\right]$. For $H_{cl}<H<9+2\Delta$, the system is in the $V$ phase and for  $H>9+2\Delta$, the system is in fully polarized phase (FP). If we make the magnetic field on $B$ sublattice as $H_B=H+\varepsilon$, $\varepsilon\ll 1$, the $Y$ and $V$ phases are replaced by distorted ones.\\
\textit{Linear spin wave theory}: In the linear spin wave approach (LSW), we perform a local rotation in the coordinate system at each sublattice at each lattice point, so that the mean-field directions of the spins point along the local $z$ axis
\begin{eqnarray}
\mathbf{S}_{j\alpha}=\left(
  \begin{array}{ccc}
    \cos\theta_{\alpha} & 0 & \sin\theta_{\alpha} \\
    0 & 1 & 0 \\
    -\sin\theta_{\alpha} & 0 & \cos\theta_{\alpha} \\
  \end{array}
\right)\tilde{\mathbf{S}}_{j\alpha}.
\end{eqnarray}
In the  harmonic approximation, we perform the Holstein-Primakoff transformation %expanded up to first order in powers of $1/S$
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
\hspace{-0.3cm}  \mathbf{\tilde{S}}_{j\alpha}^{+} \approx\sqrt{2S}a_{j\alpha}^{\dag}\left(1-\frac{a_{j\alpha}^{\dag}a_{j\alpha}}{2S}\right),\hspace{0.25cm}\mathbf{\tilde{S}}_{j\alpha}^{-} \approx\sqrt{2S}a_{j\alpha},\hspace{0.25cm}\tilde{S}^z_{j\alpha}=S-a_{j\alpha}^{\dag}a_{j\alpha},\nonumber\\
\end{eqnarray}
where $a^{\dag}_{j\alpha}(a_{j\alpha})$ are creation (annihilation) boson operators. The Hamiltonian is written in the momentum space as
\begin{eqnarray}{\label{boson}}
\mathcal{H}=\mathcal{E}_0+\frac{S}{2}\sum_{\mathbf{k}}\left(
                                              \begin{array}{cc}
                                               \phi^{\dag}(\mathbf{k}) & \bar{\phi}(-\mathbf{k}) \\
                                              \end{array}
                                            \right)
\mathcal{H}(\mathbf{k})\left(
                         \begin{array}{c}
                           \phi(\mathbf{k}) \\
                           \bar{\phi}^{\dag}(-\mathbf{k}) \\
                         \end{array}
                       \right),\nonumber\\
                       \end{eqnarray}
where $\mathcal{E}_0$ is the energy of the ground state in the mean-field approach and $\phi^{\dag}(\mathbf{k})=\left(a^{\dag}_{\mathbf{k},A},a^{\dag}_{\mathbf{k},B},a^{\dag}_{\mathbf{k},C}\right)$, where the tilde over $\phi$ means the transpose. $\mathcal{H}(\mathbf{k})$ is given by
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
\mathcal{H}(\mathbf{k})= \left(
    \begin{array}{cc}
      \mathfrak{A}(\mathbf{k}) & \mathfrak{B}(\mathbf{k}) \\
      \mathfrak{B}^{*}(-\mathbf{k}) & \mathfrak{A}^{*}(-\mathbf{k}) \\
    \end{array}
  \right),
  \end{eqnarray}
  where
  \begin{eqnarray}
  % \nonumber % Remove numbering (before each equation)
\hspace{-1.0cm} [\mathfrak{A}(\mathbf{k})]_{\alpha\alpha} = H_{\alpha}\cos\theta_{\alpha}
    -\sum_{\alpha\ne\beta} \left[3\cos(\theta_{\alpha}-\theta_{\beta})+3\Delta\cos\theta_{\alpha}\cos\theta_{\beta}\right],\nonumber \\
\hspace{-1.0cm}\left[\mathfrak{A}(\mathbf{k})\right]_{\alpha\delta}  = \left[\cos(\theta_{\alpha}-\theta_{\beta})+\Delta\sin\theta_{\alpha}\sin\theta_{\beta}+1\right]\Gamma_{\alpha\beta}(\mathbf{k}),\nonumber \\
\hspace{-1.0cm}\left[\mathfrak{B}(\mathbf{k})\right]_{\alpha\beta} = \left[\cos(\theta_{\alpha}-\theta_{\beta})+\Delta\sin\theta_{\alpha}\sin\theta_{\beta}-1\right]\Gamma_{\alpha\beta}(\mathbf{k})\nonumber\\
\hspace{-1.0cm}\left[\mathfrak{B}(\mathbf{k})\right]_{\alpha\alpha} =0,\nonumber\\
  \end{eqnarray}
where
\begin{equation}
\Gamma_{\alpha\beta}(\mathbf{k})=2\cos\left(k_x/2\right)\cos\left(\sqrt{3}k_y/2\right).
\end{equation}
 Due to time reversal symmetry, we have $\mathfrak{A}(\mathbf{k})=\mathfrak{A}^{*}(-\mathbf{k})$ and $\mathfrak{B}(\mathbf{k})=\mathfrak{B}^{*}(-\mathbf{k})$. The Hamiltonian above can be diagonalized by a paraunitary Bogoliubov transformation, finding a matrix $\mathcal{T}_{\mathbf{k}}$ such that $\omega_{\mathbf{k}}=\mathcal{T}^{\dag}_{\mathbf{k}}\mathcal{H}(\mathbf{k})\mathcal{T}_{\mathbf{k}}$, where
%\begin{equation}
$\mathcal{T}^{\dag}_{\mathbf{k}}\mathcal{H}(\mathbf{k})\mathcal{T}_{\mathbf{k}}= \hbox{diag}\left(\omega_{1,\mathbf{k}},\omega_{2,\mathbf{k}},\omega_{3,\mathbf{k}},\omega_{1,-\mathbf{k}},\omega_{2,-\mathbf{k}},\omega_{3,-\mathbf{k}}\right),$
%\end{equation}
% gives three magnon energy bands
The magnon energy bands are given by the diagonalization problem of the bosonic Hamiltonian Eq.~(\ref{boson}) \cite{colpa}
\begin{eqnarray}{\label{equa}}
\det\left(\mathcal{H}(\mathbf{k})-\omega_{\mathbf{k}}\mathbf{I}\right)=0
\end{eqnarray}
either
\begin{eqnarray}{\label{equa}}
\left|
  \begin{array}{cc}
     \mathfrak{A}(\mathbf{k})-\omega_{\mathbf{k}}) \mathbf{I} &  \mathfrak{B}(\mathbf{k}) \\
    \mathfrak{B}^{*}(-\mathbf{k}) & \mathfrak{A}^{*}(-\mathbf{k})-\omega_{\mathbf{k}} \mathbf{I} \\
  \end{array}
\right|=0,
\end{eqnarray}
in which $\mathbf{I}$ is a $m$-square unit matrix. Hence, we obtain
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  (\mathfrak{A}(\mathbf{k})-\omega_{\mathbf{k}})(\mathfrak{A}^{*}(-\mathbf{k})-\omega_{\mathbf{k}})-\mathfrak{B}(\mathbf{k})\mathfrak{B}^{*}(-\mathbf{k})=0\nonumber\\
  \omega_{\mathbf{k}}=\frac{1}{2}[(\mathfrak{A}(\mathbf{k})+\mathfrak{A}^{*}(-\mathbf{k})]\nonumber\\
  \pm\sqrt{(\mathfrak{A}(\mathbf{k})+\mathfrak{A}^{*}(-\mathbf{k}))^2-4(\mathfrak{A}(\mathbf{k})\mathfrak{A}^{*}(-\mathbf{k})-\mathfrak{B}(\mathbf{k})\mathfrak{B}^{*}(-\mathbf{k}))}].\nonumber\\
\end{eqnarray}
In a uniform magnetic field, the system is in the $Y$ phase, and the lowest magnon bands touch each other at $K$ and $K'$ points, generating massless Dirac-cone-like dispersions. The mean-field ground state is at $\theta_{0,\alpha=A}=\pi$, $\theta_{0,\alpha=B}=\theta_{0,\alpha=C}=\theta$, where
\begin{equation}
\theta=\cos^{-1}\left(\frac{3\Delta+h+3}{3\Delta+6}\right).
\end{equation}
are determined by minimizing the mean-field energy given by
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  \mathfrak{E}_{MF} = 3S^2\sum_{\langle\alpha,\beta\rangle}[(1+\Delta)\cos\theta_{\alpha}\cos\theta_{\beta}+\sin\theta_{\alpha}\sin\theta_{\beta}]\nonumber\\
  -\sum_{\alpha}H_{\alpha}\cos\theta_{\alpha},\nonumber\\
\end{eqnarray}
where $\alpha,\beta$ means summation over the pairs $(A,B)$, $(B,C)$, and $(C,A)$. $\theta_{\alpha}$ are determined by minimizing the mean-field energy $E_{MF}$, which leads to different solutions for $\theta_{\alpha}$, identifying the different phases nominated as $Y$, $V$, UUD and FP\cite{Qi}. For $0<H<3$, the system is in the $Y$ phase which we are mainly concerned. Moreover, we focused in the antiferromagnetic region $\Delta>0$, where there are four different phases in a uniform magnetic field $H_{\alpha}=H$.

\begin{figure}
    \centering
\includegraphics[width=12cm]{F1} \\
\caption{Representation of the triangular lattice, for the model Eq.~(\ref{model}), with the three sublattices denoted as $A$, $B$ and $C$. Below, we have the four different spin orientations on  $A$, $B$ and $C$ sublattices.}\label{fig_12} %$J_1=−1$, $J_2=1$, $J_{bq_1}=−3$, $J_{bq_2}=0.5$, (black),
\end{figure}

%\section{Results}{\label{5}}


\begin{figure}
    \centering
\includegraphics[width=5.0cm]{magnon_valey} \\
\caption{Entanglement negativity $E_{{N}}(\rho)$  vs. $T$ for different values of $\Delta$ at $H=1.1$, for the model Eq.~(\ref{model}). We obtain a small difference of the entanglement negatively for the bands $\omega^{+}_{\mathbf{k}}$ and $\omega^{-}_{\mathbf{k}}$ (solution of Eq.~(\ref{equa})). Furthermore, the entanglement negatively diverges at limit $T\rightarrow 0$ due to large rising of the quantum fluctuations near to $T=0$ where a quantum phase transition takes place.}\label{fig_12} %$J_1=−1$, $J_2=1$, $J_{bq_1}=−3$, $J_{bq_2}=0.5$, (black),
\end{figure}


\begin{figure}
    \centering
\includegraphics[width=6.0cm]{magnon_valey_delta} \\
\caption{$E_{{N}}(\rho)$  as a function of $\Delta$ for $H=1.1$, for the model Eq.~(\ref{model}). We obtain that the negativity entanglement is finite at XY limit ($\Delta=0$), rising with $\Delta$ up to isotropic limit $\Delta=1$.  }\label{fig_35}
\end{figure}



triangular-lattice XXZ model with three sublattices denoted by $A$, $B$ and $C$, where we focused in the antiferromagnetic region $\Delta>0$, four different phases present in an uniform magnetic field $H_{\gamma}=H$. We focused in the $Y$ phase, for $0<H<3$. The ground-state phase diagram of this spin model is obtained in the classical limit $S\rightarrow\infty$ with four different regions whose spin configurations are represented in Fig.~\ref{fig_12}.
Our results display  a divergence of quantum correlation $T=0$ further strong effect of the  anisotropy $\Delta$ on entanglement, generating an increasing of $E_{{N}}(\rho)$ from XY limit up to isotropic limit.





\section{non-Hermitian models}{\label{model}}


%\subsection{Two-dimensional non-interacting fermions models}%{\label{10}}
\subsection{Noninteracting fermions systems}

A linear Hamiltonian operator is symmetric if the domain of $\mathcal{H}$, $D(\mathcal{H})$ is dense in the Hilbert space and $\langle\mathcal{H}x,y\rangle=\langle x,\mathcal{H}y\rangle\forall x,y\in D(\mathcal{H})$. A linear operator $\mathcal{H}$ in an Hilbert space is self-adjunct is and only if $\mathcal{H}=\mathcal{H}^{\dag}$, that is, if $\mathcal{H}$ is symmetric (Hermitian) and $D(\mathcal{H})=D(\mathcal{H^{\dag}})$. The condition $D(\mathcal{H})=D(\mathcal{H}^{\dag})$ is frequently neglected where one treats the Hermitian (symmety) as sufficient condition to get the results of physics that are guaranteed to self-adjunct operators. Therefore, we have that the 2D non-Hermitian Lieb lattice (non-self adjunct) model of non-interacting fermions does not satisfy this condition. The Hamiltonian of this system is written as
\begin{eqnarray}{\label{nonher}}
% \nonumber % Remove numbering (before each equation)
%\hspace{-0.0cm} \mathcal{H} = \sum_{m,n}(\mu \alpha^{\dag}_{mn}\beta_{mn}+\mu \beta^{\dag}_{mn}\gamma_{mn}+i\nu \gamma_{mn}^{\dag}\alpha_{mn}\nonumber \\
%\hspace{-0.0cm}  +\mu \alpha^{\dag}_{mn}\beta_{m+1, n}+\mu \beta^{\dag}_{mn}\gamma_{m,n+1}+\hbox{H.c.})%\nonumber \\
%  +i\varepsilon\left(\alpha^{\dag}_{mn}\alpha_{mn}-\gamma^{\dag}_{mn}\gamma_{mn}\right),\nonumber\\
\hspace{-0.0cm} \mathcal{H} = \sum_{m,n}\{[\mu( \alpha^{\dag}_{mn}\beta_{mn}+ \beta^{\dag}_{mn}\gamma_{mn}
\hspace{-0.0cm}  + \alpha^{\dag}_{mn}\beta_{m+1, n}+ \beta^{\dag}_{mn}\gamma_{m,n+1})\nonumber \\+i\nu \gamma_{mn}^{\dag}\alpha_{mn}+\hbox{H.c.}]%\nonumber \\
  +i\varepsilon\left(\alpha^{\dag}_{mn}\alpha_{mn}-\gamma^{\dag}_{mn}\gamma_{mn}\right)\},\nonumber\\
\end{eqnarray}
where $\mu,\nu,\varepsilon \in\mathbb{R}$. In this case the Hamiltonian is rewritten as
\begin{equation}
\mathcal{H}=\sum_{\mathbf{k}}\Psi^{\dag}_{\mathbf{k}}\eta({\mathbf{k}})\Psi_{\mathbf{k}},
\end{equation}
with $\Psi_{\mathbf{k}}=\left(\alpha_{\mathbf{k}},\beta_{\mathbf{k}},\gamma_{\mathbf{k}}\right)^{T}$ and
\begin{eqnarray}
\eta({\mathbf{k}})=\left(
         \begin{array}{ccc}
           i\varepsilon & \mu(e^{ik_x}+1) & -i\nu \\
           \mu(e^{-ik_x}+1) & 0 & \mu(e^{ik_y}+1) \\
           i\nu & \mu(e^{-ik_y}+1) & -i\varepsilon \\
         \end{array}
       \right).\nonumber\\
\end{eqnarray}
The energy band of $\eta_{\mathbf{k}}$ are obtained by solving of the characteristic polynomial $\det\left(\eta_{\mathbf{k}}-\mathbf{I}\Lambda\right)=0$, ($\mathbf{I}$ is the identity matrix), which gives a cubic equation\cite{Xie}
\begin{eqnarray}
\Xi(\mathbf{k})^{3}+\left[\varepsilon^2-|\mu(e^{ik_x}+1)|^2-|\mu(e^{ik_y}+1)|^2-\nu^2\right]\Xi(\mathbf{k})\nonumber\\
+i\varepsilon\left[|\mu(e^{ik_y}+1)|^2-|\mu(e^{ik_x}+1)|^2\right]\nonumber\\+2|\mu(e^{ik_x}+1)||\mu(e^{ik_y}+1)|\nu\sin\left(\frac{k_x+k_y}{2}\right)=0.\nonumber\\
\end{eqnarray}
Thus, the energy bands are given by
\begin{eqnarray}{\label{nonherm}}
 \Xi(\mathbf{k})=\sqrt[3]{-\frac{\mathfrak{q}(\mathbf{k})}{2}+\sqrt{\Delta(\mathbf{k})}}
%  +\sqrt[3]{-\frac{q(\mathbf{k})}{2}-\sqrt{\left(\frac{q(\mathbf{k})}{2}\right)^2+\left(\frac{p(\mathbf{k})}{3}\right)^3}},
+\sqrt[3]{-\frac{\mathfrak{q}(\mathbf{k})}{2}-\sqrt{\Delta(\mathbf{k})}},
\end{eqnarray}
where
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  \mathfrak{q}(\mathbf{k})=i\varepsilon\left[|\mu(e^{ik_y}+1)|^2-|\nu(e^{ik_x}+1)|^2\right]\nonumber\\+2|\mu(e^{ik_x}+1)||\mu(e^{ik_y}+1)|\nu\sin\left(\frac{k_x+k_y}{2}\right), \nonumber\\
  \mathfrak{p}(\mathbf{k})=\varepsilon^2-|\mu(e^{ik_x}+1)|^2-|\mu(e^{ik_y}+1)|^2-\nu^2,\nonumber\\
\end{eqnarray}
being the discriminant of the equation above being given by
\begin{equation}{\label{discriminant}}
\Delta(\mathbf{k})=\left(\frac{\mathfrak{q}(\mathbf{k})}{2}\right)^2+\left(\frac{\mathfrak{p}(\mathbf{k})}{3}\right)^3.
\end{equation}
Whether  $\Delta<0$ we have $\Xi(\mathbf{k})\in\mathbb{R}$ and
we get three different energy bands for the coupling parameters $\varepsilon$, $\delta$ and $\nu$, that must satisfy the inequality
\begin{eqnarray}
\{i\varepsilon\left[|\mu(e^{ik_y}+1)|^2-|\mu(e^{ik_x}+1)|^2\right]\nonumber\\+2|\mu(e^{ik_x}+1)||\mu(e^{ik_y}+1)| \nu\sin\left(\frac{k_x+k_y}{2}\right)\}^2 /4\nonumber\\
+\left[\varepsilon^2-|\mu(e^{ik_x}+1)|^2-|\mu(e^{ik_y}+1)|^2-\nu^2\right]^3/27<0.\nonumber\\
\end{eqnarray}
For $\Delta\leq 0$, we can have imaginary roots of cubic equation above and therefore, this possibility are discarded. Consequently, the behavior of the discriminant induced by the coupling parameters of the non-Hermitian model generates a large effect on the Von Neumann entropy and entanglement negativity.
We obtain  the roots of $\Delta(\mathbf{k})$ are given by
\begin{eqnarray}
\hspace{-1.0cm}\Delta(\mathbf{k})=\sqrt{E^2+F^2}\left[\cos\left(\frac{\theta_2}{2}+n\pi\right)+i\sin\left(\frac{\theta_2}{2}+n\pi\right)\right],\nonumber\\
\Delta_{1,2}(\mathbf{k})=\pm\sqrt{E^2+F^2}\left[\cos\left(\frac{\theta_2}{2}\right)+i\sin\left(\frac{\theta_2}{2}\right)\right],\nonumber\\
%\Xi_2(\mathbf{k})=\sqrt{E^2+F^2}\left[\cos\left(\frac{\theta_2}{2}\right)+i\sin\left(\frac{\theta_2}{2}\right)\right],\nonumber\\
\end{eqnarray}
where $n\in\mathbb{Z}$. Hence, the energy bands are given by
\begin{eqnarray}{\label{energybands}}
\Xi(\mathbf{k})=2\sqrt[3]{c(\mathbf{k})^2+d(\mathbf{k})^2}\cos\left(\frac{\theta_1+2n\pi}{3}\right),%-i\sin\left(\frac{\theta_1+2n\pi}{3}\right)\nonumber\\
%-\cos\left(\frac{\theta_1+2n\pi}{3}\right)+i\sin\left(\frac{\theta_1+2n\pi}{3}\right)\bigg],\nonumber\\
\end{eqnarray}
 where $-\pi<\frac{\theta_1+2n\pi}{3}<\pi$. Moreover,
\begin{eqnarray}
\hspace{-0.25cm}c(\mathbf{k})=\sqrt{E^2+F^2}\cos\left(\frac{\theta_2}{2}\right)-\frac{A(\mathbf{k})}{2}\nonumber\\
\hspace{-0.25cm}d(\mathbf{k})=\sqrt{E^2+F^2}\sin\left(\frac{\theta_2}{2}\right)-\frac{B(\mathbf{k})}{2},\nonumber\\
\hspace{-0.25cm}\theta_1=\cos^{-1}\left(\frac{c(\mathbf{k})}{\sqrt{c(\mathbf{k})^2+d(\mathbf{k})^2}}\right),\hspace{0.15cm}
\tan\theta_1=d(\mathbf{k})/c(\mathbf{k}),\nonumber\\%\hspace{0.15cm}-\pi<\frac{\theta_1}{2}<\pi,\nonumber\\
\end{eqnarray}
 where
\begin{eqnarray}
\hspace{-0.0cm}\theta_2=\cos^{-1}\left(\frac{E(\mathbf{k})}{\sqrt{E(\mathbf{k})^2+F(\mathbf{k})^2}}\right),\hspace{0.5cm}\tan\theta_2=F(\mathbf{k})/E(\mathbf{k}),\nonumber\\
\hspace{-0.0cm}E(\mathbf{k})=\frac{1}{4}\left[A(\mathbf{k})^2+B(\mathbf{k})^2\right]+\frac{1}{27}\left[C(\mathbf{k})^3-C(\mathbf{k})D(\mathbf{k})^2\right],\nonumber\\
\hspace{-0.0cm}F(\mathbf{k})=-\frac{1}{27}\left[3C(\mathbf{k})^2D(\mathbf{k})+D(\mathbf{k})^3\right],\nonumber\\
\hspace{-0.0cm}\mathfrak{p}(\mathbf{k})=A(\mathbf{k})+iB(\mathbf{k}),\nonumber\\
\hspace{-0.0cm}\mathfrak{q}(\mathbf{k})=C(\mathbf{k})-iD(\mathbf{k}),\nonumber\\
\hspace{-0.0cm}A(\mathbf{k})=\varepsilon\nu^2\left(\sin 2k_x-\sin 2k_y+2\sin k_x-2\sin k_y\right)\nonumber\\+2\mu\nu^2\sin\left(\frac{k_x+k_y}{2}\right)\left[\cos(k_x+k_y)+\cos k_x+\cos k_y+1\right],\nonumber\\
\hspace{-0.0cm}B(\mathbf{k})=\varepsilon\nu^2\left(\cos 2k_y-\cos 2k_x+2\cos k_y-2\cos k_x\right)\nonumber\\+2\mu\nu^2\sin\left(\frac{k_x+k_y}{2}\right)\left[\sin(k_x+k_y)+\sin k_x+\sin k_y+1\right],\nonumber\\
\hspace{-0.0cm}C(\mathbf{k})=\varepsilon^2-\nu^2\left(\cos 2k_x+\cos 2k_y+2\cos k_x+2\cos k_y+2\right)-\mu^2\nonumber\\
\hspace{-0.0cm}D(\mathbf{k})=\nu^2\left(\sin 2k_x+\sin 2k_y+2\sin k_x+2\sin k_y\right).\nonumber
\end{eqnarray}
%Hence, we have the modes energy are
The behavior of $\Xi(\mathbf{k})$ induced by the coupling parameters of the non-Hermitian model generates a large effect on quantifiers quantum entanglement as Von Neumann entropy and entanglement negativity.

\subsection{Tight-binding model}% of fermions noninteracting}%{\label{7}}

The two-dimensional LL lattice of non-interacting fermions in a square lattice on the Lieb lattice is characterized by three atoms per unit cell. Introducing the fermion operators $\left(\alpha^{\dag}_{lm},\beta^{\dag}_{lm},\gamma^{\dag}_{lm}\right)$, the Hamiltonian in a perpendicular magnetic field in the momentum space is given by\cite{Nija}
\begin{eqnarray}
\hspace{-0.0cm}\mathcal{H}=\sum_{\mathbf{k}}\left(
                                                                                                             \begin{array}{ccc}
                                                                                                               \alpha^{\dag}_{\mathbf{k}} & \beta^{\dag}_{\mathbf{k}} & \gamma^{\dag}_{\mathbf{k}} \\
                                                                                                             \end{array}
                                                                                                           \right)
\left(
                                                                                                             \begin{array}{ccc}
                                                                                                               0 &\mathfrak{f}^{*}_{x}(\mathbf{k})  &\mathfrak{f}^{*}_{y}(\mathbf{k})  \\
                                                                                                              \mathfrak{f}^{}_{x}(\mathbf{k})  & 0 & 0 \\
                                                                                                              \mathfrak{f}^{}_{y}(\mathbf{k})  & 0 & 0 \\
                                                                                                             \end{array}
                                                                                                           \right)\left(
                                                                                                                    \begin{array}{c}
                                                                                                                      \alpha_{\mathbf{k}} \\
                                                                                                                      \beta_{\mathbf{k}} \\
                                                                                                                      \gamma_{\mathbf{k}} \\
                                                                                                                    \end{array}
                                                                                                                  \right),\nonumber\\
\end{eqnarray}
where $\mathfrak{f}_x(\mathbf{k})=\mu_x(1+e^{ik_x})$ and $\mathfrak{f}_y(\mathbf{k})=\mu_y(1+e^{ik_y})$.  We perform the Fourier transform given by\\ $\gamma_{\mathbf{k}}=\frac{1}{\sqrt{NM}}\sum_{l,m=1}^{\mathcal{N}}\gamma_{lm}e^{i\mathbf{k}\cdot \mathbf{r}_l}$, assuming that the lattice is composed of $\mathcal{N}=NM$ cells. We have the following eigenvalues of the model
\begin{equation}
\Omega_{\pm}({\mathbf{k}})=\pm2\sqrt{\mu_x^2\cos^2\left(\frac{k_{x}}{2}\right)+\mu^2_y\cos^2\left(\frac{k_{y}}{2}\right)}, \hspace{0.5cm}\Omega_{0}({\mathbf{k}})=0,
\end{equation}
where $\Omega_{\pm}$ are the energies of the upper and lower bands, respectively and $\Omega_0$ is the flat (non-dispersive) band of the LL lattice. The point $\Gamma=\left(\pi,\pi\right)$ is the  most interesting point into the Brillouin zone, where in the continuum limit or infinite lattice the three branches touch each other. The energies are given by to a massless spectrum $\Omega_{\pm}({\mathbf{k}})=\pm\sqrt{\mu_x^2k_{x}^2+\mu_y^2k_y^2}$, where $k_{x}$ and $k_y$ are measured from the point $\Gamma$. The expansion of the same functions at $R=(0,0)$ give to a parabolic dependence $\Omega_{\pm}({\mathbf{k}})=\pm\left(\frac{k_{x}^2}{2\mu_x}+\frac{k_y^2}{2\mu_y}\right)$, where $\mu_x$ and $\mu_y$ are the effective masses along the two directions. For $\mu_x=\mu_y=\mu$, we have effective mass exhibiting opposite signs generating a change of sign of the Hall effect.

\subsection{Non-Hermitian XXZ model}
The corresponding Hamiltonian is given by
\begin{equation}\label{h2}
  \mathcal{H}=-\frac{1}{2}\sum_{j}\left(\sigma_j^{+}\sigma_{j+1}^{-}+\sigma_j^{-}\sigma_{j+1}^{+}\right)+\Delta\sum_{j}\sigma_j^{z}\sigma_{j+1}^{z},
\end{equation}
where $\sigma_j$ are the Pauli's matrices. The system is in the ferromagnetic Ising phase when $\Delta > 1$. The system presents two degenerate ground states $|\uparrow\rangle$, $|\downarrow\rangle$, where the transformation $\sigma_j^z=-\sigma_j^z$ of the system breaks. There is a superposition of the two ground states whit the external magnetic field as $|\Psi\rangle=\pm\sqrt{1-\delta}|\uparrow\rangle+\sqrt{1+\delta}|\downarrow\rangle$  and when $|\delta|=1$, the two states coalesce to either.\cite{X,Shiyao}
\subsubsection{Non-Hermitian Ising model}

The model is given by
\begin{equation}
\mathcal{H}=\sum_{j=1}^{N}\sigma_{j}^z\sigma_{j+1}^z+\lambda(\sigma_{j}^x+i\delta \sigma_j^y),
\end{equation}
for periodic boundary condition $\sigma_j^{x}=\sigma_{j+L}^{x}$, $\sigma_j^{y}=\sigma_{j+L}^{y}$ and $\sigma_j^{z}=\sigma_{j+L}^{z}$. The Hamiltonian can be transformed in the form
\begin{equation}
\tilde{\mathcal{H}}=\sum_{j=1}^{N}T_{j}^zT_{j+1}^z+\lambda\sqrt{1-\delta^2}(T_{j}^x+i\delta T_j^y),
\end{equation}
which is the standard Ising model with field $\lambda\sqrt{1-\delta^2}$, with $|\delta|\ne 1$. Performing the Jordan-Wigner transformation
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  T_j^x =\frac{1}{2}-\bar{f}_jf_j, \\
  T_j^y =\frac{i}{2}\sum_{j<l}(1-2\bar{f}_jf_j)(\bar{f}_j-f_j), \\
  T_j^{z}=-\frac{1}{2}\sum_{j<l}(1-2\bar{f}_jf_j)(\bar{f}_j+f_j)
\end{eqnarray}
where $f_j$ and $\bar{f}_j$ are non-Hermitian operators: $\bar{f}_j=\mathfrak{S}_jc_j^{\dag}\mathfrak{S}_j^{-1}$, $f_j=\mathfrak{S}_jc_j\mathfrak{S}_j^{-1}$, with $[\bar{f}_j,f_j']=\delta_{j,j'}$ and $c_j^{\dag}$, $c_j$ are the creation and annihilation operators of spinless fermion. The parity of the number of such fermions is a conservative quantity such that the Hamiltonian can be expressed as
$\tilde{\mathcal{H}}=\tilde{\mathcal{H}}_{+}\mathbf{I}=\tilde{\mathcal{H}}_{-}\mathbf{I}$, where $\tilde{\mathcal{H}}_{+}=\tilde{\mathcal{H}}_{-}=-2(\bar{f}_{N}\bar{f}_1+\bar{f}_{N}f_1+\bar{f}_1f_{N}+f_1f_{N})$ and the Hamiltonian can be expressed as
\begin{eqnarray}
\bar{\mathcal{H}}=\frac{1}{4}\sum_{j=1}^{N}2[\lambda\sqrt{1-\delta^2}(1-2\bar{f}_jf_j)\nonumber\\
+(\bar{f}_j\bar{f}_{j+1}+\bar{f}_jf_{j+1}+\bar{f}_{j+1}f_{j}+f_{j+1}f_{j})].
\end{eqnarray}
Taking the Fourier transformation
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  f_j=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}f_{\mathbf{k}}e^{i{\mathbf{k}\cdot\mathbf{r}_j}}, \hspace{0.5cm}\bar{f}_j=\frac{1}{\sqrt{N}}\sum_{\mathbf{k}}f_{\mathbf{k}}e^{-i{\mathbf{k}\cdot\mathbf{r}_j}},
\end{eqnarray}
where $|\mathbf{k}|=k=2\pi(n+1/2)/N$, $n=0,1,2,...,N-1$. The Hamiltonian can be written as
\begin{equation}\label{ham}
  \bar{\mathcal{H}}_{+}=\bar{\mathcal{H}}_{-}=\sum_{0<k<\pi}\bar{\psi}_{\mathbf{k}}\bar{\mathcal{H}}_{+}^{\mathbf{k}}\psi_{\mathbf{k}},
\end{equation}
with $\bar{\psi}_{\mathbf{k}}=(\bar{f}_{\mathbf{k}},f_{-\mathbf{k}})$, $\psi_{\mathbf{k}}=(f_{\mathbf{k}},\bar{f}_{-\mathbf{k}})^{T}$. Making the non-Hermitian Bogoliubov transformation
\begin{eqnarray}
% \nonumber % Remove numbering (before each equation)
  \bar{\xi}_{\mathbf{k}}=\cos\frac{\chi_{\mathbf{k}}}{2}\bar{f}_{\mathbf{k}}+i\sin\frac{\chi_{\mathbf{k}}}{2}f_{-\mathbf{k}}, \\
  \xi_{\mathbf{k}}=\cos\frac{\chi_{\mathbf{k}}}{2}f_{\mathbf{k}}-i\sin\frac{\chi_{\mathbf{k}}}{2}\bar{f}_{-\mathbf{k}},
\end{eqnarray}
where $[\bar{\xi}_{\mathbf{k}},\xi_{\mathbf{k'}}]=\delta_{\mathbf{k},\mathbf{k}'}$ and $\chi_{\mathbf{k}}=\arctan[\sin(k)/(2\lambda\sqrt{1-\delta^2}-\cos(k))]$, where the Hamiltonian is recast in diagonal form: $\bar{\mathcal{H}}_{+}=\sum_{\mathbf{k}}\omega_{\mathbf{k}}\left(\bar{\chi}_{\mathbf{k}}\chi_{\mathbf{k}}-1/2\right)$,  with the dispersion relation of quasi-particles given by
\begin{equation}
\omega_{\mathbf{k}}=\sqrt{4\lambda^2(1-\delta^2)-4\lambda\cos(\mathbf{k})\sqrt{1-\delta^2}+1},
\end{equation}
which is a noninteracting Hamiltonian. If $|\delta|< 1$, the single-particle energy is real and $|\delta|> 1$, the system respects a complex single-particle spectrum regardless of $\mathbf{k}$.
\begin{figure}
    \centering
%\includegraphics[width=6.0cm]{nonhermitian.eps} \\
\caption{Entanglement negativity $E_{{N}}(\rho(T))$ vs. $T$ for different values of non-Hermitian parameter $\delta$. The inset display the different curves for different values of $\delta$, where the curves tends to touch each other at limit of Hermitian model where $\delta\rightarrow 0$.}\label{fig_11} %$J_1=−1$, $J_2=1$, $J_{bq_1}=−3$, $J_{bq_2}=0.5$, (black),
\end{figure}




%\subsection{Numerical results}
%\subsection{VN entropy for the non-Hermitian models}
 We must have a behavior for the VN entropy $E$ vs. $T$ only when the discriminant is $\Delta<0$, or when the three eigenvalues given by roots of the cubic equation,  Eq.~(\ref{nonherm}) are real. We obtain that all curves touch on at $T\rightarrow 0$ limit. At Dirac cone, were the spectrum is massless and the dispersion relation tends to $\Omega_{\pm}({\mathbf{k}})=\pm\sqrt{\mu_x^2k_x^2+\mu_y^2k_y^2}$. $k_x$ and $k_y$ are measured at point $\Gamma$. Expanding the same function at point $R=(0,0)$, we obtain a parabolic dependence: $\Omega_{\pm}=\pm\left(\frac{k_x^2}{2\mu_x}+\frac{k_y^2}{2\mu_y}\right)$, where $\mu_x$ and $\mu_y$ are the effective masses along the two directions,\cite{Nija} where the curves becomes more damped.




























, given by $\mathcal{Z}=\mathrm{Tr}e^{-\sum_{\mathbf{k}}\Lambda_\mathbf{k}\hat{n}_{\mathbf{k}}}=\prod_{\mathbf{k}}(1-e^{-\Lambda_{\mathbf{k}}})^{-1}$.
By using $\mathcal{Z}$ given above, one obtains
\begin{eqnarray}{\label{entangl}}
E=S_{GE}=-\mathrm{Tr}(\rho_{GE}\log_2\rho_{GE}).%=-\mathrm{Tr}\frac{e^{-\sum_{\mathbf{k}}\xi_{\mathbf{k}} \hat{n}_{\mathbf{k}}}}{Z}\ln\frac{e^{-\sum_{\mathbf{k}}\xi_{\mathbf{k}} \hat{n}_{\mathbf{k}}}}{Z}.%\nonumber\\
%=-\sum_{\mathbf{k}}\xi_{\mathbf{k}}\frac{\partial\ln Z}{\partial\xi_{\mathbf{k}}}+\ln Z.\nonumber\\
\end{eqnarray}





\subsection{Two-dimensional antiferromagnet on the Lieb lattice}

 The nearest-neighbor Heisenberg antiferromagnet on the Lieb lattice is described by the Hamiltonian
\begin{eqnarray}\label{model}
\hspace{-0.0cm}\mathcal{H}=-J_1\sum_{\langle i,j\rangle} \mathbf{S}_i\cdot \mathbf{S}_j\nonumber\\-J_2\sum_{\langle\langle i,j\rangle\rangle} \mathbf{S}_i\cdot \mathbf{S}_j
-\sum_{\langle\langle i,j\rangle\rangle}\mathbf{D}_{ij}\cdot\left(\mathbf{S}_i\wedge\mathbf{S}_j\right)
-h\sum_{i}\mathbf{S}^{z}_i,\nonumber\\
\hspace{0.0cm}\mathcal{H}=-J_1\sum_{\langle i,j\rangle} \mathbf{S}_i\cdot \mathbf{S}_j\nonumber\\
\hspace{0.0cm}-\frac{J_D}{2}\sum_{\langle\langle i,j\rangle\rangle}\left[e^{iJ_D\phi}\mathbf{S}^{+}_i\mathbf{S}^{-}_j+e^{-iJ_D\phi}\mathbf{S}^{-}_i\mathbf{S}^{+}_j+2J\mathbf{S}^{z}_i\mathbf{S}^{z}_j\right]%\nonumber\\
-h\sum_{i}\mathbf{S}^{z}_i,\nonumber\\
\end{eqnarray}
where $J=J_2/J_D$, $J_D=\sqrt{J_2^2+D^2}$, $\tan\phi=D/J_2$,  $\langle i,j\rangle$ denotes the NN bonds and $\langle\langle i,j\rangle\rangle$ denotes the NNN connections. $D$ is the strength of DM interaction, $D_{ij}=D\hat{\mathbf{k}}$. $\hat{\mathbf{k}}$ is the unit vector in the $z$ direction that is chosen as quantization axis. A representation of the LL lattice studied is made in Fig.~\ref{fig_120}, where each unit cell of the lattice presents three different sites. The lattice parameter (distance between the sites $AA$, $BB$) is given by the vectors ($a=1$): $\mathbf{a}_1=(a,0)$ and $\mathbf{a}_2=(0,a)$, that connects the sites $BA$ and $BC$, respectively. The arrows in blue represent: $SJ_De^{i\phi}$, in the arrow pointed on the direction $CA$, and $SJ_De^{-i\phi}$ pointed in the direction $AC$.

\begin{figure}
    \centering
    \begin{center}
\includegraphics[width=11.0cm]{f11.eps} \\
\end{center}
%\includegraphics[width=7.0cm]{EntagFerro2_Ferroquadrupolar33}\\
\caption{Scheme of the LL lattice.  Each unit cell presents three different sites. The lattice parameter is $a=1$: vectors $\mathbf{a}_1=(a,0)$ and $\mathbf{a}_2=(0,a)$ connects the sites $BA$ and $BC$, respectively. The arrows in blue corresponds to $SJ_De^{i\phi}$, in the arrow pointed on the direction $CA$ and $SJ_De^{-i\phi}$ pointed in the direction $AC$. }\label{fig_120} %$J_1=−1$, $J_2=1$, $J_{bq_1}=−3$, $J_{bq_2}=0.5$, (black),
\end{figure}

In the Holstein-Primakoff transformation, the spin operators are expressed in terms of creation and annihilator boson operators as
\begin{eqnarray}
\mathbf{S}_j^{+}=\sqrt{2S-\mathbf{n}_j} a_j,\hspace{0.5cm} \mathbf{S}_j^{-}=a_j^{\dag}\sqrt{2S-\mathbf{n}_j}, \hspace{0.5cm}\mathbf{S}_j^{z}=S-\mathbf{n}_j,\nonumber\\
\end{eqnarray}
 where $\mathbf{n}_j=a_j^{\dag}a_j$ is the density operator and we consider spin-$1/2$. In following, we expand the square roots up first order in powers of $1/S$
\begin{eqnarray}
\mathbf{S}_j^{+}\approx\sqrt{2S}a_j,\hspace{0.5cm} \mathbf{S}_j^{-}\approx\sqrt{2S}a_j^{\dag}, \hspace{0.5cm}\mathbf{S}_j^{z}\approx S.\nonumber\\
\end{eqnarray}
Consequently, the spin-wave Hamiltonian, neglecting the magnon-magnon interaction, is written as
\begin{equation}
\mathcal{H}_{SW}=\sum_{\mathbf{k}}a_{\mathbf{k}}^{\dag}H({\mathbf{k}})a_{\mathbf{k}},
\end{equation}
where
\begin{equation}
\hspace{-0.5cm}H({\mathbf{k}})=\left(
                 \begin{array}{ccc}
                   4J_1S & -2\cos(k_x)J_1S & -2\cos(k_y)J_1S \\
                   -2\cos(k_x)J_1S & 4J_1S+4J_2S & h({\mathbf{k}}) \\
                   -2\cos(k_x)J_1S & h({\mathbf{k}}) & 4J_1S+4J_2S \\
                 \end{array}
               \right).
\end{equation}
We perform the Fourier transform of the boson operators: $a_j=\frac{1}{\mathcal{N}}\sum_{\mathbf{k}}e^{i\mathbf{k}\cdot \mathbf{r}_j}a_{\mathbf{k}}$, where $\mathbf{r}_j$ is the relative distance of the $j$th site into the unit cell, ($j=1,2,3$), $\mathcal{N}$ is the number of unit cells. Furthermore,\\ $h({\mathbf{k}})=-4J_2S\cos(k_x)\cos(k_y)-4iDS\sin(k_x)\sin(k_y)$ and\\ $a_{\mathbf{k}}^{\dag}=\left(a_1^{\dag}(\mathbf{k}),a_2^{\dag}(\mathbf{k}),a_3^{\dag}(\mathbf{k})\right)$. Thus, we obtain three bands after to diagonalize  the Hamiltonian. For $J_2=D=0$, the spectrum of $H({\mathbf{k}})$ consists of one degenerate flat mode $\Lambda_1(\mathbf{k})=4J_1S$ and two dispersive modes $\Lambda_{1,3}(\mathbf{k})=4J_1S\mp2J_1S\sqrt{\cos^2(k_x)+\cos^2(k_y)}$, where the three bands if touch  at Dirac cone. If the NNN coupling is included but the DM interaction not, a gap will open between the two bands (lower bands) and the flat band becomes dispersive. When the DM interaction is turned on and $J_2=D$, a new gap will be opened between the upper bands while the gap between the lower bands will be closed, where they if touch on at the single Dirac cone. When $J_2\ne D$, all three bands are isolated and two gaps are opened.\cite{Cao}
A topological phase transition may happen with the tuning of $D$ and $J_2$. The Chern numbers of the three bands are topological invariant at topological phase transition, being obtained by the integral of the Berry curvature on the first Brillouin zone
\begin{equation}\label{berry}
  \Omega^m=\frac{1}{(2\pi)^2}\int\int_{BZ}\mathfrak{B}^{m}({\mathbf{k}})d^2k
\end{equation}
and
\begin{equation}
\mathfrak{B}^{m}({\mathbf{k}})=i\sum_{m'\ne m}\frac{\phi_m^{\dag}\frac{\partial \mathcal{H}_{SW}}{\partial k_x}\phi_{m'}\phi_{m'}^{\dag}\frac{\partial \mathcal{H}_{SW}}{\partial k_y}\phi_{m}-(k_x\leftrightarrow k_y)}{\left(\omega_m-\omega_{m'}\right)^2},
\end{equation}
where $\Lambda_m$ and $\phi_m$ are eigenvalues and eigenvectors of the $n$th band and there are four (different) topological phases given by the Chern numbers $(\Omega^1,\Omega^2,\Omega^3)$ of the three bands. Thus, the tuning of $J_2$ and $D$ indeed gives a topological phase transition.

\subsection{Tight-binding model on the Lieb lattice}%{\label{7}}

The two-dimensional LL lattice of non-interacting fermions in a square lattice on the Lieb lattice is characterized by three atoms per unit cell. Introducing the fermion operators $\left(a^{\dag}_{lm},b^{\dag}_{lm},c^{\dag}_{lm}\right)$, the Hamiltonian in a perpendicular magnetic field in the momentum space is given by\cite{Nija}
\begin{eqnarray}
\hspace{-0.0cm}\mathcal{H}=\sum_{\mathbf{k}}\left(
                                                                                                             \begin{array}{ccc}
                                                                                                               a^{\dag}({\mathbf{k}}) & b^{\dag}({\mathbf{k}}) & c^{\dag}({\mathbf{k}}) \\
                                                                                                             \end{array}
                                                                                                           \right)
\left(
                                                                                                             \begin{array}{ccc}
                                                                                                               0 &\mathfrak{t}^{*}_{x}(\mathbf{k})  &\mathfrak{t}^{*}_{y}(\mathbf{k})  \\
                                                                                                              \mathfrak{t}^{}_{x}(\mathbf{k})  & 0 & 0 \\
                                                                                                              \mathfrak{t}^{}_{y}(\mathbf{k})  & 0 & 0 \\
                                                                                                             \end{array}
                                                                                                           \right)\left(
                                                                                                                    \begin{array}{c}
                                                                                                                      a({\mathbf{k}}) \\
                                                                                                                      b({\mathbf{k}}) \\
                                                                                                                      c({\mathbf{k}}) \\
                                                                                                                    \end{array}
                                                                                                                  \right),\nonumber\\
\end{eqnarray}
where $\mathfrak{t}_x(\mathbf{k})=\mu_x(1+e^{ik_x})$ and $\mathfrak{t}_y(\mathbf{k})=\mu_y(1+e^{ik_y})$.  We perform the Fourier transform given by\\ $c_{\mathbf{k}}=\frac{1}{\sqrt{NM}}\sum_{l=1}^{M}\sum_{m=1}^{N}c_{lm}e^{i{\left(k_xl+k_ym\right)}}$, assuming that the lattice is composed of $\mathcal{N}=NM$ cells. We have the following eigenvalues of the model
\begin{equation}
\Lambda_{\pm}({\mathbf{k}})=\pm2\sqrt{\mu_x^2\cos^2\left(\frac{k_{\xi}}{2}\right)+\mu^2_y\cos^2\left(\frac{k_{\psi}}{2}\right)}, \hspace{0.5cm}\Lambda_{0}({\mathbf{k}})=0,
\end{equation}
where $\Lambda_{\pm}$ are the energies of the upper and lower bands, respectively and $\Lambda_0$ is the flat (non-dispersive) band of the LL lattice. The point $\Gamma=\left(\pi,\pi\right)$ is the  most interesting point into the Brillouin zone, where in the continuum limit or infinite lattice the three branches touch each other. The energies are given by to a massless spectrum $\Lambda_{\pm}({\mathbf{k}})=\pm\sqrt{\mu_x^2k_{x}^2+\mu_y^2k_y^2}$, where $k_{x}$ and $k_y$ are measured from the point $\Gamma$. The expansion of the same functions at $R=(0,0)$ give to a parabolic dependence $\Lambda_{\pm}({\mathbf{k}})=\pm\left(\frac{k_{x}^2}{2m_x}+\frac{k_y^2}{2m_y}\right)$, where $m_x$ and $m_y$ are the effective masses along the two directions. For $\mu_x=\mu_y=\mu$, we have effective mass exhibiting opposite signs generating a change of sign of the Hall effect.





\subsection{Concurrence}
In general, the bipartite entanglement is investigated using either negativity or other quantifiers as the negativity, concurrence and VN entropy. The concurrence is  expressed in terms of correlations functions of entanglement particles, where the bipartite entanglement in the model is made.
The concurrence of the qubits 1 and 2, $\varrho_{12}$ is given by the density matrix of either pure state or mixed state, being given by the square roots of the eigenvalues of $\mathfrak{R}_{1,2}=\rho_{12}(\sigma_z\otimes\sigma_z)\rho_{1,2}^{*}(\sigma_z\otimes\sigma_z)$:\cite{Hill,W}
%\begin{equation}
$C=\max\{\omega_1-\omega_2-\omega_3-\omega_4,0\}$,
%\end{equation}
where the quantities $\omega_1\geq\omega_2\geq\omega_3\geq\omega_4$ are the square roots of the eigenvalues of the operator $\mathfrak{R}_{1,2}$ and $\omega_1$ is the largest of the eigenvalues. This formula is applicable to the states of an arbitrary number of qubits.\cite{Hill,W} For the state $\varrho_{12}$,  a nonzero concurrence means that the qubits 1 and 2 are correlated or entangled and a zero concurrence corresponds to an unentangled state. The maximally entangled state corresponds to $C=1$. For $\mathcal{N}$ qubits, we apply the above transformation to each individual qubit.\cite{Vidal} \\

%

The concurrence relates with the ground state energy $U$ given by $C=1/2\max\left[0,-U/J\mathcal{N}-1\right]$, for a antiferromagnetic system and $C=1/2\max\left[0,U/3J\mathcal{N}-1\right]$ for the ones ferromagnetic.\cite{Vidal} Thus, the entanglement is determined by $\mathcal{Z}$, where $\mathcal{Z}|\langle\sigma_z\otimes\sigma_z\rangle|=|\sum_ie^{-\beta E_i}\langle i|\sigma_z\otimes\sigma_z|i\rangle|\leq|\sum_ie^{-\beta E_i}\langle i|\sigma_z\otimes\sigma_z|i\rangle\|\leq\sum_{i}e^{-\beta E_i}=1$, where we have $-1\leq U/3J\mathcal{N}\leq 1$. For the antiferromagnet, the increase of $U$ will generate a decreasing in the concurrence until it will become zero for a value as $-\mathcal{N}J$. As $U$ is given by
$U=T^2\frac{\partial}{\partial T}\log_2\mathcal{Z}$. Consequently the concurrence is given by
\begin{eqnarray}
\hspace{-0.75cm}C=-\frac{T^2}{3\mathcal{N}}\frac{\partial}{\partial T}\sum_{\mathbf{k}}\log_2\left(1+e^{-\beta\Lambda(\mathbf{k})/T}\right).
\end{eqnarray}
\begin{figure}
    \centering
\includegraphics[width=5.5cm]{tight-bind} \\
\caption{Von Neumann entropy as a function of $T$ for the dispersive bands $\Lambda_{+}(\mathbf{k})$ and for different values of couplings ($\mu_x,\mu_y$).}\label{fig_16}
\end{figure}

In Fig.~\ref{fig_15}, we obtain the concurrence $C$ as a function of the DM interaction $D$. We present the behavior of the concurrence
$C$ vs. $T$ using the spin wave approach for different dispersive bands $\Lambda_{1,3}(\mathbf{k})$. In general, $C$  decreases  for the two-dimensional Heisenberg
model with $T$ up to
dropping to zero at a threshold temperature $T_0$.\cite{galisova}
This trend of $C$ can be observed far from the
phase boundary at $T=0$.
We find a
sudden change in the concurrence with the $D$ and $J_2$ couplings. The quantum fluctuations are key measurable quantities that
can reflect on amount of entanglement in a
subsystem, where there is a link between entanglement
and bipartite fluctuations. For a non-interacting fermions system, the
entanglement can be experimentally measured through
conserved quantities such as
the particles number $\mathcal{N}$ when we are studying the
charges transfer through the system.
\begin{figure}
    \centering
\includegraphics[width=5.5cm]{conc} \\
\caption{Concurrence as a function of $T$ for different dispersive bands $\Lambda_1(\mathbf{k})$ (solid black line) that correspond to the set of values of Dzyaloshinskii-Moriya interaction strength and NNN coupling $D=0.0$, $J_2=0.0$  and band $\Lambda_3(\mathbf{k})$ (red-dashed line and blue-dotted line) that corresponds to the values $D=0.0$, $J_2=0.0$
and $D=0.2$, $J_2=0.2$ respectively. The bands touch each other at the Dirac cone for $D=0.0$, $J_2=0.0$, where a gap is opened in the two upper bands for $D\ne0.0$, $J_2\ne0.0$.}\label{fig_15}
\end{figure}


















\subsubsection{VN entropy for the antiferromagnet on the Lieb lattice}%{\label{res}}% by spin wave theory}{\label{res}}
In Fig.~\ref{fig_11}, we find the behavior of the VN entropy $E$ as a function of $T$. The behavior of $E$ is similar to thermodynamic entropy, where in the limit of high $T$ it gives only a qualitative description due to approach used.
 Furthermore, we obtain a split between the curves of the entanglement for the two magnon bands $\Lambda_1(\mathbf{k})$ (black solid-line and dashed-red-line) and $\Lambda_3(\mathbf{k})$ (green dot-dashed-line and blue dotted-line), where both lines tends if touch on in the $T\rightarrow 0$ limit. In the limit of Dirac cone in the first Brillouin zone, $J_2=0$, $D=0$, where the gap between the lower two bands is closed, the curves for $\Lambda_{1,3}(\mathbf{k})$ if becomes more damped however, without if touch on for high $T$.
In Fig.~\ref{fig_12}, we display $E$ as a  function of the  Dzyaloshinskii-Moriya interaction for a value of $J_2$ held fixed, such as $J_2=0.0$. At $D=0$, we are in the Dirac cone limit and we obtain a large discrepancy in  the two curves.  We find that $E$ decreases with the  increase of the strength of the couplings $D$ and $J_2$, obtaining the same behavior for $E$ vs. $J_2$, for $D$ held fixed in $D=0.0$, due to shape of Eq.~(\ref{entangl}).%Figs.~\ref{fig_13} and \ref{fig_14}.
\begin{figure}
    \centering
\includegraphics[width=5.5cm]{lieb.eps}\\
\caption{VN entropy as a function of $T$ for different dispersive bands $\Lambda_1(\mathbf{k})$ (solid black line and red-dashed line) that correspond to the set of values of Dzyaloshinskii-Moriya interaction strength, and NNN coupling $D=0.2$, $J_2=0.2$ and band $\Lambda_3(\mathbf{k})$ (green-dot-dashed line and blue dotted line) that corresponds to the values $D=0.0$, $J_2=0.0$. The bands if touch each other at  $D=0.0$, $J_2=0.0$, where a gap is opened in the two bands for $D\ne0.0$, $J_2\ne0.0$.}\label{fig_11}
\end{figure}







The nearest-neighbor Heisenberg antiferromagnet on the Lieb lattice is described by the Hamiltonian $\mathcal{H'}=\mathcal{H}+\mathcal{H}_{SI}+\mathcal{H}_{DM}$, where
\begin{equation}\label{model}
\mathcal{H}=\sum_{\langle i,j\rangle} J_{ij}\mathbf{S}_i\otimes \mathbf{S}_j,
%\mathcal{H}=J\sum_{\langle i,j\rangle}\left[\cos\theta (S_i\cdot S_j)+\sin\theta (S_i\cdot S_j)^2\right],
\end{equation}
where $\langle i,j\rangle$ denotes the NN bonds and $\theta_{ij}=\theta_i-\theta_j$ is the angle between two neighboring spins and
\begin{equation}
\otimes=\left(
  \begin{array}{ccc}
    \cos\theta_{ij} & 0 & -\sin\theta_{ij} \\
    0 & 1 & 0 \\
    \sin\theta_{ij} & 0 & \cos\theta_{ij} \\
  \end{array}
\right).
\end{equation}
The spins are assumed forming an ordered state with spins forming a coplanar  spin structure $120^{\circ}$. In Fig.~\ref{fig_120}, we display a representation of the lattice considered with the different spin configurations. The primitive vectors of the kagome lattice and numbering of site within each unit cell are shown in the Figure.
 In unities used here, we have $\hbar=c=k_B=1$, where $c$ is the light velocity $c$ and $k_B$ is the Boltzmann constant. Consequently, the dispersion relation energy $\xi_{\mathbf{k}}$ and temperature $T$ are given in $J_1$ units.
\begin{figure}
    \centering
\hspace{-0.3cm}\includegraphics[width=10.5cm]{f11.eps} \\
%\includegraphics[width=7.0cm]{EntagFerro2_Ferroquadrupolar33}\\
\caption{Different spin configurations, $\mathbf{k}=0$ (left side) and $\mathbf{k}=\sqrt{3}\times\sqrt{3}$ (right side). The primitive vectors and numbering of sites within the unit cell are shown in each figure.}\label{fig_120} %$J_1=−1$, $J_2=1$, $J_{bq_1}=−3$, $J_{bq_2}=0.5$, (black),
\end{figure}

We consider the extension of the Heisenberg model on the kagome lattice to the XXZ model with anisotropy of easy-plane, $0\leq\Delta\leq 1$, due to degeneracy among the $120^{\circ}$ coplanar states this model remains the same as the Heisenberg model for us to extend the parameter space and to study the effect of quantum fluctuations in the ground-state without degeneracy of the classical ground state.\cite{zhi} We introduce the Holstein-Primakoff representation given by
\begin{equation}
S_j^{-}=a_j^{\dag}\sqrt{2S-a_j^{\dag}a_j},\hspace{0.5cm}S_j^{+}=\left(S_j^{-}\right)^{\dag},\hspace{0.5cm}S_j^{z}=S-a_j^{\dag}a_j.
\end{equation}
Performing the Fourier transformation\\ $a_j=\frac{1}{\sqrt{N}}\sum_{\bf{k}}a_{\bf{k}}e^{i{\bf{k}}\cdot {\bf{r}}_j}$, the spin-wave dispersion relation for the XXZ model is given by
\begin{equation}
\varepsilon_{\bf{k}}=2JS\sqrt{1-\Delta\gamma_{\bf{k}}-\frac{1}{4}(1-\Delta)(1\pm\sqrt{1+8\gamma_{\bf{k}}})}
\end{equation}
for the dispersive modes and
\begin{equation}
\varepsilon_{\bf{k}}=2JS\sqrt{\frac{3}{2}(1-\Delta)},
\end{equation}
for the flat mode, which is at finite energy. $\gamma_{\bf{k}}\equiv c_1c_2c_3$, $\sum_{m=1}^{3}c_m^2=1+2\prod_{m=1}^{3}c_m$ and $c_m=\cos\left(\mathbf{k}\cdot\delta_m/2\right)$, with $\delta_1=(1,0)$, $\delta_2=\left(\frac{1}{2},\frac{\sqrt{3}}{2}\right)$, $\delta_3=\delta_2-\delta_1$ in units of interatomic distance.
\subsection{Single-ion anisotropy}
The easy-plane anisotropy can be also generated with addition of a positive single-ion anisotropy term
\begin{equation}
 \mathcal{H}_{SI}=D\sum_i\left(S_i^z\right)^2,
\end{equation}
where this term gives zero contribution to the classical energy and does not contribute
to the degeneracy lifting of $120^{\circ}$ and the dispersion relation is given by
\begin{equation}\label{single}
\varepsilon_{\bf{k}}=2JS\sqrt{1+D-\gamma_{\bf{k}}-\frac{D}{4}(1\pm\sqrt{1+8\gamma_{\bf{k}}})},
\end{equation}
which should be compared to the XXZ results.
\subsection{Dzyaloshinskii-Moriya interaction}
The antisymmetric Dzyaloshinskii-Moriya spin coupling is given by\cite{Dzyaloshinskii,Moriya}
\begin{equation}\label{DM interaction}
 \mathcal{H}_{DM}=\sum_{\langle ih\rangle}\mathbf{D}_{ij}\cdot(\mathbf{S}_{i}\wedge\mathbf{S}_{j}).
\end{equation}
We consider $\mathbf{D}_{ij}=(0,0,D_z)$ which is a consequence of bond gauge, in witch the two types of triangles are circumvented oppositely.\cite{Zhitomirsky} The dispersion relation for this case is given by
\begin{equation}
\varepsilon_{\bf{k}}=2JS\sqrt{1+d_m-\gamma_{\bf{k}}-\frac{d_m}{4}(1\pm\sqrt{1+8\gamma_{\bf{k}}})},
\end{equation}
which is similar to the single-ion results.

%\bibitem{Kartsev} Alexey Kartsev, Mathias Augustin, Richard F. L. Evans, Kostya S. Novoselov, Elton J. G. Santos, npj Computational Materials 6, (2020) 150.
%\bibitem{AS} A. Smerald, Nic Shannon, Phys. Rev. B 88, (2013) 184430.
%\bibitem{Hu} W-J. Hu, H-H Lai, S-S. Gong, R. Yu, E. Dagotto, Q. Si, Phys. Rev. Research 2, (2020) 023359.
%\bibitem{Lauchli} A. L\"{a}uchli, F. Mila, K. Penc, Phys. Rev. Lett. 97, (2006) 087205.
%\bibitem{Lauchli2} A. L\"{a}uchli, F. Mila, K. Penc, Phys. Rev. Lett. 97, (2006) 229901.
%\bibitem{Lauchli3} A. Joshi, M. Ma, F. Mila, D. N. Shi, F. C. Zhang,  Phys. Rev. B 60, (1999) 6584.
%\bibitem{lslima2017} L. S. Lima, J. Magn. Magn. Mater. 428, (2017) 448.
%\bibitem{lslima2021} L. S. Lima, J. Magn. Magn. Mater. 525, (2021) 167657.
%Where the different biquadratic couplings present in spin-1 antiferromagnets favor to the nematic ordering.\cite{Lauchli,Lauchli2,Lauchli3,lslima2017,lslima2021,lslima20212}%\cite{Lauchli,Lauchli2,Lauchli3,lslima2017,lslima20212,lslima2021}
%This compound is parent of the  iron-based superconductors or iron pnictides which exhibits a rich variety of unusual phases. Although the N\'eel antiferromagnetic order is found
%in the cuprates superconductors, the iron pnictides  display a  collinear antiferromagnetic
%order, as well however, in the cuprates superconductors the magnetism of parent compounds is well described
%by a two-dimensional nearest-neighbor Heisenberg model while in the iron pnictides this type
%of magnetic interaction is not well
%described by models like the Heisenberg model.

\textit{SU(3) flavor-wave theory (SBMF):} For the calculation of the von Neumann entropy Eq.~(\ref{entangl}), we use the SU(3) flavor-wave formalism (SBMF) proposed in Refs.~\cite{papanicolau,Wang} to study the XXZ model at phase of large anisotropy. An extension of the method to study the model Eq.~\ref{model} was performed in Refs.~\cite{Luo,piresbic}. Three boson operators $t_x$, $t_y$ and $t_z$ are defined acting in the vacuum state $\mid v\rangle$ as
\begin{eqnarray}
t_x^{\dag}|v\rangle=|x\rangle,\hspace{0.75cm}t_y^{\dag}|v\rangle=|y\rangle,\hspace{0.75cm}t_z^{\dag}|v\rangle=|z\rangle,
\end{eqnarray}
being this method adequate to study the nematic order. Moreover, in terms of operators $t_{\alpha,\beta,\gamma}$, ($\alpha,\beta,\gamma=x,y,z$),  the spin operators $S_j^{\alpha,\beta,\gamma}$  are written as
\begin{eqnarray}\label{SB}
S_j^{\alpha}  = -i(t_{\beta}^{\dag}t_{\gamma}-t_{\gamma}^{\dag}t_{\beta}).
\end{eqnarray}
\noindent
In following, we define other boson operators $u_{\mu}^{\dag}$;  $\mu=1,2$  by
\begin{equation}
u_{1}^{\dag}=\frac{1}{\sqrt{2}}(t_{x}^{\dag}+it_{y}),\hspace{0.75cm}u_{2}^{\dag}=\frac{1}{\sqrt{2}}(t_{x}^{\dag}-it_{y})\\
\end{equation}
with the additional constraint condition $\sum_{\mu}u_{\mu}^{\dag}u_{\mu}=1$.
\begin{figure}
    \centering
%\includegraphics[width=7.0cm]{EntagFerro3.eps} \\
%\includegraphics[width=7.0cm]{EntagFerro2_Ferroquadrupolar22}\\
\caption{von Neumann entropy as a function of the next-nearest-neighboring interaction $J_2>0$ (AFM), for the values $J_1$, $J_{bq_1}$ and $J_{bq_2}$ held fixed at points $\Gamma=(0,0)$, $M=(\pi,0)$, and
$X=(\pi,\pi)$.}\label{fig_10}
\end{figure}
The resulting Hamiltonian can be written in the diagonal form via the following Bogoliubov transformation
\begin{equation}
u_{\mu\mathbf{k}}=\chi_{\mu\mathbf{k}}\alpha_{\mu\mathbf{k}}-\rho_{\mu\mathbf{k}}\beta^{\dag}_{\mu\mathbf{k}},
\end{equation}
where $\chi_{\mathbf{k}}=\left(\Pi_{\mathbf{k}}+\xi_{\mathbf{k}}\right)/2\xi_{\mathbf{k}}$ and $\rho_{\mathbf{k}}=\left(\Pi_{\mathbf{k}}-\xi_{\mathbf{k}}\right)/2\xi_{\mathbf{k}}$,
$\mu=1,2$. Thus, the ground state is determined with the dispersion relation given by
\begin{equation}
\xi_{\mathbf{k}}=\sqrt{\Pi_{\mathbf{k}}^2-\Theta_{\mathbf{k}}^2},
\end{equation}
where
\begin{eqnarray}
\Pi_{\mathbf{k}}=\lambda + 4t^2\left( J_1\gamma_{\mathbf{k}}+ J_2\tilde{\gamma}_{\mathbf{k}}\right), \\ \Theta_{\mathbf{k}}=4t^2\left[(J_1-J_{bq_1})\gamma_{\mathbf{k}}+ (J_2-J_{bq_2})\tilde{\gamma}_{\mathbf{k}}\right],\\ \tilde{\gamma}_{\mathbf{k}}=\cos k_x\cos k_y.
\end{eqnarray}
Furthermore, the ground state energy is given by
  $\mathcal{E}_{\mathbf{k}}=\int_{BZ}d\mathbf{k}(\xi_{\mathbf{k}}-\Pi_{\mathbf{k}})+C$,
where $C$ is an additive constant. From minimization of $\mathcal{E}_{\mathbf{k}}$, we obtain the integral equation \\ $t^2=2-\int_{BZ}d\mathbf{k}\Pi_{\mathbf{k}}/\xi_{\mathbf{k}}$, where the integral is performed into the first Brillouin zone (BZ).
The ferroquadrupolar phase spontaneously breaks the spin rotational symmetry of the system where there is a gapless Golsdstone mode in the dispersion relation $\xi_{\mathbf{k}}$ at point $\mathbf{k} = (0,0)$.


\textit{Analysis by DMRG in two dimensions:} Density matrix renormalization group (DMRG) is a well known numerical technique suited to treat the one-dimensional spin-1/2 Heisenberg model \cite{white1,white2}. However, any finite two-dimensional lattice can be mapped in a one-dimensional lattice where the sites of the lattice are numbered and therefore long range interactions are introduced. Since in mean field theories, the dynamics of the operators $\mathcal{S}(\mathcal{r})$ is omitted, the variational principle  assumes an expectation value of the operator $\langle\vec{S}(\vec{r})\rangle$ where we neglect the fluctuations and drop higher-order terms.
In Fig.~\ref{fig_5}, we present the von Neumann entropy as a function of $J_{bq_2}$ using DMRG. The calculations were performed for different lattice sizes ($L=16$, $L=64$). We obtain an influence of $J_{bq_2}$  on quantum entanglement.  Thus, the analysis by DMRG in Fig.~\ref{fig_5} seems to confirm the decreasing of the concurrence $C$ with the coupling $J_{bq_2}$ and hence, the dependence of the critical line between the N\'eel and quadrupolar phases with the quantum fluctuations.%even though a quantitative agreement could not be reached. 