Ground State, Magnetization Process and Bipartite Quantum Entanglement of a Spin-1/2 Ising–Heisenberg Model on Planar Lattices of Interconnected Trigonal Bipyramids

Abstract

The ground state, magnetization scenario and the local bipartite quantum entanglement of a mixed spin-$1/2$ Ising–Heisenberg model in a magnetic field on planar lattices formed by identical corner-sharing bipyramidal plaquettes is examined by combining the exact analytical concept of generalized decoration-iteration mapping transformations with Monte Carlo simulations utilizing the Metropolis algorithm. The ground-state phase diagram of the model involves six different phases, namely, the standard ferrimagnetic phase, fully saturated phase, two unique quantum ferrimagnetic phases, and two macroscopically degenerate quantum ferrimagnetic phases with two chiral degrees of freedom of the Heisenberg triangular clusters. The diversity of ground-state spin arrangement is manifested themselves in seven different magnetization scenarios with one, two or three fractional plateaus whose values are determined by the number of corner-sharing plaquettes. The low-temperature values of the concurrence demonstrate that the bipartite quantum entanglement of the Heisenberg spins in quantum ferrimagnetic phases is field independent, but twice as strong if the Heisenberg spin arrangement is unique as it is two-fold degenerate.

1. Introduction

Quantum entanglement has been attracting a lot of attention in the last few years mainly due to its crucial role in the development of quantum computers, superdense coding, quantum communication, quantum teleportation, as well as quantum information theory [1,2,3]. The application potential of this unique phenomenon also exceeds into the quantum biology [4,5] and quantum metrology [6,7].

In quantum theory, quantum entanglement provides a novel platform for exploring long-range quantum correlations, quantum phase transitions as well as exotic properties of many-body systems [8,9,10,11]. The low-dimensional Heisenberg spin models, involving quantum fluctuations between spins, play a significant role in this regard because they have been proven to be ideal candidates for a rigorous investigation of the entangled states under the influence of the external stimuli such as magnetic field (homogeneous or inhomogeneous) and/or temperature [12,13,14,15,16,17,18,19,20,21,22,23,24]. Moreover, many analytical and numerical calculations have been performed to examine the tuning of the quantum and thermal bipartite entanglement by varying the exchange anisotropy parameter [19,20,21,22,23,24,25,26,27,28], the uniaxial single-ion anisotropy [16,17], the Dzyaloshinskii–Moriya interaction (spin-orbit coupling) [18,19,20,26,27], the next-nearest-neighbour interaction [13,14,29], as well as by introducing impurities into the system [28,30].

However, the rigorous investigation of the bipartite entanglement in the pure Heisenberg models represents a complex task, which is considerably limited due to a non-commutability of spin operators in the Hamiltonian. This computational problem makes the rigorous study of the phenomenon in general inaccessible across whole parameter space of the systems. On the other hand, replacing some of the Heisenberg spins with three spin components by the Ising ones with only one (z-) component at the nodal lattice sites is the alternative way to exactly examine the entanglement in various simpler mixed-spin Ising–Heisenberg models by using the standard transfer-matrix method [31] and/or the concept of generalized mapping transformations [32,33,34,35]. Taking into account the fact that the finite Heisenberg clusters formed by three-component Heisenberg spins are indirectly coupled with each other through the intermediate one-component Ising spin(s), one finds that the eigenstates of two adjacent Heisenberg clusters are separable. Thus, any quantity measuring the local bipartite entanglement in the considered mixed-spin model can be rigorously calculated for each quantum Heisenberg cluster separately.

To date, the bipartite entanglement has been rigorously examined in several one- (1D) and two-dimensional (2D) mixed-spin Ising–Heisenberg models formed by the identical Heisenberg dimers or triangular clusters which interact with each other via the intermediate nodal Ising spin(s) [36,37,38,39,40,41,42,43,44,45]. The investigations brought a deeper insight into the thermal and magnetic-field-driven changes of the phenomenon [36,37,38,39,40,41,42,43], the impact of the model’s parameters on the phenomenon [36,39,40,41,42,43,44,45], as well as the evolution of the phenomenon near and above second-order (continuous) phase transitions [44,45] without any artefacts arising from approximations. Despite their simplicity and the general opinion that the simpler mixed-spin Ising–Heisenberg systems involving isolated local quantum correlations are artificial models, some of the results were in a very good correspondence with ones obtained for more complex Heisenberg counterparts [40] and also with experiments [38,46,47].

In the present paper, we will rigorously solve a spin-$1/2$ Ising–Heisenberg model in a longitudinal magnetic field on 2D lattices formed by identical corner-sharing trigonal bipyramidal plaquettes. Our recent studies [45,48] of the model without magnetic field on the particular lattice with four inter-connected bipyramidal units have shown that this quantum mixed-spin model represents a suitable playground for a rigorous study of various unconventional physical phenomena such as the macroscopic degeneracy of the spontaneous long-range order caused by chiral spin degrees of freedom, the spin frustration, and the bipartite entanglement. The aforementioned findings motivated us to extend the investigation of the model also to the effect of the longitudinal magnetic field. The goals of the present paper are to shed a light on the nature of ground states invoked by the applied field, to identify the actual fractional plateaus in the zero-temperature magnetization process, to find out a general formula describing how the values of these plateaus depend on the current number of interconnected bipyramidal plaquettes and, finally, to quantify the bipartite quantum entanglement between the Heisenberg spins in individual ground states.

In addition to the academic interest, our investigation of the spin-$1/2$ Ising-Heisenberg model on 2D lattices formed by interconnected trigonal bipyramids is motivated by the existence of a class of geometrically frustrated structures, namely cobaltates YBaCo${}_4$O${}_7$ (Y denotes a rare-earth ion) [49] and anion-radical salts (MDABCO${}^{+}$)(C${}_{60}^{•-}$) (MDABCO${}^{+}$ represents N-methyldiazabicyclooctanium cation, C${}_{60}^{•-}$ is a radical anion) [50], in which one can clearly identify corner-sharing trigonal bipyramidal clusters. Although the mentioned compounds do not represent a precise experimental realization of the magnetic structure proposed in the present paper, we hope that a targeted design of the magnetic material with a magnetic structure of interconnected trigonal bipyramids is feasible. The targeted chemical synthesis involving highly anisotropic spin carriers such Dy${}^{3+}$ or Co${}^{2+}$ magnetic ions and anion-radical salts could possibly afford desiring such a quantum mixed-spin system. The findings presented in this paper could serve as a motivation for chemists to achieve this goal.

The outline of the paper is as follows: in Section 2, a magnetic structure of the investigated model is described and the most important steps of its rigorous treatment combining the analytical and numerical approaches are clarified. In Section 3, we present the most interesting numerical results for the ground state and the magnetization process of the model. The section also includes an analysis of the bipartite quantum entanglement in the individual ground states. Finally, the summary of the most important findings are presented in Section 4.

2. Model and Its Rigorous Treatment

We consider a mixed spin-$1/2$ Ising–Heisenberg model in a longitudinal magnetic field on 2D lattices consisting of identical corner-sharing trigonal bipyramidal plaquettes, as is schematically depicted in Figure 1 for one particular lattice with four such plaquettes. In this figure, the common vertices of plaquettes (white circles) are occupied by the Ising spins $\sigma =1/2$ that interact with other spins solely through their z-components. The rest ones (red circles), forming internal equilateral triangles oriented perpendicularly to the plaquette axes, are occupied by the Heisenberg spins $S=1/2$ that are coupled to each other via x-, y-, and z-components. Assuming q bipyramidal plaquettes share a common vertex, the total Hamiltonian of the mixed spin-$1/2$ Ising–Heisenberg model can be written as a sum of plaquette (five-spin cluster) Hamiltonians $\widehat{\mathcal{H}}={\sum }_{j=1}^{Nq/2}{\widehat{\mathcal{H}}}_j$, where N labels the total number of the nodal lattice sites occupied by the Ising spins (we consider the thermodynamic limit $N\to \infty$). Each plaquette Hamiltonian ${\widehat{\mathcal{H}}}_j$ contains all exchange interactions realized within the jth Ising–Heisenberg trigonal bipyramid and Zeeman terms that describe the influence of the applied external magnetic field on magnetic moments of the individual spins:

$$\begin{array}{cc}\hfill {\widehat{\mathcal{H}}}_j=& -J_H\sum _{k=1}^3\left[\Delta ({\widehat{S}}_{j,k}^x{\widehat{S}}_{j,k+1}^x+{\widehat{S}}_{j,k}^y{\widehat{S}}_{j,k+1}^y)+{\widehat{S}}_{j,k}^z{\widehat{S}}_{j,k+1}^z\right]-J_I\sum _{k=1}^3{\widehat{S}}_{j,k}^z({\widehat{\sigma }}_j^z+{\widehat{\sigma }}_{j+1}^z)\hfill \\ \hfill & -H_H\sum _{k=1}^3{\widehat{S}}_{j,k}^z-\frac{H_I}{q}({\widehat{\sigma }}_j^z+{\widehat{\sigma }}_{j+1}^z).\hfill \end{array}$$

(1)

In the above, ${\widehat{S}}_{j,k}^{\alpha }$ ($\alpha =x,y,z$) and ${\widehat{\sigma }}_j^z$ are spatial components of the spin-$1/2$ operator of the Heisenberg spin from the jth triangle and z-component of the Pauli operator with the eigenvalues $\pm 1/2$ at the jth nodal lattice site, respectively, which satisfy the periodic boundary conditions ${\widehat{S}}_{j,4}^{\alpha }\equiv {\widehat{S}}_{j,1}^{\alpha }$ and ${\widehat{\sigma }}_{Nq/2+1}^z\equiv {\widehat{\sigma }}_1^z$. The parameter $J_H$ marks the XXZ Heisenberg interaction within the Heisenberg triangles, $\Delta$ is the exchange anisotropy parameter in this interaction, and $J_I$ labels the Ising-type interaction between the nearest-neighbouring Ising and Heisenberg spins. The last two terms $H_H$ and $H_I$ in the second line of Equation (1) are Zeeman terms, which account for the magnetostatic energy of the Heisenberg and Ising spins in an applied longitudinal magnetic field, respectively.

Figure 1. A schematic representation of the j-th trigonal bipyramidal plaquette and the mixed spin-$1/2$ Ising–Heisenberg model on the particular 2D lattice with four ($q=4$) corner-sharing bipyramidal plaquettes. White circles label lattice sites occupied by the Ising spin $\sigma =1/2$ and red circles denote lattice sites occupied by the Heisenberg spin $S=1/2$. Black dashed lines illustrate the Ising-type interaction $J_I$ between the Ising and Heisenberg spins and red solid lines indicate XXZ Heisenberg exchange interaction $J_H(\Delta )$ between the Heisenberg spins in the plaquette.

Figure 1. A schematic representation of the j-th trigonal bipyramidal plaquette and the mixed spin-$1/2$ Ising–Heisenberg model on the particular 2D lattice with four ($q=4$) corner-sharing bipyramidal plaquettes. White circles label lattice sites occupied by the Ising spin $\sigma =1/2$ and red circles denote lattice sites occupied by the Heisenberg spin $S=1/2$. Black dashed lines illustrate the Ising-type interaction $J_I$ between the Ising and Heisenberg spins and red solid lines indicate XXZ Heisenberg exchange interaction $J_H(\Delta )$ between the Heisenberg spins in the plaquette.

As we have shown in our recent work on the zero-field case of the model [45], it is convenient for further calculations to introduce the composite spin operators:

$$\begin{array}{c}\hfill {\widehat{\mathbf{t}}}_j=\sum _{k=1}^3{\widehat{\mathbf{S}}}_{j,k},{\widehat{t}}_j^{\alpha }=\sum _{k=1}^3{\widehat{S}}_{j,k}^{\alpha }(\alpha =x,y,z),\end{array}$$

(2)

which determine the total spin of the Heisenberg triangular clusters and its spatial components, respectively. From the definition of the latter operators, one can easily obtain the spin identity ${\left({\widehat{t}}_j^{\alpha }\right)}^2=3/4+2{\sum }_{k=1}^3{\widehat{S}}_{j,k}^{\alpha }{\widehat{S}}_{j,k+1}^{\alpha }$. This, in combination with the identity for the square of the total composite spin operator ${\left({\widehat{\mathbf{t}}}_j\right)}^2={\widehat{\mathbf{t}}}_j·{\widehat{\mathbf{t}}}_j={\left({\widehat{t}}_j^x\right)}^2+{\left({\widehat{t}}_j^y\right)}^2+{\left({\widehat{t}}_j^z\right)}^2$, allows one to find the following two relations for the Heisenberg spin operators from the same triangular cluster:

$$\begin{array}{c}\hfill \sum _{k=1}^3\left({\widehat{S}}_{j,k}^x{\widehat{S}}_{j,k+1}^x+{\widehat{S}}_{j,k}^y{\widehat{S}}_{j,k+1}^y\right)=\frac{1}{2}\left[{\left({\widehat{\mathbf{t}}}_j\right)}^2-{\left({\widehat{t}}_j^z\right)}^2\right]-\frac{3}{4},\sum _{k=1}^3{\widehat{S}}_{j,k}^z{\widehat{S}}_{j,k+1}^z=\frac{1}{2}{\left({\widehat{t}}_j^z\right)}^2-\frac{3}{8}.\end{array}$$

(3)

Bearing in mind the above relations and the definition of the z-component of the composite spin operator ${\widehat{t}}_j^z$ listed in Equation (2), the plaquette Hamiltonian (1) can be expressed in the alternative form:

$$\begin{array}{cc}\hfill {\widehat{\mathcal{H}}}_j=& \frac{3J_H}{8}(2\Delta +1)-\frac{J_H\Delta }{2}{\left({\widehat{\mathbf{t}}}_j\right)}^2+\frac{J_H}{2}(\Delta -1){\left({\widehat{t}}_j^z\right)}^2-J_I{\widehat{t}}_j^z({\widehat{\sigma }}_j^z+{\widehat{\sigma }}_{j+1}^z)\hfill \\ \hfill & -H_H{\widehat{t}}_j^z-\frac{H_I}{q}({\widehat{\sigma }}_j^z+{\widehat{\sigma }}_{j+1}^z).\hfill \end{array}$$

(4)

It is easy to prove that the operators ${\left({\widehat{\mathbf{t}}}_j\right)}^2$, ${\widehat{t}}_j^z$ appearing in Equation (4) satisfy the commutation relations $\left[{\widehat{\mathcal{H}}}_j,{\left({\widehat{\mathbf{t}}}_j\right)}^2\right]=0$ and $\left[{\widehat{\mathcal{H}}}_j,{\widehat{t}}_j^z\right]=0$, which implies that they are both conserved quantities with well defined quantum spin numbers $t_j(t_j+1)$ and $t_j^z=\{-t_j,-t_j+1,\dots ,t_j\}$ for $t_j=\{3/2,1/2\}$, respectively. In this regard, Equation (4) represents a fully diagonal form of the plaquette Hamiltonian (1), which implies that the corresponding energy eigenvalues can be expressed in terms of the respective quantum spin numbers:

$$\begin{array}{cc}\hfill E_{t_j,t_j^z}=& \frac{3J_H}{8}(2\Delta +1)-\frac{J_H\Delta }{2}{t}_j(t_j+1)+\frac{J_H}{2}(\Delta -1){\left(t_j^z\right)}^2-J_I{t}_j^z({\sigma }_j^z+{\sigma }_{j+1}^z)\hfill \\ \hfill & -H_H{t}_j^z-\frac{H_I}{q}({\sigma }_j^z+{\sigma }_{j+1}^z).\hfill \end{array}$$

(5)

At this calculation stage, the partition function of the considered spin-$1/2$ Ising–Heisenberg model can be partially factorized due to commuting character of different plaquette Hamiltonians and written in terms of the eigenvalues (5) of the plaquette Hamiltonian:

$$\mathcal{Z}=\mathrm{Tr}\exp \left(-\beta \widehat{\mathcal{H}}\right)=\sum _{\left\{{\sigma }_n\right\}}\prod _{j=1}^{Nq/2}{\mathrm{Tr}}_j\exp \left(-\beta {\widehat{\mathcal{H}}}_j\right)=\sum _{\left\{{\sigma }_n\right\}}\prod _{j=1}^{Nq/2}\sum _{t_j,t_j^z}{g}_{t_j}\exp \left(-\beta E_{t_j,t_j^z}\right).$$

(6)

Here, $\beta =1/\left(k_{\mathrm{B}}T\right)$ ($k_{\mathrm{B}}$ is the Boltzmann’s constant and T is the absolute temperature of the system), and the summation symbol ${\sum }_{\left\{{\sigma }_n\right\}}$ denotes a summation over all possible spin configurations of the Ising spins, the product symbol ${\prod }_{j=1}^{Nq/2}$ runs over all trigonal bipyramids, and the double summation symbol ${\sum }_{t_j,t_j^z}$ runs over all possible values of the quantum numbers $t_j$, $t_j^z$ of the composite spins. Finally, $g_{t_j}$ is the degeneracy factor, which takes the value 1 for the quantum number $t_j=3/2$ and 2 for the quantum number $t_j=1/2$. After performing double summations over $t_j$ and $t_j^z$, one gains the effective Boltzmann’s weight $w({\sigma }_j^z,{\sigma }_{j+1}^z)$, which depends only on the Ising spin states ${\sigma }_j^z$, ${\sigma }_{j+1}^z$, and, thus, it can be replaced by a simpler but equivalent expression using the generalized decoration-iteration mapping transformation [32,33,34,35]:

$$\begin{array}{cc}\hfill w({\sigma }_j^z,{\sigma }_{j+1}^z)=& \sum _{t_j,t_j^z}{g}_{t_j}\exp \left(-\beta E_{t_j,t_j^z}\right)=2\exp \left[\frac{\beta H_I}{q}({\sigma }_j+{\sigma }_{j+1})-\frac{\beta J_H}{4}\right]\hfill \\ \hfill & \times \{\left[\exp \left(\beta J_H\Delta \right)+2\exp \left(-\frac{\beta J_H\Delta }{2}\right)\right]cosh\left[\frac{\beta J_I}{2}({\sigma }_j^z+{\sigma }_{j+1}^z)+\frac{\beta H_H}{2}\right]\hfill \\ \hfill & +\exp \left(\beta J_H\right)cosh\left[\frac{3\beta J_I}{2}({\sigma }_j^z+{\sigma }_{j+1}^z)+\frac{3\beta H_H}{2}\right]\}\hfill \\ \hfill =& A\exp \left[\beta J_{\mathrm{eff}}{\sigma }_j^z{\sigma }_{j+1}^z+\frac{\beta H_{\mathrm{eff}}}{q}({\sigma }_j^z+{\sigma }_{j+1}^z)\right].\hfill \end{array}$$

(7)

The novel effective parameters A, $J_{\mathrm{eff}}$, and $H_{\mathrm{eff}}$ emerging on the right-hand side of Equation (7) are determined by ’self-consistency’ of the algebraic approach used:

$$A=\sqrt[4]{w_{+}{w}_{-}{w}_0^2},J_{\mathrm{eff}}=k_{\mathrm{B}}Tln\left(\frac{w_{+}{w}_{-}}{w_0^2}\right),H_{\mathrm{eff}}=\frac{k_{\mathrm{B}}Tq}{2}ln\left(\frac{w_{+}}{w_{-}}\right).$$

(8)

Here, $w_{\pm }=w(\pm 1/2,\pm 1/2)$ and $w_0=w(\pm 1/2,\mp 1/2)$. After substituting Equation (7) into Equation (6), one obtains the rigorous equivalence between the partition function $\mathcal{Z}$ of the spin-$1/2$ Ising–Heisenberg model given by the Hamiltonian (1) and the partition function ${\mathcal{Z}}_{IM}$ of the effective spin-$1/2$ Ising model on the corresponding q-coordinated 2D lattice given by the Hamiltonian ${\mathcal{H}}_{IM}=-J_{\mathrm{eff}}{\sum }_{\langle j,n\rangle }^{Nq/2}{\sigma }_j^z{\sigma }_n^z-H_{\mathrm{eff}}{\sum }_{j=1}^N{\sigma }_j^z$:

$$\mathcal{Z}(T,J_H,J_I,\Delta ,H_H,H_I)=A^{Nq/2}{\mathcal{Z}}_{IM}(T,J_{\mathrm{eff}},H_{\mathrm{eff}}).$$

(9)

The mapping relation (9) represents the crucial result of the rigorous solution of the considered 2D spin-$1/2$ Ising–Heisenberg model in an external magnetic field because of all important physical quantities clarifying a ground-state arrangement, magnetization process, and quantum bipartite entanglement between the Heisenberg spins, namely, the local magnetization $m_I=\langle {\widehat{\sigma }}_j^z\rangle$ and $m_H=\langle {\sum }_{k=1}^3{\widehat{S}}_{j,k}^z\rangle$ per nodal Ising spin and Heisenberg triangular cluster, respectively, the total magnetization m per bipyramidal plaquette, as well as the pair correlation functions $C_{HH}^{xx\left(yy\right)}=\langle {\widehat{S}}_{j,k}^x{\widehat{S}}_{j,k+1}^x\rangle =\langle {\widehat{S}}_{j,k}^y{\widehat{S}}_{j,k+1}^y\rangle$, $C_{HH}^{zz}=\langle {\widehat{S}}_{j,k}^z{\widehat{S}}_{j,k+1}^z\rangle$ and $C_{IH}^{zz}=\langle {\widehat{\sigma }}_j^z{\widehat{S}}_{j,k}^z\rangle =\langle {\widehat{\sigma }}_{j+1}^z{\widehat{S}}_{j,k}^z\rangle$, can be directly derived from the formula for the Gibbs free energy $\mathcal{G}=-k_{\mathrm{B}}Tln\mathcal{Z}$ by means of the differential calculus:

$$\begin{array}{ccc}\hfill m_I& =& -\frac{1}{2N}\frac{\partial \mathcal{G}}{\partial H_I},m_H=-\frac{2}{Nq}\frac{\partial \mathcal{G}}{\partial H_H},m=\frac{q{m}_H+2m_I}{q},\hfill \end{array}$$

(10a)

$$\begin{array}{ccc}\hfill C_{HH}^{zz}& =& \frac{2}{3Nq{J}_H}\left(\Delta \frac{\partial \mathcal{G}}{\partial \Delta }-\frac{\partial \mathcal{G}}{\partial J_H}\right),C_{HH}^{xx\left(yy\right)}=-\frac{1}{3Nq{J}_H}\frac{\partial \mathcal{G}}{\partial \Delta },C_{IH}^{zz}=-\frac{1}{3Nq}\frac{\partial \mathcal{G}}{\partial J_I}.\hfill \end{array}$$

(10b)

The final analytical expressions of all the physical quantities listed in Equation (10a,b) depend on the on-site magnetization $m_{IM}={\langle {\sigma }_j^z\rangle }_{IM}$ and the pair correlation function $C_{IM}^{zz}={\langle {\sigma }_j^z{\sigma }_{j+1}^z\rangle }_{IM}$ of the effective 2D q-coordinated spin-$1/2$ Ising lattice with the temperature-dependent nearest-neighbour interaction $J_{\mathrm{eff}}$ in the temperature-dependent magnetic field $H_{\mathrm{eff}}$. Because an exact solution for the 2D spin-$1/2$ Ising model in an external magnetic field still belongs to unresolved issues of condensed matter physics, one has to resort to some numerical algorithm applicable to the 2D Ising lattices to gain the accurate results for $m_{IM}$ and $C_{IM}^{zz}$. In the present paper, we will employ the classical Monte Carlo (MC) simulations implementing the standard Metropolis algorithm [51,52] for the effective spin-$1/2$ Ising lattice of a sufficiently large linear size L.

3. Discussion of the Numerical Results

In this section, we will proceed to a discussion of the most interesting numerical results for the 2D spin-$1/2$ Ising–Heisenberg model in an external magnetic field with the antiferromagnetic Ising-type interaction $J_I<0$ between the Ising and Heisenberg spins. For simplicity, we will assume that the local magnetic fields acting in the Ising and Heisenberg spins are identical $H_I=H_H=H$. The absolute value of the interaction $J_I$ will be used as an energy unit for defining a relative strength of the Heisenberg interaction $J_H/\left|J_I\right|$ and the magnetic field $H/|J_I|$.

3.1. Ground-State Phase Diagrams

First, we take a look at possible magnetic ground-state arrangement of the model, which can be determined by a systematic inspection of the eigenvalues (5) of the plaquette Hamiltonian (1) for all possible combinations of quantum spin numbers $t_j$, $t_j^z$ entering therein. The typical ground-state phase diagrams are depicted in Figure 2 in the $J_H/|J_I|-H/|J_I|$ parameter plane for two representative values of the exchange anisotropy $\Delta =0.5$ and 2 by assuming four different numbers q of corner-sharing trigonal bipyramidal plaquettes forming 2D lattices. As one can see from Figure 2a, the ground-state phase diagram of the model with the easy-axis exchange anisotropy $\Delta =0.5$ contains four different ground states. Specifically, two ground states are macroscopically degenerate quantum ferrimagnetic phases ${|-1/2;1/2\rangle }_{R,L}$ and ${|1/2;1/2\rangle }_{R,L}$, which differ from each other only by the orientation of Ising spins with respect to the applied magnetic field as indicated by the corresponding eigenvectors and eigenenergies per plaquette:

$$\begin{array}{ccc}\hfill {|\pm 1/2;1/2\rangle }_{R,L}& =& \prod _{j=1}^{Nq/2}{\left|\pm \right.〉}_{{\sigma }_j^z}\otimes {\left|1/2,R\mathrm{or}L\right.〉}_{{▵}_j},\hfill \end{array}$$

(11a)

$$\begin{array}{ccc}\hfill E_{{|\pm 1/2;1/2\rangle }_{R,L}}& =& \frac{J_H}{4}+\frac{J_H\Delta }{2}\mp \frac{J_I}{2}-\frac{(q\pm 2)H}{2q}.\hfill \end{array}$$

(11b)

The state vector ${\left|1/2,R\mathrm{or}L\right.〉}_{{▵}_j}$ in Equation (11a) describes a quantum superposition of three different up-up-down spin states of the j-th Heisenberg triangular cluster with two opposite (Right- and Left-hand side) chiral degrees of freedom:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\left|1/2,R\right.〉}_{{▵}_j}& =\frac{1}{\sqrt{3}}{\left(\left|\uparrow \uparrow \downarrow \right.〉+{\mathrm{e}}^{\frac{2\pi \mathrm{i}}{3}}\left|\uparrow \downarrow \uparrow \right.〉+{\mathrm{e}}^{\frac{4\pi \mathrm{i}}{3}}\left|\downarrow \uparrow \uparrow \right.〉\right)}_{{▵}_j},\hfill \\ \hfill {\left|1/2,L\right.〉}_{{▵}_j}& =\frac{1}{\sqrt{3}}{\left(\left|\uparrow \uparrow \downarrow \right.〉+{\mathrm{e}}^{\frac{4\pi \mathrm{i}}{3}}\left|\uparrow \downarrow \uparrow \right.〉+{\mathrm{e}}^{\frac{2\pi \mathrm{i}}{3}}\left|\downarrow \uparrow \uparrow \right.〉\right)}_{{▵}_j},\hfill \end{array}\end{array}$$

(12)

The two-fold degeneracy of each Heisenberg triangle results in the field-independent macroscopic degeneracy $2^{Nq/2}$ of the phases ${|-1/2;1/2\rangle }_{R,L}$ and ${|1/2;1/2\rangle }_{R,L}$, which is obviously highly sensitive to the current number q of interconnected trigonal bipyramidal plaquettes (Heisenberg triangular clusters). The direct relation between the number of plaquettes sharing a common vertex and the macroscopic degeneracy of the phases ${|-1/2;1/2\rangle }_{R,L}$, ${|1/2;1/2\rangle }_{R,L}$ is also reflected in a current value of the residual entropy per nodal Ising spin observed in both the phases. Specifically, it proportionally grows with q [53]:

$$\begin{array}{c}\hfill \frac{{\mathcal{S}}_{\mathrm{res}}}{N{k}_{\mathrm{B}}}=\underset{N\to \infty }{lim}\frac{1}{N}ln{2}^{Nq/2}\approx 0.347q.\end{array}$$

(13)

The other two ground states are the classical ferrimagnetic phase $|-1/2;3/2\rangle$ and the fully saturated phase $|1/2;3/2\rangle$. These two phases again differ from each other only by the orientation of the Ising spins with respect to the applied magnetic field, while the Heisenberg spins are fully polarized into the magnetic field direction without any quantum correlations between their x- and y-components in both phases:

$$\begin{array}{ccc}\hfill |\pm 1/2;3/2\rangle & =& \prod _{j=1}^{Nq/2}{\left|\pm \right.〉}_{{\sigma }_j^z}\otimes {\left|\uparrow \uparrow \uparrow \right.〉}_{{▵}_j},\hfill \end{array}$$

(14a)

$$\begin{array}{ccc}\hfill E_{|\pm 1/2;3/2\rangle }& =& -\frac{3J_H}{4}\mp \frac{3J_I}{2}-\frac{(3q\pm 2)H}{2q}.\hfill \end{array}$$

(14b)

The uniqueness of the classical spin arrangements in the phases $|-1/2;3/2\rangle$ and $|1/2;3/2\rangle$ given by Equation (14a) is reflected in the zero entropy per Ising spin $\mathcal{S}/\left(N{k}_{\mathrm{B}}\right)=0$ in parameter regions corresponding to these phases.

It is obvious from Figure 2a that the classical ferrimagnetic phase $|-1/2;3/2\rangle$ can be detected in the whole parameter region with the ferromagnetic Heisenberg coupling $J_H/\left|J_I\right|>0$ and partially also in the region with the antiferromagnetic Heisenberg interaction $J_H/\left|J_I\right|<0$. By contrast, the macroscopically degenerate quantum phases ${|-1/2;1/2\rangle }_{L,R}$ and ${|1/2;1/2\rangle }_{L,R}$ are stable solely for the antiferromagnetic Heisenberg couplings $J_H/\left|J_I\right|<0$. Finally, the saturated phase $|1/2;3/2\rangle$ represents the actual ground state at high enough magnetic fields regardless of whether the ferro- or antiferromagnetic Heisenberg interaction $J_H/\left|J_I\right|$ is considered.

On the other hand, the ground-state phase diagram corresponding to the model with the easy-plane anisotropy $\Delta =2$ is a little more complex (see Figure 2b). It contains two more ground states in addition to the previous four, namely, the unique quantum ferrimagnetic phases $|-1/2;1/2\rangle$ and $|1/2;1/2\rangle$ with the Heisenberg triangular clusters in a symmetric quantum superposition of three possible up-up-down spin states but an opposite orientation of the Ising spins:

$$\begin{array}{ccc}\hfill |\pm 1/2;1/2\rangle & =& \prod _{j=1}^{Nq/2}{\left|\pm \right.〉}_{{\sigma }_j^z}\otimes \frac{1}{\sqrt{3}}{\left(\left|\uparrow \uparrow \downarrow \right.〉+\left|\uparrow \downarrow \uparrow \right.〉+\left|\downarrow \uparrow \uparrow \right.〉\right)}_{{▵}_j},\hfill \end{array}$$

(15a)

$$\begin{array}{ccc}\hfill E_{|\pm 1/2;1/2\rangle }& =& \frac{J_H}{4}-J_H\Delta \mp \frac{J_I}{2}-\frac{(q\pm 2)H}{2q}.\hfill \end{array}$$

(15b)

As shown in Figure 2b, both the quantum ferrimagnetic phases $|-1/2;1/2\rangle$ and $|1/2;1/2\rangle$ emerge in the ground-state phase diagram exclusively in the parameter region of the ferromagnetic Heisenberg interaction $J_H/\left|J_I\right|>0$. Naturally, the zero entropy per Ising spin $\mathcal{S}/\left(N{k}_{\mathrm{B}}\right)=0$ solely can be detected in their stability regions due to unique quantum spin arrangement given by Equation (15a).

In addition to their location in the zero-temperature $J_H/|J_I|-H/|J_I|$ parameter plane, it is also possible to understand from Figure 2 how the individual phases develop depending on the number q of corner-sharing bipyramidal plaquettes. Namely, the quantum ferrimagnetic phases ${|-1/2;1/2\rangle }_{R,L}$, $|-1/2;1/2\rangle$ and the classical one $|-1/2;3/2\rangle$ are gradually extended to stronger magnetic fields with increasing number q of the corner-haring plaquettes. Moreover, the classical phase $|-1/2;3/2\rangle$ simultaneously spreads to the regions of stronger antiferromagnetic (ferromagnetic) Heisenberg interactions $J_H/\left|J_I\right|<0$ ($J_H/\left|J_I\right|>0$). The remaining three phases ${|1/2;1/2\rangle }_{R,L}$, $|1/2;1/2\rangle$ and $|1/2;3/2\rangle$ faithfully follow the evolution of the adjacent ones ${|-1/2;1/2\rangle }_{R,L}$, $|-1/2;1/2\rangle$, $|-1/2;3/2\rangle$: the quantum phases ${|1/2;1/2\rangle }_{R,L}$, $|1/2;1/2\rangle$ are gradually shifted to stronger antiferromagnetic and ferromagnetic Heisenberg interactions $J_H/\left|J_I\right|<0$ and $J_H/\left|J_I\right|>0$, respectively, and the saturated one $|1/2;3/2\rangle$ is shifted to stronger magnetic fields.

3.2. Magnetization Process

The rich ground-state phase diagrams depicted in Figure 2 suggest various zero-temperature magnetization scenarios of the studied spin-$1/2$ Ising–Heisenberg model either with one, two or three different plateaus at fractional values of the saturation magnetization $m_{sat}=(3q+2)/\left(2q\right)$. In accordance with the definition of the total magnetization m per plaquette listed in Equation (10a), the values of these plateaus are given by a current number q of the trigonal bipyramids sharing a common vertex:

$$\frac{m}{m_{sat}}=\frac{q-2}{3q+2},\frac{q+2}{3q+2},\frac{3q-2}{3q+2}.$$

(16)

The first (lowest) magnetization plateau at $m/m_{sat}=(q-2)/(3q+2)$ can be identified at low magnetic fields in the stability regions of the macroscopically degenerate quantum ferrimagnetic phase ${|-1/2;1/2\rangle }_{R,L}$ and the unique quantum ferrimagnetic phase $|-1/2;1/2\rangle$. The second one at $m/m_{sat}=(q+2)/(3q+2)$ is a result of the spin arrangements present in the quantum phases ${|1/2;1/2\rangle }_{R,L}$ and $|1/2;1/2\rangle$, and therefore it can be found at moderate magnetic fields. The highest fractional plateau at $m/m_{sat}=(3q-2)/(3q+2)$ relates to the classical ferrimagnetic phase $|-1/2;3/2\rangle$.

To illustrate the above statements, two three-dimensional (3D) plots of the isothermal magnetization curves for the particular version of the lattice with four ($q=4$) corner-sharing bipyramidal plaquettes are depicted in Figure 3 at the sufficiently low temperature $k_{\mathrm{B}}T/\left|J_I\right|=1.5\times {10}^{-3}$. The plots are the outcomes of MC simulations for $100\times 100$ nodal lattice sites (Ising spins), which corresponds to 19,800 corner-sharing trigonal bipyramidal plaquettes. The adequate numerical accuracy was achieved by $12\times {10}^4$ MC steps per node. For easy reference, the interaction ratio $J_H/\left|J_I\right|$ is fixed to the same ranges and the anisotropy parameter $\Delta$ to the same values as were used in Figure 2. It can be understood from a comparison of these plots with corresponding ground-state phase diagrams in Figure 2 that the displayed low-temperature magnetization curves faithfully reflect up to seven different types of zero-temperature magnetization scenarios with the real $1/7$-, $3/7$-, and/or $5/7$-plateaus satisfying the general formulas listed in Equation (16):

i.

  ${|-1/2;1/2\rangle }_{R,L}{-|1/2;1/2\rangle }_{R,L}-|1/2;3/2\rangle$,

ii.

 ${|-1/2;1/2\rangle }_{R,L}{-|1/2;1/2\rangle }_{R,L}-|-1/2;3/2\rangle -|1/2;3/2\rangle$,

iii.

${|-1/2;1/2\rangle }_{R,L}-|-1/2;3/2\rangle -|1/2;3/2\rangle$,

iv.

 $|-1/2;3/2\rangle -|1/2;3/2\rangle$,

v.

  $|-1/2;1/2\rangle -|-1/2;3/2\rangle -|1/2;3/2\rangle$,

vi.

 $|-1/2;1/2\rangle -|1/2;1/2\rangle -|-1/2;3/2\rangle -|1/2;3/2\rangle$,

vii.

$|-1/2;1/2\rangle -|1/2;1/2\rangle -|1/2;3/2\rangle$

(see the magnetization curves of different colors). Steep continuous rises between different fractional plateaus as well as between fractional plateaus and the saturation magnetization indicate a presence of the discontinuous magnetization jumps that exist at the critical fields corresponding to the first-order phase transitions only at zero temperature. In agreement with the ground-state analysis performed in Section 3.1, the first three magnetization processes i.–iii., which contain the macroscopically degenerate quantum phases ${|\pm 1/2;1/2\rangle }_{R,L}$, can be observed for the easy-axis exchange anisotropy $\Delta =0.5$ and also the easy-plane exchange anisotropy $\Delta =2$, but only in the parameter region of the antiferromagnetic Heisenberg interactions $J_H/\left|J_I\right|<0$. On the other hand, the magnetization scenario iv., reflecting the single field-induced transition from the classical ferrimagnetic phase $|-1/2;3/2\rangle$ to the saturated one $|1/2;3/2\rangle$, appears for both the antiferromagnetic ($J_H/\left|J_I\right|<0$) and ferromagnetic ($J_H/\left|J_I\right|>0$) Heisenberg couplings. The last three magnetization processes v.–vii., which involve unique quantum ferrimagnetic phases $|\pm 1/2;1/2\rangle$, emerge in the parameter region of the ferromagnetic Heisenberg couplings $J_H/\left|J_I\right|>0$ under the condition $\Delta >1$. It should be noted for completeness that the steep staircase dependences of all magnetization curves plotted in Figure 3 are generally gradually smoothing upon increasing of temperature due to a thermal activation of excited states, until they completely disappear.

3.3. Quantum Bipartite Entanglement

The discussion in the last subsection will be devoted to a bipartite quantum entanglement of the Heisenberg spins in the individual ground states. It is obvious from the plaquette Hamiltonian (1) that the spins may be quantum-mechanically entangled only within the Heisenberg triangular clusters in individual plaquettes. Those from different plaquettes can never be entangled due to the Ising spin at their common vertices.

In general, a degree of the bipartite quantum entanglement between the Heisenberg spins at k-th and $(k+1)$-th vertex of the j-th plaquette can be quantified by the quantity referred to as concurrence [54]. For the studied 2D spin-$1/2$ Ising–Heisenberg model, the concurrence can be simply calculated from the local magnetization $m_H$ of the Heisenberg triangular cluster and the corresponding pair correlation functions $C_{HH}^{xx\left(yy\right)}$, $C_{HH}^{zz}$ defined by Equation (10a,b) through the following formula [55,56]:

$${\mathcal{C}}_{k,k+1}=\max \left\{0,4|C_{HH}^{xx\left(yy\right)}|-2\sqrt{{\left(\frac{1}{4}+C_{HH}^{zz}\right)}^2-{\left(\frac{m_H}{3}\right)}^2}\right\}.$$

(17)

Of course, the identical XXZ exchange coupling $J_H(\Delta )$ in a given Heisenberg triangle results in the same degree of the bipartite entanglement of the spin pairs. This is reflected in the same values of the corresponding concurrences:

$${\mathcal{C}}_{1,2}={\mathcal{C}}_{2,3}={\mathcal{C}}_{3,1}=\mathcal{C}.$$

(18)

The global picture on a degree of the bipartite quantum entanglement between the Heisenberg spin pairs in the individual ground-state phases is illustrated in Figure 4, which shows the low-temperature ($k_{\mathrm{B}}T/\left|J\right|=1.5\times {10}^{-3}$) density map of the concurrence $\mathcal{C}$ of the spin-$1/2$ Ising–Heisenberg model on the regular 2D lattice with four corner-sharing bipyramidal plaquettes in the $J_H/|J_I|-H/|J_I|$ plane for the fixed value of the exchange anisotropy parameter $\Delta =2$. The plotted data have again been obtained by MC simulations performed for the lattice of 19,800 corner-sharing bipyramidal plaquettes, whereas $2\times {10}^7$ MC steps per nodal Ising spin were used to achieve the accuracy better than ${10}^{-3}$. It is clear from a direct comparison of Figure 4 with the corresponding ground-state phase diagram depicted in Figure 2b that the Heisenberg spins forming triangular clusters are quantum-mechanically entangled only if the macroscopically degenerate phases ${|\pm 1/2;1/2\rangle }_{R,L}$ and the unique phases $|\pm 1/2;1/2\rangle$ are ground states. Due to the macroscopic degeneracy caused by two possible chiral degrees of freedom of each triangular cluster, the bipartite entanglement of the Heisenberg spins in former two phases is half weaker than that in latter ones. This is also proven by the corresponding zero-temperature asymptotic values of the concurrence ${\mathcal{C}}_{{|\pm 1/2;1/2\rangle }_{R,L}}=1/3$ and ${\mathcal{C}}_{|\pm 1/2;1/2\rangle }=2/3$. On the other hand, the remaining white region with the zero concurrence ${\mathcal{C}}_{|\pm 1/2;3/2\rangle }=0$ confirms completely non-entangled arrangements of the Heisenberg spins in the classical ferrimagnetic phase $|-1/2;3/2\rangle$ and the fully saturated phase $|1/2;3/2\rangle$.

To get a deeper insight onto a role of pair correlations between the Heisenberg spins in their bipartite quantum entanglement, the concurrence $\mathcal{C}$ as function of the magnetic field $H/|J_I|$ and the corresponding dependencies of the pair correlation functions $C_{HH}^{xx\left(yy\right)}$, $C_{HH}^{zz}$ are plotted in Figure 5 for the anisotropy parameter $\Delta =2$ and two selected interaction ratios $J_H/\left|J_I\right|=-4$ and 6 at the temperature $k_{\mathrm{B}}T/\left|J\right|=1.5\times {10}^{-3}$. The variations are completed by low-temperature dependences of the local magnetization $m_I$ and $m_H$ to facilitate identification of the current ground-state spin arrangement. We note for completeness that all the curves are results of the MC simulations for the lattice of $100\times 100$ nodal Ising spins by using $2\times {10}^7$ MC steps per node to achieve accuracy better than ${10}^{-3}$.

Figure 5a captures the sequence of field-induced phase transitions ${|-1/2;1/2\rangle }_{R,L}{-|1/2;1/2\rangle }_{R,L}-|1/2;3/2\rangle$. Evidently, the nonzero concurrence $\mathcal{C}=1/3$, which can be found at the magnetic fields $H/|J_I|<9$ due to stability of the macroscopically degenerate quantum phases ${|\pm 1/2;1/2\rangle }_{R,L}$, is a result of the negative pair correlation functions $C_{HH}^{xx\left(yy\right)}=C_{HH}^{zz}=-1/12$ and the reduced local magnetization $m_H=1/2$. Identical values of the transverse and longitudinal correlation functions and their minus sign clearly indicate that the macroscopically degenerate unsaturated bipartite entanglement of the Heisenberg spins from the same XXZ triangle comes from antiferromagnetic $xx\left(yy\right)$ correlations of these spins, which are of the same strength and character as those in z-axis direction.

A different situation can be found in Figure 5b, which illustrates the sequence of field-induced phase transitions $|-1/2;1/2\rangle -|1/2;1/2\rangle -|1/2;3/2\rangle$. Here, the low-temperature concurrence $\mathcal{C}=2/3$, which can be observed at the magnetic fields $H/|J_I|<7$ due to the presence of the unique quantum phases $|\pm 1/2;1/2\rangle$, comes from the positive value of the transverse pair correlation function(s) $C_{HH}^{xx\left(yy\right)}=1/6$, the negative longitudinal correlation function $C_{HH}^{zz}=-1/12$, and the reduced local magnetization $m_H=1/2$. It is easy to understand from the values of $C_{HH}^{xx\left(yy\right)}$ and $C_{HH}^{zz}$ that the origin of the local quantum bipartite entanglement of the Heisenberg spins peculiar to the unique quantum phases $|\pm 1/2;1/2\rangle$ lies in ferromagnetic $xx\left(yy\right)$ correlations and these are twice as strong as the antiferromagnetic ones along z-axis.

4. Conclusions

In the present work, we have comprehensively studied the ground-state properties, possible magnetization scenarios, and the local bipartite quantum entanglement of the Heisenberg spins in the individual quantum ground states of the mixed spin-$1/2$ Ising–Heisenberg model in a longitudinal magnetic field on 2D lattices formed by identical corner-sharing bipyramidal plaquettes. The numerical results have been obtained by combining the exact analytical approach called the decoration-iteration mapping transformation [32,33,34,35] with numerical Monte Carlo simulations utilizing the Metropolis algorithm [51,52].

It has been demonstrated that the ground-state phase diagram of the investigated quantum mixed-spin model qualitatively does not depend on its lattice topology (the number q of corner-sharing plaquettes). In general, it involves in total six different ground states, namely the standard ferrimagnetic phase, fully saturated phase, two unique quantum ferrimagnetic phases and two macroscopically degenerate quantum ferrimagnetic phases with two chiral degrees of freedom of the Heisenberg spins forming triangular clusters in plaquettes. It is also proven that the diversity of the ground-state phase diagram gives rise to seven different magnetization scenarios with one, two or up to three fractional plateaus. The magnitudes of these plateaus are determined by the current number q of the corner-sharing plaquettes.

Other interesting findings are concerned with the bipartite quantum entanglement, which has been quantified by the concept of the concurrence. We have verified that the Heisenberg spins of the same XXZ triangular cluster of a given plaquette can be entangled either due to stability of the unique quantum ferrimagnetic phases, where they are in a symmetric quantum superposition of three possible up-up-down states, or due to macroscopically degenerate quantum ferrimagnetic phases characterized by two chiral degrees of freedom of each Heisenberg triangle. The strength of the entanglement in all the phases does not depend on the applied magnetic field. Moreover, the corresponding values of concurrence clearly indicate that the entanglement of the Heisenberg spins is twice as strong when their arrangement is unique that when it is two-fold degenerate. Thus, it can be concluded that the macroscopic degeneracy of the Heisenberg triangles proportionally reduces the bipartite quantum entanglement of their spins.

Following our recent paper [45], dealing with the spin-$1/2$ Ising–Heisenberg model on the planar lattice formed by trigonal bipyramids without a magnetic field, there is strong indication that the bipartite entanglement between the Heisenberg spins observed in the unique ferrimagnetic and macroscopically degenerated ferrimagnetic ground-state phases in the present paper will also persist at finite temperatures. Moreover, the interesting thermally induced reentrant behavior of the phenomenon can be expected near first- and second-order phase transitions of the system. Our future investigation will continue in this direction.