Low temperature pseudo-phase transition in an extended Hubbard diamond chain
Onofre Rojas, Jordana Torrico, L. M. Veríssimo, M. S. S. Pereira, S. M. de Souza, M. L. Lyra
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Low temperature pseudo-phase transition in an extended Hubbard diamond chain
Onofre Rojas and S. M. de Souza
Departamento de Física, Universidade Federal de Lavras, 37200-900, Lavras - MG, Brazil
Jordana Torrico,
Departamento de Física, Universidade Federal de Minas Gerais,C. P. 702, 30123-970, Belo Horizonte - MG, Brazil.
L. M. Veríssimo, M. S. S. Pereira, and M. L. Lyra
Instituto de Física, Universidade Federal de Alagoas, 57072-970 Maceió - AL, Brazil
We consider the extended Hubbard diamond chain with an arbitrary number of particles driven bychemical potential. The interaction between dimer diamond chain and nodal couplings is consideredin the atomic limit (no hopping), while the dimer interaction includes the hopping term. Wedemonstrate that this model exhibits a pseudo-transition effect in the low-temperature regime. Here,we explore the pseudo-transition rigorously by analyzing several physical quantities. The internalenergy and entropy depict sudden, although continuous, jumps which closely resembles discontinuousor first-order phase transition. At the same time, the correlation length and specific heat exhibitastonishing strong sharp peaks, quite similar to a second-order phase transition. We associatethe ascending and descending part of the peak with power-law "pseudo-critical"exponents. Wedetermine the pseudo-critical exponents in the temperature range where these peaks are developed,namely ν = 1 for the correlation length and α = 3 for the specific heat. We also study the behaviorof the electron density and isothermal compressibility around the pseudo-critical temperature. I. INTRODUCTION
In recent investigations of several decorated one-dimensional models with short-range interactions, thefirst derivative of free energy, such as entropy, inter-nal energy, and magnetization, shows a steep like func-tion of temperature but still with continuous changewhich is quite similar to a first-order phase transi-tion behavior. The second-order derivative of free en-ergy, like the specific heat and magnetic susceptibility,resembles a typical second-order phase transition be-havior at finite temperature. This peculiar behaviordrew attention to a more careful study, as consideredin reference [1]. In reference[2], an additional discus-sion of the above phenomenology focused in the behav-ior of correlation function for arbitrarily distant spinsaround the pseudo-transition. Similar pseudo-transitionswere shown to take place in the Ising-Heisenberg dia-mond chain[3, 4] and even in the pure Ising diamondchain[5]. It also has been explored in the one-dimensionaldouble-tetrahedral model, where the nodal sites are oc-cupied by localized Ising spins and alternate with apair of delocalized mobile electrons within a triangu-lar plaquette[6]. Similarly, ladder model with alternat-ing Ising-Heisenberg coupling[7], as well as the triangu-lar tube model with Ising-Heisenberg coupling [8] de-pict pseudo-transition signatures. A universal charac-ter of pre-asymptotic pseudo-critical exponents has beendemonstrated[9, 10]. These pseudo-transitions takingplace at finite-temperatures in one-dimensional modelsystems with short-range interactions are of a distinctnature from the true phase-transition exhibited in thepresence of long-range couplings for which the correlationlenght diverges, but e.g the specific heat can be without divergence [11]. In all the above model systems present-ing a pseudo-transition at finite temperature, one of thecouplings was assumed to be Ising-like in order to allowfor the exact calculation of the thermodynamic quanti-ties. Examples of pseudo-transitions taking place in one-dimensional systems of interacting electrons without theassumption of an Ising-like nature of relevant couplingsare still missing.The Hubbard model is one of the simplest models thatdescribe more accurately strongly interacting electronsystems, which have attracted a great deal of interest overthe past decades related to the possible emergence of ge-ometrical frustration properties [12, 13]. Magnetic frus-tration in highly correlated electron models arises due tothe geometric structure of the lattice, which induces sys-tem failures to satisfy simultaneously conflicting local re-quirements. The geometric frustration of Hubbard modelhas been extensively studied, particularly in the diamondchain structure, as considered by Derzhko et al.[14–16],where frustration for a particular class of lattice was dis-cussed. Montenegro-Filho and Coutinho-Filho[17] alsoconsidered the doped AB Hubbard chain, both in theweak coupling and the infinite- U limit (atomic limit)where quite interesting phases were identified as a func-tion of hole doping away from the half-filled band, as wellas 1/3-plateau magnetization, Kosterlitz-Thouless transi-tion, and Luttinger liquid[18]. Further, Gulacsi et al.[19]also discussed the diamond Hubbard chain in a magneticfield and a wide range of properties such as flat-bandferromagnetism, correlation-induced metallic and half-metallic processes. The thermodynamics of the Hubbardmodel on a diamond chain in the atomic limit was dis-cussed in reference [20]. Furthermore, frustrated quan-tum Heisenberg double-tetrahedral and octahedral chainsat high magnetic fields was discussed in[21]. Fermionicentanglement due to spin frustration was investigated inhybrid diamond chain with localized Ising spins and mo-bile electrons[22].On the other hand, generally rigorous analysis of theHubbard model is a challenging task. Only in a partic-ular case it is possible to obtain exact results [23]. Ear-lier in the seventies, Beni and Pincus[24] focused in theone-dimensional Hubbard model. Later Mancini [25, 26]discussed several additional properties of extended one-dimensional Hubbard model in the atomic limit, obtain-ing the chemical potentials plateaus of the particle den-sity, as a function of the on-site Coulomb potential at zerotemperature. Earlier, the spinless versions of the Hub-bard model on diamond chain also was investigated[27],as well as Lopes and Dias[28] performed a detailed inves-tigation using the exact diagonalization approach. TheIsing-Hubbard diamond chain has been investigated inreference [29]. In addition, experimental data regard-ing the / magnetization plateau in azurite[30, 31]were reproduced in several theoretical model systemssuch as the Ising-Heisenberg diamond chain[32–34]. Thequantum block-block entanglement was investigated inthe one-dimensional extended Hubbard model by ex-act diagonalization[35]. When the absolute value of thenearest-neighbor Coulomb interaction becomes small, theeffects of the hopping term and the on-site interactioncannot be neglected. The experimental observation of thedouble peaks both in the magnetic susceptibility and spe-cific heat [36–38] can be described accurately by the ex-tended Hubbard diamond chain model without the hop-ping of electrons or holes between the nodal sites.From an experimental point of view, the diamondchain structure is also motivating. Recently, thecompound Cu (CH COO) (OH) · O has been syn-thetized [39], which exhibits a unique one-dimensional di-amond chain structure. There are other compounds suchas Cu Cl (H O) · C SO trimer chain system[40],as well as the well known natural mineral azurite Cu (CO ) (OH) [37] which are well represented by 1Ddiamond chains. These and similar compounds wouldbe ideal physical systems on which pseudo-transitions atfinite temperature could be searched.Here, we advance in the study of pseudo-phase tran-sitions in quasi one-dimensional systems by presenting adetailed exact study of the thermodynamic properties ofthe diamond chain Hubbard model in the atomic limit.The present article is organized as follow: In sec.2, werevisit the extended Hubbard model on diamond chainstructure[20] with nodal sites considered in the atomiclimit. In sec.3 we present our main findings where we fo-cus in the existence of a pseudo-critical temperature byexploring the behavior of the correlation length. In sec.4(sec.5) we analyze first (second) derivative physical quan-tities in the vicinity of the pseudo-transition. Finally, sec.6 is devoted to our conclusions and perspectives. i + 1 i a b c U t V t UU V V Figure 1: Schematic representation of the extended Hubbardmodel on the diamond chain. Onsite Coulomb repulsion in-teraction is denoted by U and nearest neighbor repulsion in-teraction is represented by V and V , t stands for the electronhopping term. II. THE EXTENDED HUBBARD MODEL
In this section, we revisit the model considered in ref-erence [20]. Some results that will be used in the fol-lowing section are updated and summarized. The modelillustrated in fig.1, consider the hopping term t betweensites a and b . Additionally, there is an onsite Coulombrepulsion interaction U and nearest neighbor repulsioninteraction V between a and b , whereas V correspondsthe coupling between of nodal sites c with sites a and b (as labeled in the last block of fig.1). We also assumethat this model has an arbitrary particle density. Thus,the system will be described by including a chemical po-tential denoted by µ . The Hamiltonian of the proposedmodel can be expressed by: H = N X i =1 H i,i +1 , (1)with N being the number of unit cells (sites a , b and c ),and H i,i +1 is given by H i,i +1 = − t X σ = ↓ , ↑ (cid:16) a † i,σ b i,σ + b † i,σ a i,σ (cid:17) − µ (cid:0) n ai + n bi + n ci (cid:1) + U (cid:0) n ai, ↑ n ai, ↓ + n bi, ↑ n bi, ↓ + n ci, ↑ n ci, ↓ (cid:1) V n ai n bi + V ( n ai + n bi )( n ci + n ci +1 ) , (2)with a i,σ , and b i,σ ( a † i,σ and b † i,σ ) being the Fermi annihi-lation (creation) operators for electrons, while σ standsfor the electron spin, and n αi,σ stands for the numberoperator, with α = { a, b, c } . Using this number opera-tor, we also define conveniently the following operators n αi = n αi, ↑ + n αi, ↓ .In order to contract and symmetrize the Hamiltonian(2), we can define properly the following operators p i,i +1 = ( n ci + n ci +1 ) , q i,i +1 = ( n ci, ↑ n ci, ↓ + n ci +1 , ↑ n ci +1 , ↓ ) . (3)Using these operators, we can rewrite the Hamiltonian(2), which becomes as follows: H i,i +1 = − t X σ = ↓ , ↑ (cid:16) a † i,σ b i,σ + b † i,σ a i,σ (cid:17) − µ p i,i +1 − ( µ − V p i,i +1 ) (cid:0) n ai + n bi (cid:1) + V n ai n bi + U (cid:0) n ai, ↑ n ai, ↓ + n bi, ↑ n bi, ↓ (cid:1) + U q i,i +1 . (4)It is worth mentioning that this model already was in-vestigated in reference [20], for arbitrary number of elec-trons. Here we consider in each site the following basis { , ↑ , ↓ , ⇌ } .The eigenvalues and eigenvectors of the dimer plaque-tte are summarized in Table I, which are valid in generalfor arbitrary values of the Hamiltonian (1) parameters.It is worth to mention that all analysis performedthroughout this work will be done in the thermodynamiclimit. Finite size effects can be evaluated in a similar wayto that put forward in reference [41]. A. Phase Diagram
In order to analyze some relevant features, we will focusin the more interesting case when the particle-hole sym-metry is satisfied which follows the restriction V = V / .Under this condition, for instance, we can analyze thehalf-filled band case, which occurs under the followingrestriction for the chemical potential µ = U/ V .It is worthy to mention that, the Hamiltonian (2), has64 eigenvalues per diamond plaquette. The zero tem-perature phase diagram analysis was already discussedin reference [20]. Below we just give some ground-stateenergies relevant to the following analysis of the pseudotransition features.Fig.2a illustrates the zero temperature phase diagramin the plane t − µ , for fixed U = 1 , V = 0 . and <µ/U < . . There is a phase corresponding to a dimerantiferromagnetic ( AF M ) state, | AF M i = N Y i =1 | S ( − ) ll i| i . (5)In this case, there are two electrons with opposite spinsin the dimer sites, while nodal sites are empty. Thus, theelectron density per unit cell is ρ = 2 / . The correspond-ing eigenvalue is given by E AF M = − µ + V − t tan( θ ) . (6)There is also a phase corresponding to a dimer andnodal frustrated state ( F RU ), | F RU i = N Y i =1 1 √ (cid:16)(cid:12)(cid:12)(cid:12) σ i E − (cid:12)(cid:12)(cid:12) σ i E(cid:17) | τ i i . (7)This state also has two electrons: one electron is in adimer site while the other occupies a nodal site, both t/U µ / U t/U µ / U (a)(b) AFM FRU FRU FRU FRU FRU FRU AFM U = 1.0 ; V = 0.1
Figure 2: Zero temperature phase diagram in the plane t − µ ,assuming fixed U = 1 , V = 0 . . (a) Shows the region of µ/U values where ρ = 2 / phases appear; (b) The range of µ/U values where ρ = 4 / phases appear. with arbitrary spin orientation. The electron density ρ =2 / . The corresponding residual entropy in units of k B is S = ln(4) , whereas the ground-state energy becomes E F RU = − µ − t + V. (8)The vertical red line corresponds to t c = ( U − V ) / . . We will focus in this phase boundary between | AF M i and | F RU i . The corresponding interface resid-ual entropy is S = ln(4) , as we will verify ahead.Other surrounding states in the phase diagram are: | F RU i = N Y i =1 1 √ (cid:16)(cid:12)(cid:12)(cid:12) σ E − (cid:12)(cid:12)(cid:12) σ E(cid:17) | i (9)with n = 1 particle per unit cell and electron density ρ = 1 / . | F RU i = N Y i =1 1 √ (cid:16)(cid:12)(cid:12)(cid:12) ⇌ E − (cid:12)(cid:12)(cid:12) ⇌ E(cid:17) | τ i i (10)with n = 3 electrons per unit cell or ρ = 1 . The residualentropy in the frustrated phases (9) and (10) is given by S = ln(2) . m M ab Eigenvalues g Eigenvectors λ , = D , | S i = (cid:12)(cid:12)(cid:12)(cid:12) (cid:29) . λ ( ± )0 ,σ = D ± t | S ( ± )0 σ i = √ (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) σ (cid:29) ∓ (cid:12)(cid:12)(cid:12)(cid:12) σ (cid:29)(cid:19) λ σ,σ = D + V | S σσ i = (cid:12)(cid:12)(cid:12)(cid:12) σσ (cid:29) λ (1) , ⇌ = D + U | S (1)0 , ⇌ i = √ (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ⇌ (cid:29) − (cid:12)(cid:12)(cid:12)(cid:12) ⇌ (cid:29)(cid:19) λ (+) ll = D + V + 2 t cot( θ ) 1 | S (+) ll i = √ (cid:26) cos( θ ) (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ⇌ (cid:29) + (cid:12)(cid:12)(cid:12)(cid:12) ⇌ (cid:29)(cid:19) − sin( θ ) (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ↑↓ (cid:29) + (cid:12)(cid:12)(cid:12)(cid:12) ↓↑ (cid:29)(cid:19)(cid:27) λ ( − ) ll = D + V − t tan( θ ) 1 | S ( − ) ll i = √ (cid:26) sin( θ ) (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ⇌ (cid:29) + (cid:12)(cid:12)(cid:12)(cid:12) ⇌ (cid:29)(cid:19) + cos( θ ) (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ↑↓ (cid:29) + (cid:12)(cid:12)(cid:12)(cid:12) ↓↑ (cid:29)(cid:19)(cid:27) λ (2) ↓ , ↑ = D + V | S (2) ↓ , ↑ i = √ (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ↑↓ (cid:29) − (cid:12)(cid:12)(cid:12)(cid:12) ↓↑ (cid:29)(cid:19) . λ ( ± ) ⇌ ,σ = D + U + 2 V ± t | S ( ± ) ⇌ ,σ i = √ (cid:18)(cid:12)(cid:12)(cid:12)(cid:12) ⇌ σ (cid:29) ∓ (cid:12)(cid:12)(cid:12)(cid:12) σ ⇌ (cid:29)(cid:19) λ ⇌ , ⇌ = D + 2 U + 4 V | S ⇌ , ⇌ i = (cid:12)(cid:12)(cid:12)(cid:12) ⇌⇌ (cid:29) Table I: The first column means the number of particles in dimer plaquette, the second column describes the dimer plaquettemagnetization, third column reports the eigenvalues, while fourth column corresponds to the eigenvalues degeneracy. Thefifth column corresponds to the eigenvectors of dimer plaquette. Here D m = − ( µ − V p i,i +1 ) m − µ p i,i +1 + U q i,i +1 and cot (2 θ ) = U − V t . Similarly, in fig.2b we show the phase diagram for thesame set of fixed parameters but in the range . <µ/U < . where ρ = 4 / phases appear.In this region we observe also another dimer antiferro-magnetic ( AF M ) state given by | AF M i = N Y i =1 | S ( − ) ll i| ⇌ i , (11)where two electrons with opposite spins are located inthe dimer sites, and the other pair of electrons is locatedin the nodal site. The respective ground-state energybecomes E AF M = − µ + 5 V + U − t tan( θ ) . (12)There is also a dimer and nodal frustrated state( F RU ), | F RU i = N Y i =1 1 √ (cid:16)(cid:12)(cid:12)(cid:12) ⇌ σ i E − (cid:12)(cid:12)(cid:12) σ i ⇌ E(cid:17) | τ i i , (13)with corresponding ground-state energy E F RU = − µ − t + 5 V + U, (14)and electron density ρ = 4 / .The vertical red line corresponds to t = ( U − V ) / . . The phase boundary between | AF M i and | F RU i on which the residual entropy is S = ln(4) , as we willverify ahead.The additional phase states illustrated in this diagramis composed by n = 5 : | F RU i = N Y i =1 1 √ (cid:16)(cid:12)(cid:12)(cid:12) ⇌ σ i E − (cid:12)(cid:12)(cid:12) σ i ⇌ E(cid:17) | ⇌ i , (15)with ρ = 5 / . This frustrated phase has residual entropy S = ln(2) .Further details of the ground state phase diagrams canbe found in reference[20]. Pseudo-transitions are iden-tified in the close vicinity of the AF M − F RU and AF M − F RU phase-boundaries. B. Thermodynamics
In order to use the decoration transformation ap-proach, we write the Boltzmann factors for the extendedHubbard model on diamond chain as follows: w n ci , n ci +1 =e − β D + 4 (cid:16) e − β D + e − β ( D + U +2 V ) (cid:17) cosh( βt )+ e − β D (cid:0) e − βU + 3e − βV (cid:1) + e − β ( D +2 U +4 V ) + e − β ( D + V ) (cid:16) e − βϑ + t/ + e − βϑ − t/ (cid:17) , (16)with D m = − ( µ − V p i,i +1 ) m − µ p i,i +1 + U q i,i +1 . (17)Here, ϑ ± is defined by ϑ ± = U − V ± p ( U − V ) + 16 t t . (18)To solve the effective Hubbard model with up tofour-body coupling, we can use the transfer matrixmethod[42], similarly to that used in reference [24, 27].Therefore the symmetric Hamiltonian by exchanging i → i + 1 and i + 1 → i , leads to a symmetric transfer matrixwhich can be expressed by: W = w , w , w , w , w , w , w , w , w , w , w , w , w , w , w , w , , (19)where the elements of W is given by Eq. (16). Let usdefine a convenient notation, w , ( x ) = 1 + 2 x (cid:18) x z y (cid:19) (cid:18) γ + γ (cid:19) ++ x (cid:18) z + 1 y + 1 yzς + ςyz (cid:19) + x y z , (20)with x = e βµ , y = e βU , z = e βV , γ = e βt and ς = e β √ ( U − V ) +16 t . All other Boltzmann factors couldbe expressed in terms of w , ( x ) defined by Eq. (20) asfollows: w n ,n ( x ) = x ( n + n ) / y ⌊ n / ⌋ + ⌊ n / ⌋ w , (cid:16) xz n + n (cid:17) , (21)by ⌊ . . . ⌋ we mean the floor function for any real number.In order to carry out a reduced transfer matrix, weuse the symmetry of the system. The proposed Hamilto-nian is invariant with respect to electron spin orientationin nodal sites (site c). Therefore, the reduced transfermatrix becomes V = w , √ w , w , √ w , w , √ w , w , √ w , w , . (22)The determinant of the reduced transfer matrix be-comes a cubic equation of the form det( V − Λ) = (cid:0) Λ + a Λ + a Λ + a (cid:1) = 0 , (23)where the coefficients become a =2 w , w , + 2 w , w , + 2 w , w , − w , w , w , − w , w , w , ,a =2 w , w , + 2 w , w , + w , w , − w , − w , − w , ,a = − w , − w , − w , . (24) Consequently, the roots of the algebraic cubic equationmay be expressed as follow Λ j = 2 p Q cos (cid:16) φ − πj (cid:17) − a , j = 0 , , , (25)with φ = arccos (cid:18) R √ Q (cid:19) , (26) Q = a − a , (27) R = 9 a a − a − a . (28)It is verified that Q > in appendix A, which impliesthat all three roots must be different and real. We alsoanalyze which eigenvalues must be the largest and thelowest one. So, it is enough to restrict < φ < π inthe cubic solution without loosing its general solution, asdiscussed in appendix A. Other intervals just exchangethe cubic root solutions, as illustrated in table II. In theinterval < φ < π , the eigenvalues are ordered as Λ > Λ > Λ . III. PSEUDO-CRITICAL TEMPERATURE
In this section we will analyze the anomalous thermo-dynamic property of the present extended diamond chainHubbard model in the atomic limit[20]. In order to studythe thermodynamic properties, we will use the exact freeenergy f = − β ln(Λ ) as a starting point. Therefore, wewill proceed our discussion of thermodynamic propertiesas a function of temperature and chemical potential. Par-ticularly, we will analyze entropy, internal energy, corre-lation length, specific heat, electron density and isother-mal compressibility. We aim to study physical quantitiesaround the pseudo-critical temperature T p . We stressthat the present definition of "pseudo-critical" is differ-ent to that defined by Saito [43] to describe the critical-like behavior in approaching the spinodal point near thefirst order transition. Therefore, by using a perturbationapproach, we can find the eigenvalues of transfer matrix(22), as discussed in Appendix B.In general, pseudo-transitions can be manipulated andanalyzed using perturbation techniques, as detailed inappendix B. For this purpose, we consider as unper-turbed matrix V defined in appendix eq.(B1). Theeigenvalues of the matrix B1, are given by (B2-B4), wherewe can observe the largest and second-largest eigenval-ues are given by (B2) and (B3), respectively. How-ever, in low temperature region we have the condition w , ≪ { w , , w , , w , } . This is so because the lead-ing term of w , comes from the lowest energy contribu-tion of w , that includes only excited states, while theleading term of w , , w , and w , is given by the cor-responding ground state energies of the system. There-fore, w , becomes exponentially small as compared tothe other elements in the low-temperature regime. When(B2) and (B3) become eventually equal, we can get apseudo-critical temperature from (B5), which reads u (0)0 = u (0)1 . (29)Similar as discussed in [1], the relation (29) can induceus to believe that there is a true phase transition at finitetemperature. However, the condition (29) does not meanthat the transfer matrix eigenvalues satisfy Λ = Λ : thecondition (29) is satisfied only when the matrix (B8) is ig-nored and for w , = 0 . Therefore, the first order pertur-bative corrections (see eqs.(B9) and (B11)) are u (1)0 > and u (1)1 < , which implies that the approximate solu-tion given in (B12) must satisfy Λ > Λ .Returning to the relation (29), we get the followingequation w , = 2 w , for w , > w , , (30)and this one corresponds to the vicinity of the boundaryline between AF M and F RU . The pseudo-critical tem-perature T p is illustrated in fig.3a where we plot T p × − (red line) as a function of chemical potential µ , for fixedparameters t = 0 . , U = 1 and V = 0 . . In fig.3b wereport T p for other values of t .Similarly, the relation (29) leads to the following equa-tion w , = 2 w , , for w , < w , , (31)which holds in the close vicinity of the F RU − AF M boundary-line.In fig.3a, we illustrate the pseudo-critical temperature T p × − as a function of chemical potential µ . The con-tinuous line corresponds to the pseudo-transition that oc-curs at low chemical potentials for which there are nearly2 electrons per unit cell (near the boundary between AF M and F RU phases). The dashed line accounts forthe second pseudo-transition appearing at larger chemi-cal potential for which there are 4 electrons per unit cellin the ground-state (near the boundary between AF M and F RU phases). In this latter case, the chemical po-tential was conveniently shifted to µ − V − U = µ − . (top scale) to allow representing both pseudo-transitionsin the same frame.In panel (b) the pseudo-critical temperatures T p × − as a function of t , for µ = 0 . (red line) and µ = 1 . (blue line) are illustrated.The pseudo-critical tempera-ture vanishes linearly T p ∝ ( t − t c ) as one approaches theboundary-line t c = ( U − V ) / from above, i.e., withinone of the AF M ground-states.
A. Correlation length
Since the eigenvalues are non-degenerate, the correla-tion length ξ can be obtained by using the largest Λ µ T P x10 - t T P x10 - µ = 0.18 µ = 1.28 (a)(b) U = 1.0 ; V = 0.1 µ - 1.1 Figure 3: (a) Pseudo-critical temperature T p × − for fixedparameters t = 0 . , U = 1 and V = 0 . near AF M and F RU phases boundary (solid line) and near the AF M and F RU phases boundary (dashed line). In the latter the chem-ical potential was shifted to µ − . (top scale). (b) Pseudo-critical temperature T p × − as a function of t for µ = 0 . (solid line) and µ = 1 . (dashed line). and second largest Λ eigenvalues given from (25). It issimply written as ξ = (cid:20) ln (cid:18) Λ Λ (cid:19)(cid:21) − . (32)In fig. 4 we illustrate the correlation length as a func-tion of temperature. Panel (a) reports data for severalvalues of t and assuming fixed µ = 0 . , U = 1 and V = 0 . . Here we observe how the peak becomes morepronounced when t → t c = ( U − V ) / . (approach-ing the F RU − AF M ground-state phase-boundary).For larger t the height of the peak becomes lower andbroader. For t = 0 . we already observe a strong peakaround T p , which was computed with high precision tobe T p = 1 . × − . Similarly, panel (b)reports for the same parameters set used in (a) but for µ = 1 . (in the vicinity of the F RU − AF M phase-boundary).There are two situations where pseudo-transition oc-curs. The first one satisfy the condition w , ∼ w , butin the perturbation regime for which w , ≪ ( w , − w , ) . Under this condition, we have the following re- -4 -3 -2 -1 T -1 ξ ( T ) t = 0.303t = 0.305t = 0.307t = 0.310t = 0.320 -4 -3 -2 -1 T -1 ξ ( T ) (a)(b) µ = 0.18 µ = 1.28 Figure 4: Correlation length as a function of temperaturein logarithmic scale for fixed U = 1 and V = 0 . . (a) Forseveral values of t , and fixed µ = 0 . (in the vicinity of the F RU − AF M phase-boundary). (b) For several values of t , and fixed µ = 1 . (in the vicinity of the F RU − AF M phase-boundary). sult, Λ Λ → ( w , w , , w , > w , w , w , , w , < w , . (33)Consequently, in the regime where the perturbationapproach stands, the correlation length in the close vicin-ity of the pseudo-critical temperature can be expressedby ξ ( τ ) = c ,ξ | τ | − + O ( τ ) . (34)Here τ = ( T − T p ) /T p and the coefficient is given by c ,ξ = ˜ w , T p ˜ v , (35)with ˜ w , is w , evaluated at T = T p , and ˜ v = (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( w , − w , ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) T p . (36)Similarly, the second pseudo-transition occurs when w , ∼ w , but in the perturbation regime w , ≪ -4 -3 -2 -1 | τ | -2 -1 ξ ( T ) τ > 0 (T > T P ) τ < 0 (T < T P ) ξ α |τ| -1 -4 -3 -2 -1 | τ | -2 -1 ξ ( T ) (a)(b) µ = 0.18 µ = 1.28 Figure 5: Correlation length as a function of τ , for fixed pa-rameters U = 1 , V = 0 . , t = 0 . . Red line correspondsfor τ > ( T > T p ) and dashed blue line denotes τ < ( T < T p ), while straight dash-doted line reports ξ ∝ | τ | − .(a) For µ = 0 . . (b) For µ = 1 . . ( w , − w , ) , that provides us Λ Λ → ( w , w , , w , > w , w , w , , w , < w , . (37)Analogously, the correlation length close to the pseudo-transition, can be expressed around T p , resulting in ξ ( τ ) = c ,ξ | τ | − + O ( τ ) , (38)where the coefficient is given by c ,ξ = ˜ w , T p ˜ v . (39)with ˜ w , evaluated at T = T p , and ˜ v = (cid:12)(cid:12)(cid:12)(cid:12) ∂ ( w , − w , ) ∂T (cid:12)(cid:12)(cid:12)(cid:12) T p . (40)These pre-asymptotic power-laws fail for very small τ ,because the perturbation conditions can not be satisfiedat T p . Therefore, the correlation length actually depictsa pronounced peak and not a true divergence.In fig.5a the correlation length is depicted as a func-tion of τ in logarithmic scale, assuming fixed parameters U = 1 , V = 0 . , t = 0 . and µ = 0 . . The solidcurve represents τ > ( T > T p ) while the dashed curvedenotes τ < ( T < T p ). The straight dash-doted linereports the asymptotic limit ξ ∝ | τ | − , where we observeclearly the critical exponent ν = 1 in the pseudo-criticalregime. Notice that when T is very close to T p the powerlaw critical exponent fails, evidencing the ultimate non-singular behavior of the thermodynamic quantities. Asimilar plot is depicted in fig. 5b, assuming the same setof parameters but for chemical potential µ = 1 . . Hereagain we observe manifestly the critical exponents ν = 1 .Note that the power-law behavior holds for nearly twodecades.In order to have a clear picture of the range of tem-peratures around T p on which the pseudo-critical regimeholds, we can estimate crossover boundaries between thenon-singular and power-law regimes, as well as betweenthe power-law regimes and the low and high temperatureones. The non-singular regime in the close vicinity of T p emerges as the perturbation analysis fails. In the vicinityof the F RU − AF M ground-state transition, crossoverlines can be built using w , = ( w , − w , ) , whichseparates the perturbative from the non-perturbativeregimes. A similar crossover line can be written in thevicinity of the F RU − AF M ground-state transition[ w , = ( w , − w , ) ]. The power-law regime holdswhile the correlation length remains much larger the lat-tice spacement. We will consider that crossover lines fromthe power-law to the low and high-temperature regimessatisfies ξ = 10 . Eqs. (32) and (33) can be used to drawthe respective crossover lines.In fig.6 we show the above crossover lines as a func-tion the hopping parameter t in the close vicinity’s of the F RU − AF M and F RU − AF M boundary lines foran illustrative set of the other Hamiltonian parameters.Firstly notice that the width of the non-singular regimedecreases quite fast as one approaches t c = ( U − V ) / .While the pseudo-transition temperature T p decreaseslinearly, the non-singular temperature range decays muchfaster as e − a/ ( t − t c ) , with a being a constant at t c . Thenon-singular regime widens as we further depart from t c . On the other hand, the crossover lines to the highand low-temperature regimes are weakly dependent on t .Therefore, in the close vicinity of t c , the power-law regimeextends over a temperature range of the order of a fewpercents of T p . The crossover lines delimiting the non-singular and the high and low temperature regimes even-tually meet. After this point, no pre-asymptotic power-law regime can be identified. IV. FIRST-ORDER LIKEPSEUDO-TRANSITION
It is interesting to analyze some physical quantities,which are obtained as a first derivative of the free energy t -0.10-0.050.000.050.10 ( T - T P ) / T P t -0.15-0.10-0.050.000.050.100.15 ( T - T P ) / T P pseudo-criticalpseudo-criticalhigh-temperaturelow-temperature non-singular µ =1.28 (b)(a) µ =0.18non-singularpseudo-criticallow-temperaturehigh-temperaturepseudo-critical Figure 6: Distinct temperature regimes as a function of thehopping parameter t in the close vicinity of the F RU − AF M (a) and F RU − AF M (b) ground-state transitions.Here we used U = 1 , V = 0 . , (a) µ = 0 . and (b) µ = 1 . .The non-singular regime becomes exponentially narrow as t → t c = ( U − V ) / . . The power-law regime holdsfor τ = ( T − T p ) /T p of the order of a few percents in the closevicinity of t c . It dies away far from t c when the perturbationand large correlation length conditions become incompatible. and exhibit an almost step-like behavior. A. Entropy
In fig.7a we illustrate the density plot of entropy( S = − ∂f∂T ) in the plane t − µ , for fixed parameters V = 0 . , U = 1 and T = 0 . , for µ = [0 , . . Thisplot is depicted in the same scale of the phase diagramillustrated in fig.2a. We definitely observe a pseudo-transition in the boundary of quasi-antiferromagnetic qAF M and quasi-frustrated qF RU phases character-ized by a sharp boundary, while electron density remainsconstant ρ = 2 / . The regions are mostly governed bythe zero temperature phases but, due to thermal fluc-tuations, the zero temperature phases become "quasi"long-ranged because of the lack of actual spontaneouslong-range order at any finite temperature. A similar µ S µ t qAF M qF RU qF RU qF RU .
26 0 . . . . . . . t . . . . . qAF M qF RU qF RU qF RU . . . . . . . . . (a) (b) Figure 7: Density plot of entropy in the plane t − µ , assumingfixed V = 0 . , U = 1 and T = 0 . . (a) For µ = [0 , . . (b)For µ = [1 . , . . density plot is illustrated in fig.7b, assuming the sameparameters but in the interval of µ = [1 . , . . We canobserve a steepy entropy change in the boundary betweenthe qAF M and qF RU phases where the electron den-sity ρ = 4 / per cell remains unaltered.In fig.8a we report the magnitude of entropy as a func-tion of temperature in semi-logarithmic scale, assumingfixed parameters t = 0 . , µ = 0 . , U = 0 . . Clearlyone can observe a nearly step behavior at T p = 1 . × − .For temperatures below T p the entropy is almost null,which corresponds to the quasi-antiferromagnetic phase( qF AM ). On the other hand, the entropy has a plateauregion ( qF RU ) with S = ln(4) for T > T p . As the tem-perature is further increased, the entropy behaves likein standard models. Therefore we observe a strong con-tinuous change of entropy around the pseudo-transition,which typically resembles a discontinuous (first-order)phase transition. However, we stress that it remains an-alytical at T p . -4 -3 -2 -1 T S ( T ) -4 -3 -2 -1 T -0.565-0.560-0.555-0.550 u ( T ) (a)(b) U = 1.0 ; V = 0.1t = 0.303 ; µ = 0.18 Figure 8: (a) Entropy as a function of temperature in semi-logarithmic scale, for fixed t = 0 . , µ = 0 . , U = 1 and V = 0 . . (b) Internal energy as a function of temperature insemi-logarithmic scale, for fixed t = 0 . , µ = 0 . , U = 1 and V = 0 . . B. Internal energy
Other quantity we discuss here is the internal energy u ( T ) = T S + f , which can also be obtained after a first-derivative of the free-energy.In fig.8b, the internal energy as a function of tem-perature in semi-logarithmic scale is depicted, for fixed t = 0 . , µ = 0 . , U = 1 and V = 0 . . Internal energyis continuous as a function of temperature, although italso shows a steep variation around the pseudo-criticaltemperature T p . For temperatures below T p , the inter-nal energy leads to u = − . , which corresponds tothe AF M ground state energy E AF M . For T > T p ,the internal energy is mostly dominated by the F RU ground-state energy E F RU = − . . C. The electron density
Another important quantity we explore is the thermalaverage electron density ρ = − (cid:16) ∂f∂µ (cid:17) per unit cell.In fig.9 the electron density as function of temperatureis reported in semi-logarithmic scale, for fixed t = 0 . , U = 1 , V = 0 . and µ = 0 . . At low temperatures,0 -4 -3 -2 -1 T ρ ( T ) -4 -3 -2 -1 T ρ ( T ) T T (a)(c) µ = 1.28 µ = 0.12 (b)(d) Figure 9: Electron density as a function of temperature insemi-logarithmic scale, for fixed t = 0 . , U = 1 and V =0 . . (a) For µ = 0 . . (b) Zooming panel (a) around T p . (c)For µ = 1 . . (d) Zooming panel (c) around T p . the electron density remains almost constant ρ = 2 / ,up to roughly around T ∼ . , and then the electrondensity decreases to ρ ≈ . , for higher temperaturethe electron density increases with temperature ( for fur-ther details see reference [20]). Therefore apparently nopseudo-critical behavior at T p is evidenced in ρ . In fig.9bwe plot the electron density as a function of temperaturearound pseudo-critical temperature. We observe a tinydepression at T p reminiscent of the pseudo-transition dueto thermal excitations to states with a single electron perunit cell (see phase-diagram).Similarly, the electron density is depicted in fig.9c asfunction of temperature in semi-logarithmic scale, assum-ing same set of parameters, but for µ = 1 . . Once againelectron density remains almost constant ρ = 4 / , fortemperatures below T ∼ . , while for higher temper-ature there is a peak reaching ρ ≈ . . In fig.9d wepresent a zooming plot of panel (c) around the pseudo-critical temperature, where a tiny peak reminiscent ofthe pseudo-critical temperature is also present refletingthermal excitations to states with 5 electrons per unitcell (see phase diagram). Surely we cannot expect anystrong change in the electron density around T p becausethe two competing grounds state have equivalent electrondensities. V. SECOND-ORDER LIKEPSEUDO-TRANSITION
Now let us turn our attention to the second-orderderivative of the free energy. The following quantitiesexhibit trends quite similar to second-order phase tran-sition. It is important to reinforce that there is no sin-gularity at T p but just a rather sharp peak. -3 -2 -1 T -6 -4 -2 C ( T ) t = 0.303t = 0.305t = 0.307t = 0.310t = 0.320 -3 -2 -1 T -6 -4 -2 C ( T ) -1 t = 0.303 (a)(b) µ = 0.18 µ = 1.28 t = 0.303 Figure 10: Specific heat as a function of temperature in log-arithmic scale for U = 1 and V = 0 . . (a) For several valuesof t and µ = 0 . (in the vicinity of the F RU − AF M phase-boundary). Inset: zooming around T p . (b) For severalvalues of t and µ = 1 . (in the vicinity of the F RU − AF M phase-boundary). A. The specific heat
In fig. 10, we display the specific heat C = T (cid:16) ∂ S ∂T (cid:17) asa function of temperature in logarithmic scale assumingfixed U = 1 , V = 0 . . In panel (a) we depict for severalvalues of t and µ = 0 . (close to the F RU − AF M ground-state transition). Here we see how the height ofthe peak increases and the peak becomes sharper when t → . , while for t larger the height of the peak be-comes lower and broader. At low temperatures, we ob-serve clearly a huge sharp peak quite similar to a second-order phase transition divergence around T p . However,a zooming look as provided in the inner plot around T p , evidences the rounded nature of the peak at the pseudo-critical temperature. Panel (b) reports the specific heatfor the same set of parameters used in panel (a) but for µ = 1 . (close to the F RU − AF M ground-state phasetransition).Now let us analyze the nature of the peak around T p ,looking for some critical exponent universality. For thispurpose, we consider the specific heat around pseudo-critical temperature in asymptotic limit (but not veryclose to T p ), which can be expressed according the dis-cussion in reference [9] as1 -4 -3 -2 -1 | τ | -8 -6 -4 -2 C ( τ ) τ > 0 (T > T P ) τ < 0 (T < T P )C α | τ | -3 -4 -3 -2 -1 | τ | -8 -6 -4 -2 C ( τ ) (a)(b) µ = 0.18 µ = 1.28 Figure 11: Logarithmic specific heat as a function of ln( | τ | ) ,for fixed parameters U = 1 , V = 0 . , t = 0 . . Solid linecorresponds for τ > ( T > T p ) and dashed line denotes τ < ( T < T p ). The straight dash-doted line reports ξ ∝ | τ | − . (a)For µ = 0 . . (b) For µ = 1 . . C ( τ ) = T (cid:18) ∂ S ∂τ (cid:19) (cid:18) ∂τ∂T (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) T p = 2 c f τ − , (41)where the coefficient c f is given by c f = ( ˜ w , c ,ξ , w , ∼ w , ˜ w , c ,ξ , w , ∼ w , . (42)Therefore the specific heat has the following pre-asymptotic expression when approaching the pseudo-critical temperature C ( τ ) ∝ | τ | − α , (43)with critical exponent α = 3 . The Hubbard diamondchain in the atomic limit also satisfy the pseudo-criticalexponent found in reference [9]. However, it is worthstressing that this pre-asymptotic regime is only validaround the ascending and descending parts of the peak,and surely fails very close to the peak top where theperturbation condition can not prevail.In order to confirm the above result. We report infig.11a the C ( τ ) as a function of τ , assuming fixed pa-rameters U = 1 , V = 0 . , t = 0 . and µ = 0 . (near the AF M − F RU phase-boundary). Continuous linerepresents data above T p , while dashed line correspondsdata below T p . The dash-doted line describes the asymp-totic function, with critical exponent α = 3 . The straightline with angular coefficient α = 3 fits accurately dataover nearly two decades. Very close to T p the power-law behavior breaks down, reflecting the actual analyticbehavior of the specific heat. In panel (b) we plot the spe-cific heat for µ = 1 . (near the AF M − F RU phaseboundary), showing similar trends. B. Compressibility
Another interesting quantity to be discuss is theisothermal electron compressibility[42] κ T = 1 ρ (cid:16) ∂ρ∂µ (cid:17) T ,as a function of Hamiltonian parameters, temperatureand electron density.In fig.12a the isothermal compressibility is shown asa function of temperature in semi-logarithmic scale, weobserve a tiny peak at the pseudo-critical temperature.Magnifying around the pseudo-critical temperature (seefig. 12b) it illustrated a double peak with local minimumat T p . Since the pseudo-transition results from competingground-states with the same electron density, we actuallycannot expect any giant peak of the electron compress-ibility at T p because that would result in a pronouncedelectron density change. VI. SUMMARY AND CONCLUSIONS
In summary, we considered the extended Hubbard di-amond chain restricted to the atomic limit with an ar-bitrary number of particles driven by chemical potential.The interaction between dimer diamond chain and nodalcouplings are taken in the atomic limit (no hopping),while dimer interaction includes the hopping term. Weshowed that this model exhibits a pseudo-transition ef-fect in the low-temperature region. The internal energyand entropy were shown to change quite abruptly in avery narrow range of temperatures around the pseudo-transition when the physical parameters are properlytuned in a close vicinity of special ground-state phasetransitions. The correlation length and specific heat dis-play pronounced peaks with well defined power-laws ina well defined temperature range in the vicinity of thepseudo-transition point. The pseudo-critical exponentsassociated with the correlation length and specific heatwere shown to be ν = 1 and α = 3 , respectively. Theseare the same pseudo-exponents reported to hold is spinmodels with Ising-like interactions, pointing towards auniversal behavior of pseudo-transitions[9, 10]. We alsodemonstrated that the electron density and respectiveelectronic compressibility displays reminiscent signaturesof the pseudo-transition. The present results add to thegeneral understanding of the remarkable phenomenon of2 -4 -3 -2 -1 T κ T ( T ) T κ T ( T ) (a)(b) t=0.303 ; µ = 0.18 Figure 12: Isothermal compressibility as a function of temper-ature in semi-logarithmic scale, for fixed t = 0 . , µ = 0 . , U = 1 and V = 0 . . (b) Magnified around T p . pseudo-transitions taking place at finite temperatures inone-dimensional equilibrium systems having structuredinteractions within the relevant unit cell. VII. ACKNOWLEDGEMENTS
This work was supported by CAPES (Coordenação deAperfeiçoamento de Pessoal de Nível Superior), CNPq(Conselho Nacional de Desenvolvimento Científico e Tec-nológico), FAPEAL (Fundação de Apoio à Pesquisa doEstado de Alagoas), and FAPEMIG (Fundação de Apoioà Pesquisa do Estado de Minas Gerais).
Appendix A: Real roots of cubic equation
Transfer matrix is a symmetric matrix, so its eigenval-ues are guaranteed to be real values. Therefore, here weverify the cubic equation (23) roots must be different andreal numbers, given by (25).Obviously, we can convince that both Q and R are realnumbers from (27) and (28).Now let verify the real number Q is positively definedby using the cubic equation coefficients (24), thus aftersome algebraic manipulation we obtain the following ex- φ Largest 2nd Largest Lowest h , π i Λ Λ Λ h π, π i Λ Λ Λ h π, π i Λ Λ Λ h π, π i Λ Λ Λ h π, π i Λ Λ Λ h π, π i Λ Λ Λ Table II: Cubic root solutions are tabulated in decreasing or-der by intervals in φ . pression Q = (cid:0) w , − w , − w , (cid:1) + ( w , − w , ) + w , + w , + w , , (A1)which is strictly a positive number ( Q > ), since allBoltzmann factors are positive. According to a cubicequation property, this condition is enough to concludethat three the eigenvalues of transfer matrix will be dif-ferent and real numbers.Now the question is to determine which are the largestand the lowest eigenvalues? Since all eigenvalues mustbe real and different because Q > . Which implies thatthe solution (25), must satisfy the following restriction φ = ± nπ , with n = { , , , . . . } . Therefore, we canchoose φ conveniently such that < φ < π . Thus wehave the following trigonometric relation relation cos (cid:16) φ (cid:17) > cos (cid:16) φ − π (cid:17) > cos (cid:16) φ − π (cid:17) , (A2)if we multiply all terms of inequalities by √ Q > , andadding − a to all expressions (note that a < ), thuswe conclude that Λ > Λ > Λ . (A3)Note that if φ = 0 , in principle we would have Λ = Λ but this condition is forbidden because Q > . Similarlyfor φ = π we have Λ = Λ , but again this conditioncannot be satisfied since Q > . Therefore, we haveidentified which eigenvalues is the largest one and thelowest one.Of course we can choose other intervals equivalently,and verify how the eigenvalues are ordered. In table IIis reported the eigenvalues for each intervals. The openintervals of φ guaranties all eigenvalues must be real anddifferent values, in order to satisfy the condition Q > .The three cubic root solutions are exchanging periodi-cally which depends on interval of φ , making a bit puzzleto identify which eigenvalues is the largest one. In tableII it is reported how the eigenvalues are ordered. In eachinterval the solutions are equivalent, because other inter-vals simply exchange the eigenvalues with no relevance.Therefore, here we conveniently choose the first interval φ ∈ h , π i , so the eigenvalues must be ordered as follow Λ > Λ > Λ .3Even more, according to Perron-Frobenius theorem,the largest eigenvalue must be non-degenerate and posi-tive. Appendix B: Perturbative Transfer MatrixCorrection
In order to study the pseudo-transitions property inlow-temperature region, we need to analyze the transfermatrix in the low temperature region. In principle, cubicroot solution could be a bit cumbersome task to identifyas an easy handling solutions. A simple strategy to ob-tain a reasonable solution could be considering the trans-fer matrix (19), as a sum o two matrix V = V + ζ V .The first term is given by V = w , w , w , w , w , , (B1)This matrix can be considered as the unperturbed trans-fer matrix, whose eigenvalues are u (0)0 = [ w , + w , + s ] , (B2) u (0)1 =2 w , , (B3) u (0)2 = [ w , + w , − s ] , (B4)where s = q d + 4 w , and d = w , − w , .In the limit of w , → , the unperturbed solution caneventually satisfy the following relation u (0)0 = u (0)1 . (B5)This condition leads to a pseudo-transition, which we cansimplify: w , = 2 w , , for w , > w , , (B6) w , = 2 w , , for w , < w , . (B7) The above condition is essential to find pseudo-criticaltemperature T p .The second term of transfer matrix is given by V = √ w , √ w , √ w , √ w , . (B8)Close to the pseudo-critical temperature, the elementsof matrix V become smaller than the elements of matrix V . Therefore, we can find the transfer matrix eigenval-ues by using a perturbative approach.Consequently, the solution of the transfer matrix aftera standard perturbative manipulation up to first orderterm gives us the following root corrections, u (1)0 = [( s − d ) w , + 2 w , w , ] s ( s − d ) (cid:16) u (0)0 − w , (cid:17) , (B9) u (1)1 = − u (1)0 − u (1)2 , (B10) u (1)2 = [( s + d ) w , − w , w , ] s ( s + d ) (cid:16) u (0)2 − w , (cid:17) . (B11)As a consequence of the perturbative correction, thetransfer matrix eigenvalues becomes Λ j = u (0) j + ζ u (1) j + O ( ζ ) , j = { , , } . (B12)Assuming the elements of matrix V are small, we canobtain an approximate result of (25), fixing ζ = 1 . Hence,it is evident that the largest eigenvalue of the transfermatrix will be Λ , even when u (0)0 = u (0)1 , since u (1)0 > and u (1)1 < . [1] S.M. de Souza and O. Rojas, Sol. State. Comm. , 131(2018).[2] I.M. Carvalho, J. Torrico, S.M. de Souza, O. Rojas, andO. Derzhko, Ann. Phys. , 45 (2019).[3] J. Torrico, M. Rojas, S.M. de Souza, O. Rojas, and N.S.Ananikian, Eur. Phys. Lett. , 50007 (2014).[4] J. Torrico, M. Rojas, S.M. de Souza, and O. Rojas, Phys.Lett. A , 3655 (2016).[5] J. Strečka, Acta Phys. Pol. A , 610 (2020).[6] L. Galisova and J. Strecka, Phys. Rev. E , 022134(2015).[7] O. Rojas, J. Strečka and S.M. de Souza, Sol. State.Comm. , 68 (2016).[8] J. Strecka, R.C. Alecio, M.L. Lyra, and O. Rojas, J. Magn. Magn. Mater. , 124 (2016).[9] O. Rojas, J. Strecka, M.L. Lyra, and S.M. de Souza,Phys. Rev. E , 042117 (2019).[10] T. Hutak, T. Krokhmalskii, O. Rojas, S.M. de Souza,and O. Derzhko, Phys. Lett. A 387, 127020 (2021).[11] I. Daruka and Z. Gulacsi, Phys. Rev. E , 5403 (1998).[12] I.S. Hagemann, Q. Huang, X.P.A. Gao, A.P. Ramirez,and R.J. Cava, Phys. Rev. Lett. , 894 (2001).[13] J. van Lierop and D.H. Ryan, Phys. Rev. Lett. , 4390(2001).[14] O. Derzhko, A. Honecker and J. Richter, Phys. Rev. B , 054403 (2009).[15] O. Derzhko, J. Richter, A. Honecker, M. Maksymenko,and R. Moessner, Phys. Rev. B , 014421 (2010). [16] O. Derzhko and J. Richter, Phys. Rev. B , 045152(2014).[17] R.R. Montenegro-Filho and M.D. Coutinho-Filho, Phys.Rev. B , 125117 (2006).[18] R.R. Montenegro-Filho, F.S. Matias and M.D. Coutinho-Filho, Phys. Rev. B , 035137 (2020).[19] Z. Gulácsi, A. Kampf and D. Vollhardt, Phys. Rev. Lett. , 026404 (2007); Prog. Theor. Phys. Supp. , 1(2008).[20] O. Rojas, S.M. de Souza and N.S. Ananikian, Phys. Rev.E , 061123 (2012).[21] O. Krupnitska, Phys. Rev. B , 064403 (2020).[22] J. Torrico, M. Rojas, M.S.S. Pereira, J. Strečka, and M.L.Lyra, Phys. Rev. B , 014428 (2016).[23] B. Nachtergaele, J.P. Solovej and J. Yngvason, Con-densed Matter Physics and Exactly Soluble Models , Se-lecta of E.H. Lieb, (Springer, Berlin-Heidelberg, 2004).[24] G. Beni and P. Pincus, Phys. Rev. B , 2963 (1974).[25] F. Mancini, Eur. Phys. J. B , 527 (2005).[26] F. Mancini and F.P. Mancini Phys. Rev. E , 061120(2008).[27] O. Rojas and S.M. de Souza, Phys. Lett. A , 1295(2011).[28] A.A. Lopes and R.G. Dias, Phys. Rev. B , 085124(2011).[29] B.M. Lisnii, Low Temp. Phys. , 296 (2011).[30] K.C. Rule, A.U.B. Wolter, S. Sullow, D.A. Tennant, A.Bruhl, S. Kohler, B. Wolf, M. Lang, and J. Schreuer,Phys. Rev. Lett. , 117202 (2008).[31] H. Kikuchi, Y. Fujii, M. Chiba, S. Mitsudo, T. Idehara,T. Tonegawa, K. Okamoto, T. Sakai, T. Kuwai, and H.Ohta, Phys. Rev. Lett. , 227201 (2005); H. Kikuchi, Y. Fujii, M. Chiba, S. Mitsudo, T. Idehara, T. Tonegawa, K.Okamoto, T. Sakai, T. Kuwai, K. Kindo, A. Matsuo, W.Higemoto, K. Nishiyama, M. Horvati?, and C. Bertheir,Prog. Theor. Phys. Suppl. , 1 (2005).[32] L. Canova, J. Strecka and M. Jascur, J. Phys.: Condens.Matter , 4967 (2006).[33] N.S. Ananikian, L.N. Ananikyan, L.A.Chakhmakhchyan, and O. Rojas, J. Phys.: Condens.Matter , 256001 (2012).[34] L. Chakhmakhchyan, N. Ananikian, L. Ananikyan, andC. Burdik, J. Phys: Conf. Series , 012022 (2012).[35] Sh.-S. Deng, Sh.-J. Gu, and H.-Q. Lin, Phys. Rev. B ,045103 (2006).[36] M.S.S. Pereira, F.A.B.F. de Moura and M.L. Lyra, Phys.Rev. B , 024402 (2008); Phys. Rev. B , 054427(2009).[37] H. Jeschke, I. Opahle, H. Kandpal, R. Valenti, H. Das,T. Saha-Dasgupta, O. Janson, H. Rosner, A. Bruhl, B.Wolf, M. Lang, J. Richter, S. Hu, X. Wang, R. Peters, T.Pruschke, and A. Honecker, Phys. Rev. Lett. , 217201(2011).[38] B. Gu and G. Su, Phys. Rev. B , 174437 (2007); O.Derzhko and J. Richter, Eur. Phys. J. B , 23 (2006).[39] M. Cui, N. Wang, S. Zhang, and Z. He, Crystal Growth& Design , 547 (2019).[40] A. Honecker and A. Lauchli, Phys. Rev. B , 174407(2001).[41] O. Rojas, Chi. Jour. Phys. , 157 (2021).[42] R.J. Baxter, Exactly Solved Models in Statistical Mechan-ics , (Academic Press, New York, 1982).[43] Y. Saito, Prog. Theor. Phys.59