Thermodynamics of symmetric spin--orbital model: One- and two-dimensional cases

TL;DR

This paper investigates the thermodynamic properties of symmetric spin-orbital models in one- and two-dimensional systems, revealing features akin to phase transitions driven by quantum entanglement, even without long-range order.

Contribution

It provides a detailed analysis of thermodynamic behavior in spin-orbital models using a self-consistent approach, highlighting entanglement effects in phase transition-like features.

Findings

Thermodynamic features appear without long-range order.

Quantum entanglement influences thermodynamic behavior.

Behavior resembles phase transition signatures.

Abstract

The specific heat and susceptibilities for the two- and one-dimensional spin--orbital models are calculated in the framework of a spherically symmetric self-consistent approach at different temperatures and relations between the parameters of the system. It is shown that even in the absence of the long-range spin and orbital order, the system exhibits the features in the behavior of thermodynamic characteristics, which are typical of those manifesting themselves at phase transitions. Such features are attributed to the quantum entanglement of the coupled spin and orbital degrees of freedom.

Equations

$$\hat{\mathbf{H}}=J\sum _{<\mathbf{i},\mathbf{j}>}{\mathbf{S}}_{\mathbf{i}}{\mathbf{S}}_{\mathbf{j}}+I\sum _{<\mathbf{i},\mathbf{j}>}{\mathbf{T}}_{\mathbf{i}}{\mathbf{T}}_{\mathbf{j}}+K\sum _{<\mathbf{i},\mathbf{j}>}\left({\mathbf{S}}_{\mathbf{i}}{\mathbf{S}}_{\mathbf{j}}\right)\left({\mathbf{T}}_{\mathbf{i}}{\mathbf{T}}_{\mathbf{j}}\right),$$

$$\hat{\mathbf{H}}=\hat{J}+\hat{I}+\hat{K},$$

$$⟨{\hat{S}}_{\mathbf{i}}⟩=⟨{\hat{T}}_{\mathbf{i}}⟩=0$$

$$⟨{\hat{S}}_{\mathbf{i}}^{\alpha }{\hat{S}}_{\mathbf{j}}^{\beta }⟩=0,⟨{\hat{T}}_{\mathbf{i}}^{\alpha }{\hat{T}}_{\mathbf{j}}^{\beta }⟩=0,⟨{\hat{S}}_{\mathbf{i}}^{\alpha }{\hat{T}}_{\mathbf{j}}^{\beta }⟩=0.$$

$$G_{\mathbf{q}}={⟨S_{\mathbf{q}}^z\mid S_{-\mathbf{q}}^z⟩}_{\omega },$$

$$R_{\mathbf{q}}={⟨T_{\mathbf{q}}^z\mid S_{-\mathbf{q}}^z⟩}_{\omega }.$$

$$G_{\mathbf{q}}=\frac{F_{ac}(\mathbf{q})}{{\omega }^2-{\omega }_{ac}^2(\mathbf{q})}+\frac{F_{opt}(\mathbf{q})}{{\omega }^2-{\omega }_{opt}^2(\mathbf{q})},$$

$$R_{\mathbf{q}}=\frac{F_{ac}(\mathbf{q})}{{\omega }^2-{\omega }_{ac}^2(\mathbf{q})}-\frac{F_{opt}(\mathbf{q})}{{\omega }^2-{\omega }_{opt}^2(\mathbf{q})},$$

$$c_r=⟨{\hat{S}}_{\mathbf{i}}^z{\hat{S}}_{\mathbf{i}+\mathbf{r}}^z⟩,r=g,d,2g$$

$$m_0=⟨S_{\mathbf{i}}^z{T}_{\mathbf{i}}^z⟩,m_g=⟨S_{\mathbf{i}}^z{T}_{\mathbf{i}+\mathbf{g}}^z⟩,$$

$$G_{\mathbf{q}}={⟨S_{\mathbf{q}}^z\mid S_{-\mathbf{q}}^z⟩}_{\omega },$$

$$R_{\mathbf{q}}={⟨T_{\mathbf{q}}^z\mid S_{-\mathbf{q}}^z⟩}_{\omega },$$

$$G_{\mathbf{q}}=\frac{F_{ac}(\mathbf{q})}{{\omega }^2-{\omega }_{ac}^2(\mathbf{q})}+\frac{F_{opt}(\mathbf{q})}{{\omega }^2-{\omega }_{opt}^2(\mathbf{q})},$$

$$R_{\mathbf{q}}=\frac{F_{ac}(\mathbf{q})}{{\omega }^2-{\omega }_{ac}^2(\mathbf{q})}-\frac{F_{opt}(\mathbf{q})}{{\omega }^2-{\omega }_{opt}^2(\mathbf{q})},$$

$$F_{ac}=\frac{F_1+F_2}{2},F_{opt}=\frac{F_1-F_2}{2},$$

$$F_1=-8J{c}_g(1-{\gamma }_{\mathbf{q}})-M{m}_0,F_2=M{m}_0,$$

$${\omega }_{ac}^2(\mathbf{q})=W_1+W_2,{\omega }_{opt}^2(\mathbf{q})=W_1-W_2,$$

$$W_1$$

$$W_2$$

$$m_0=⟨S_{\mathbf{i}}^z{T}_{\mathbf{i}}^z⟩,m_g=⟨S_{\mathbf{i}}^z{T}_{\mathbf{i}+\mathbf{g}}^z⟩,$$

$${\omega }_{opt}(\mathbf{\Gamma })\ge {\omega }_{ac}(\mathbf{\Gamma })=0,{\omega }_{ac}(\mathbf{Q})\ge {\omega }_{opt}(\mathbf{Q})\ge 0.$$

$$F_{ac}=\frac{F_1+F_2}{2},F_{opt}=\frac{F_1-F_2}{2},$$

$$F_1=-4J{c}_g(1-{\gamma }_{\mathbf{q}})-M{m}_0,F_2=M{m}_0$$

$$W_1$$

$$W_2$$

Full text

I Introduction

The problem concerning the formation of entangled quantum states plays an important role in many fields of physics being especially attractive in the research related to quantum computations and their possible implementations Bengts06_Book ; Benent07_Book . The entanglement of states corresponding to different spin projections is most often considered. No less fascinating problems also arise in two-spin systems, where the entanglement is perceived as some correlated state of two spin variables You15_NJP .

The two-spin models themselves usually appear in the description of specific features of transition metal compounds with the coupled spin and orbital degrees of freedom; that is why such models are often referred to as spin–orbital ones Kugel82_UFN ; Oles09_APPA ; Oles12_JPCM . Unusual effects related to the spin–orbital correlations and the corresponding quantum entanglement are widely discussed in the current literature. In particular, the possibility of extraordinary spin–orbital quantum states and transitions between them was pointed out Belemu17_PRB ; Belemu18_NJP ; Brzezi18_JSNM .

It is interesting that an unusual behavior of spinorbital correlations can manifest itself also in systems without any spin or orbital order owing to their low dimensionality. For example, the symmetric two-dimensional model is characterized by vanishing spin–orbital correlations at certain threshold values of the temperature or parameters of the model Kaga14_JL_R . At the same time, the behavior of spin and orbital correlation functions themselves does not demonstrate any peculiarities.

his work is aimed at finding out the thermodynamic manifestations of such unusual specific features of spin–orbital correlations in the systems under study. Below, we demonstrate that, although long-range order is absent because of the low dimensionality of the system, a steep threshold type increase in spinorbital correlations at a certain intersubsystem exchange coupling constant is accompanied by clearly pronounced features in the thermodynamic characteristics resembling those typical of a phase transition. We relate these features to the formation of an entangled spinorbital state.

We study the one-dimensional (1D) and two-dimensional (2D) cases using the spherically symmetric self-consistent approach, which provides quite reliable results for low-dimensional systems  Kaga14_JL_R ; Bara11_TMP_R . Note that the specific features in the correlation functions and in the thermodynamic parameters manifest themselves in different ways depending on the dimensionality of the system under study. In particular, in the accepted approach, a reentrant (with respect to the temperature) transition to the state with nonzero intersubsystem correlations is possible in the 1D case in contrast to the 2D one.

II Model and method

We consider a symmetrical version of the quantum spin–orbital model (symmetrical Kugel–Khomskii model) on a square lattice (or a linear chain). The corresponding Hamiltonian has the form

[TABLE]

[TABLE]

where $<\mathbf{i},\mathbf{j}>$ denotes the summation over the nearest neighbor bonds on the square lattice or linear chain; $\widehat{\mathbf{S}}_{\mathbf{i}}$ and $\widehat{\mathbf{T}}_{\mathbf{i}}$ are respectively the spin and pseudospin operators for $S=1/2$ and $T=1/2$, respectively (the latter operators correspond to the orbital degrees of freedom).

We study the case of antiferromagnetic (AFM) spin–spin and pseudospin–pseudospin interactions, $J=I>0$, as well as of the negative exchange coupling between subsystems, $K<0$ (below, we always use the energy units corresponding to $J=I=1$). In Ref. Pati98_PRL, , it was demonstrated within the symmetric spin–orbital model that the interplay of spin and pseudospin (orbital) degrees of freedom is the most clearly pronounced at such relation between the parameters and results in the “entanglement” of these degrees of freedom You12_PRB ; Lundgr12_PRB .

According to the Mermin–Wagner theorem Mermin66_PRL , the long-range order is not possible for the decoupled subsystems ($K=0$) in the 2D and 1D cases at any nonzero temperature $T\neq 0$. It is natural to suppose that switching on the intersubsystem coupling, $K\neq 0$, leads to an additional contribution of thermal and quantum fluctuations, and the long-range order is still missing. Below, we limit ourselves only to the case of nonzero temperatures.

Thus, we adopt the following assumptions:

• •

(i) All lattice sites in the system under study are equivalent (the translational symmetry is not broken).

• •

(ii) All single-site averages are equal to zero (the long-range order is absent)

[TABLE]

• •

(iii) Correlation functions involving different spin and pseudospin components ($\alpha\neq\beta$) also vanish (SU(2) symmetries in the spin and pseudospin spaces are not broken)

[TABLE]

Conditions (i)–(iii) are automatically satisfied in the quantum spherically symmetric self-consistent approach (see, e.g., Refs. Bara11_TMP_R, ; Kondo72_PTP, ; Shimah91_JPSJ, ), which we use below.

Similarly to Ref. Kaga14_JL_R, , we consider the spin–spin and spin–pseudospin retarded Green’s functions

[TABLE]

[TABLE]

Obviously, $\left\langle T_{\mathbf{q}}^{z}\mid T_{\mathbf{-q}}^{z}\right\rangle_{\omega}=\left\langle S_{\mathbf{q}}^{z}\mid S_{\mathbf{-q}}^{z}\right\rangle_{\omega}$ since we consider the symmetric case $I=J$.

The calculations employing the conventional spherically symmetric self-consistent scheme (see, e.g., Refs. Bara11_TMP_R, ; Kaga14_JL_R, ) including the approximations for $\widehat{K}$ used in Ref. Kaga14_JL_R, and the dominant contribution of the on-site intersubsystem exchange coupling, result in the following expressions for $G_{\mathbf{q}}$ and $R_{\mathbf{q}}$:

[TABLE]

[TABLE]

Rather lengthy expressions for the numerators in the Green’s functions and for the spectra corresponding to the optical and acoustic branches of excitations are given in Appendix.

Expressions (6) and (7) involve spin–spin correlation functions. In the two-dimensional case,

[TABLE]

are the spin–spin correlation functions for the first (side of the square), $c_{g}\equiv c_{1}$, second (diagonal), $c_{d}\equiv c_{2}$, and third (doubled side), $c_{2g}\equiv c_{3}$, nearest neighbors, respectively. In the one-dimensional case, the diagonal correlation functions are, of course, absent.

The on-site and intersite spin–pseudospin correlation functions are

[TABLE]

respectively.

In the further numerical procedure, all correlation functions $c_{r}$, $m_{0}$, and $m_{g}$ are self-consistently determined in terms of the Green’s functions $G$ and $R$.

III 2D case. Results and discussion

An important qualitative result for the 2D case has been already reported in Ref. Kaga14_JL_R . Namely, at a fixed temperature and a small intersubsystem coupling parameter $|K|$, the spin–pseudospin correlation functions $m_{0}$ and $m_{1}$ vanish. However, when $|K|$ increases to a threshold value $K_{c}$, the correlation functions $m_{0}$ and $m_{1}$ become nonzero and increase according to a power law. A similar situation occurs at fixed $K$ with the decrease in the temperature. The correlation functions $m_{0}$ and $m_{1}$ start to increase in a power-law manner below the threshold temperature $T_{c}$.

In both cases, the intrasubsystem (spin–spin) correlation functions $c_{r}$ change only slightly and do not exhibit any anomalies at the critical point. For the sake of convenience, the corresponding plots taken from Ref. Kaga14_JL_R, are shown in Appendix.

Naturally, the arising intersubsystem correlations should appreciably affect the thermodynamic characteristics. In the framework of our approach, the energy is completely determined by one- and two-site correlation functions. Therefore, one could expect that the steep increase in the spin–pseudospin correlation functions in the entangled state leads to anomalies of the specific heat.

Figure 1 illustrates the evolution of the specific heat (a) as a function of $K$ at fixed and (b) as a function of $T$ at fixed $K$. As expected, at entering the range with intersubsystem correlations, the specific heat undergoes a jump. The higher the critical temperature $T_{c}$, the larger the jump. In Fig. 1b, we can also see that all curves have a common asymptotic behavior at high temperatures and $C\to 0$ at $T\to 0$ (the Nernst theorem is satisfied).

Figure 2 shows the behavior of spin–spin $\chi_{ss}$ and spin–pseudospin $\chi_{st}$ susceptibilities along the same paths in the phase diagram: (a) as functions of $K$ at fixed $T$ and (b) as functions of $T$ at fixed $K$. Similarly to the specific heat, the susceptibility exhibits a jump at critical points.

Outside the transition range, the susceptibility only slightly depends on $T$ and $K$. The behavior of the spin–spin susceptibility at small $|K|$ is in good agreement with that characteristic of the well-studied Heisenberg model Shimah91_JPSJ .

Thus, in the 2D case, i.e., for the square lattice, the state with entangled spin and pseudospin degrees of freedom arises at any temperature with the increase in the absolute value of the intersubsystem coupling constant $K$.

IV 1D case. Results and discussion

The calculations in the one-dimensional case are similar to those discussed above for the 2D case (see the corresponding expressions in Appendix).

Here, we adopt a mean-field approach implying the absence of decay of the spin excitations in the Green’s functions given by Eqs. (6) and (7). Nevertheless, due to underlying conditions (i)–(iii) (see Section II), such approach in the 1D case gives the results in good agreement with the exact solution or, in the situation with frustrations when the exact solution does not exist, with the numerical results for finite chains Suzuki94_JPSJ ; Junger04_PRB ; Hartel08_PRB .

Figure 3 shows phase diagrams for both the 1D and 2D cases. We can see that the one-dimensional case is qualitatively different from the two-dimensional one. Indeed, in contrast to the 2D case, the line separating the regions with zero and nonzero spin–pseudospin correlations in the 1D case begins at a nonzero intersubsystem exchange coupling constant (at $K\approx-2.44$ and $T=0$). Moreover, the phase boundary is characterized by a nonmonotonic behavior suggesting the possibility of a reentrant transition to the region with entanglement under variation of the temperature.

Of course, this nonmonotonicity manifests itself in other results. In Figs. 4a and 4b for the 1D chain (analogous to Figs. B1 and B3 for the 2D case presented in Appendix), we illustrate the evolution of spin–spin and spin–pseudospin correlation functions changing along the horizontal and vertical paths shown in Fig. 3(as a function of $K$ at fixed $T$ and as a function of $T$ at fixed $K$).

The plots in Fig. 4 are not as pictorial as similar ones for the 2D case shown in Fig. B1 because the curves corresponding to the temperatures at both sides from the inflection point in the phase boundary overlap. An additional qualitative difference from the 2D case is a kink in the spin–spin correlation function $c_{1}$ arising at the $K_{c}$ point simultaneously with the beginning of the power-law growth of on-site ($m_{0}$) and intersite ($m_{1}$) spin–pseudospin correlation functions.

The most clearly pronounced feature in Fig. 4b in the form of a “bubble” at $K=-2.3$ is related just to the reentrant transition (at entering the temperature range of entanglement and leaving it). In Fig. 4a, we can also see that the spin–spin correlation function $c_{1}$ exhibits a kink at the $K_{c}$ point and decreases steeply in absolute value with the increase in $|K|$ within the range of entanglement simultaneously with an increase in the spin–pseudospin correlation functions.

In the 2D case, only the intersubsystem correlations exhibit an anomalous behavior. Therefore, it is intuitively obvious that the simultaneous existence of anomalous features in both intra- and intersubsystem types of correlation functions in the 1D case can underlie substantial difference between these two cases, at least, in the behavior of the specific heat.

Indeed, this is seen in Fig. 5 (analogous to Fig 1 for the 2D case). Here, as well as in the two-dimensional case, we see a jump in the specific heat at entering the range of entanglement. However, in contrast to the two-dimensional case, the sign of such a jump depends on the transition temperature. This jump is positive and negative in the lower and upper part of the 1D phase boundary (see Fig. 3). Such a behavior is inconsistent with the usual concepts of the specific heat jump accompanying transitions to a more ordered state and can be attributed to the interplay between different-sign changes in intra- and intersubsystem correlations. The fast increase in the spin–pseudospin correlation functions dominates over the decrease in spin–spin correlation functions at low temperatures and vice versa at temperatures above that corresponding to the inflection point at the phase boundary.

In the 1D case, the asymptotic behaviors of the specific heat plots in the limits of high and low temperatures are the qualitatively the same as in 2D at all values of $K$: all curves at high temperatures tend to a common asymptotic curve and $C\to 0$ for all plots at $T\to 0$.

Finally, in Fig. 6 (analogous to Fig. 2 for the 2D case), we illustrate the behavior of spin–spin $\chi_{ss}$ and spin–pseudospin $\chi_{st}$ susceptibilities along the horizontal and vertical paths shown in Fig. 3. The plots in Fig. 6a, as those in the similar figure for the correlation functions, are not pictorial enough because of the overlap of curves corresponding to the temperatures lying on both sides of the inflection point in the phase boundary. Nevertheless, we can see the fast decrease in the spin–spin correlation function with increasing $|K|$.

The two-side steps in Fig 6b, as well as the “bubble” for the correlation function $m_{1}$ shown in Fig. 4b, are obviously due to the temperature-induced reentrant transition to the entanglement range and back.

In the rest aspects, the behavior of both susceptibilities is qualitatively the same as that in the two-dimensional case. Both the spin–spin $\chi_{ss}$ and spin–pseudospin $\chi_{st}$ susceptibilities exhibit a jump at the critical points.

V Conclusions

The results obtained in this work clearly demonstrate the possibility of a transition to the state with quantum entanglement in spin–orbital systems even in the absence of long-range order. Such a transition characterized by the threshold onset of spin–orbital correlations is clearly similar to usual phase transitions (which manifests themselves in the features of the specific heat and susceptibilities). In fact, it is a phase transition in some quantum fluid. Therefore, such phenomena can be expected in other systems with the quantum entanglement. At the chosen signs of the model parameters ($I,J>0$, $K<0$), the entangled state occurs in wide ranges of the temperature and model parameters. The possibility of quantum entanglement in spin–orbital models is widely discussed in the literature (see, e.g., Refs. Oles12_JPCM, ; Lundgr12_PRB, ; Pati98_PRL, ; Itoi00_PRB, ; Chen07_PRB, )) but usually at other relations between the parameters (in particular, when $K>0$ at the same signs of $I$ and $J$). Such choice gives rise to entanglement only within a rather limited region of the phase diagram. We emphasize that the studied features of the thermodynamic characteristics seem to be of a purely quantum nature. Indeed, the analysis of the exact solution of the one-dimensional spin–orbital model for the case of classical spins Kugel80_FNT does not reveal the same features in the behavior of the spin–orbital correlation function (its sharp vanishing at certain values of the temperature and model parameters) and the corresponding peculiarities in the thermodynamic characteristics.

Acknowledgments

This work was supported by the Russian Foundation for Basic Research, project nos. 17-52-53014 and 19-02-00509. K.I. Kugel also acknowledges the support of the Russian Foundation for Basic Research, project nos. 17-02-00323 and 17-02-00135.

Appendix A Expressions for the Green’s functions

A.1 General expressions

Spin-spin and spin–pseudospin retarded Green’s function are

[TABLE]

[TABLE]

Obviously, we have $\left\langle T_{\mathbf{q}}^{z}\mid T_{\mathbf{-q}}^{z}\right\rangle_{\omega}=\left\langle S_{\mathbf{q}}^{z}\mid S_{\mathbf{-q}}^{z}\right\rangle_{\omega}$, because we consider symmetric case $I=J$.

The self-consistent spherically symmetric approach leads to the following expressions for $G_{\mathbf{\mathbf{q}}}$ and $R_{\mathbf{q}}$:

[TABLE]

[TABLE]

A.2 2D case

The numerators for acoustic and optical branches are

[TABLE]

[TABLE]

and the excitations spectra

[TABLE]

[TABLE]

here $c_{r}=\langle\widehat{S}_{\mathbf{i}}^{z}\widehat{S}_{\mathbf{i+r}}^{z}\rangle,\quad r=g,d,2g$ are spin-spin correlation functions, respectively for first (side of the square) $c_{g}\equiv c_{1}$, second (diagonal) $c_{d}\equiv c_{2}$ and third (doubled side) $c_{2g}\equiv c_{3}$ nearest neighbors, $\widetilde{c}_{r}=\alpha_{r}{c}_{r}$ are correlation functions with vertex corrections. The lattice sum for square case $\gamma_{\mathbf{q}}=\frac{1}{4}\sum_{\mathbf{g}}e^{i\mathbf{qg}}=\frac{1}{2}(\cos(q_{x})+\cos(q_{y}))$.

Hereinbefore on-site $m_{0}$ and intersite $m_{g}\equiv m_{1}$ spin-pseudospin correlation functions are

[TABLE]

and for the intersubsystem vertex corrections $\widetilde{m}_{0}=\alpha_{ST}^{0}m_{0}$, $\widetilde{m}_{g}=\alpha_{ST}^{g}m_{g}$ we adopted the approximation $\alpha_{ST}^{0}=\alpha_{ST}^{g}=1$. For simplicity, we use the notation $M=8Km_{0}$.

Note, that the following relations for the symmetrical points $\mathbf{\Gamma}=(0,0)$ and $\mathbf{Q}=(\pi,\pi)$ in the Brillouin zone are always satisfied

[TABLE]

A.3 1D case.

[TABLE]

[TABLE]

[TABLE]

now $M=4Km_{0}$, $\gamma_{\mathbf{q}}$ is one-dimensional, other notations are the same as for 2D case.

Appendix B Three figures from Ref. 10

In Figs. B1–B3 almost invisible nonzero values of $m_{0}$ and $m_{1}$ just before the transition ($T\gtrsim T_{c}$ for Fig. B1 and $T\gtrsim T_{c}$ for Fig. B2) show the accuracy of the self-consistent calculations.