## 2 The Physical Model and Its Master Equation  
The incorporating Hamiltonian of the Heisenberg XYZ chain with Dyzaloshinskii-Moriya (DM) interaction is given by  
\[\hat{H} = \sum_{k=1}^{n}[J_{x}\hat{\sigma}_{k}^{x}\hat{\sigma}_{k+1}^{x}+J_ {y}\hat{\sigma}_{k}^{y}\hat{\sigma}_{k+1}^{y}+J_{z}\hat{\sigma}_{k}^{z}\hat{ \sigma}_{k+1}^{z}+\vec{D}.(\vec{\sigma}_{k}\times\vec{\sigma}_{k+1})], \tag{1}\]  
where \(J_{i}\) (\(i=x,y,z\)) are the symmetric spin-spin interaction coupling constants, \(\vec{\sigma}_{k}\) is the Pauli operators vector for any qubit, and \(\vec{D}\) is the DM vector. Despite its simplicity, the model describes the magnetism of several real physical systems based on the DM interaction which arises from spin-orbit coupling [24; 25]. Moreover, it was successful in explaining the ferromagnetic and antiferromagnetic behaviors. Besides, the models of the XYZ chain have been used to accomplish more important resources of the quantum information, like quantum coherence [26] and quantum correlations [27; 28; 29]. According to experimental demonstration [30], the system of spin-\(\frac{1}{2}\) pyrochlores R\({}_{2}\)V\({}_{2}\)O\({}_{7}\) can realize the DM interaction in the anisotropic XYZ chain. Here, we consider the case of the two-qubit of the spin chain with the DM interaction along \(z\)-axis. Thus, the Hamiltonian can be expressed as  
\[\hat{H} = J_{x}\hat{\sigma}_{1}^{x}\hat{\sigma}_{2}^{x}+J_{y}\hat{\sigma}_{ 1}^{y}\hat{\sigma}_{2}^{y}+J_{z}\hat{\sigma}_{1}^{z}\hat{\sigma}_{2}^{z}+ \chi_{z}(\hat{\sigma}_{1}^{x}\hat{\sigma}_{2}^{y}-\hat{\sigma}_{1}^{y}\hat{ \sigma}_{2}^{x}), \tag{2}\]  
where \(\chi_{z}\) represents the strength of DM interaction. This Hamiltonian is exactly diagonalizable and satisfies the eigenvalue problem: \(\hat{H}|\Psi_{k}\rangle=E_{k}|\Psi_{k}\rangle(k=1-4)\), with four eigenstates \(|\Psi_{k}\rangle\) corresponds to four eigenvalues \(E_{k}\). In the basis two-qubit states: \(\{|\varpi_{1}\rangle=|0_{A}0_{B}\rangle,\,|\varpi_{2}\rangle=|0_{A}1_{B}\rangle, \,|\varpi_{3}\rangle=|1_{A}0_{B}\rangle,\,|\varpi_{4}\rangle=|1_{A}1_{B}\rangle\}\), the eigenstates \(|\Psi_{k}\rangle\) of Equation (1) are given by \[\left( \begin{array}{c}|\Psi_{1}\rangle\\ |\Psi_{2}\rangle\\ |\Psi_{3}\rangle\\ |\Psi_{4}\rangle\end{array} \right)=M\left( \begin{array}{c}|\sigma_{1} \rangle\\ |\sigma_{2}\rangle\\ |\sigma_{3}\rangle\\ |\sigma_{4}\rangle\end{array} \right),\ M=\left( \begin{array}{cccc}\frac{1 }{\sqrt{2}}&0&0&\frac{1}{\sqrt{2}}\\ 0&\frac{1}{\sqrt{2}}&\frac{u}{\sqrt{2}}&0\\ 0&\frac{1}{\sqrt{2}}&-\frac{u}{\sqrt{2}}&0\\ \frac{1}{\sqrt{2}}&0&0&-\frac{1}{\sqrt{2}}\end{array} \right), \tag{3}\]  
and the corresponding eigenvalues are  
\[E_{1} = J_{x}-J_{y}+J_{z},\quad E_{2}=J_{z}+\sqrt{4\chi_{z}^{2}+S^{2}}, \tag{4}\] \[E_{3} = J_{z}-\sqrt{4\chi_{z}^{2}+S^{2}},\quad E_{4}=-J_{x}+J_{y}+J_{z},\]  
with \(\mu=(S+2i\chi_{z})/\sqrt{4\chi_{z}^{2}+S^{2}}=u+iv\) and \(S=J_{x}+J_{y}\). The principal purpose here is to investigate the robustness of the generated non-local correlations due to the unitary two-qubit interaction and against the effects of the phase damping, in the phase-damping model. In this model, (1) the phase is damped without energy relaxation [18; 31; 32] and (2) the reservoir is coupled to the two-qubit systems via \(\sigma_{i}^{+}=|1_{i}\rangle\langle 0_{i}|\) and \(\sigma_{i}^{-}=|0_{i}\rangle\langle 1_{i}|\) (\(i=A,B\)). The master equation of the model is given by  
\[\frac{d\hat{\rho}(t)}{dt} = -i[\hat{H},\hat{\rho}(t)]\] \[+ \frac{1}{2}\sum_{i=A,B}\gamma_{i}[2|1_{i}\rangle\langle 1_{i}| \hat{\rho}(t)|1_{i}\rangle\langle 1_{i}|-|1_{i}\rangle\langle 1_{i}|\hat{\rho}(t)- \hat{\rho}(t)|1_{i}\rangle\langle 1_{i}|],\]  
where \(\gamma_{i}\) are the phase damping rates, we take \(\gamma_{i}=\gamma\). To find analytical solutions for the master equation of Equation (5), we use the eigenstates representation (ESR) method, that is based on the eigenstates and the eigenvalues of the Hamiltonian of Equation (1).