Hubbard pair cluster with elastic interactions. Studies of thermal expansion, magnetostriction and electrostriction
AAccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 In Memory of Professor Leszek Wojtczak
Hubbard pair cluster with elastic interactions. Studies of thermalexpansion, magnetostriction and electrostriction
T. Balcerzak, K. Szałowski ∗ Department of Solid State Physics, Faculty of Physics and Applied Informatics,University of Ł´od´z, ul. Pomorska 149 / Abstract
The pair cluster (dimer) is studied within the framework of the extended Hubbard model andthe grand canonical ensemble. The elastic interatomic interactions and thermal vibrational en-ergy of the atoms are taken into account. The total grand potential is constructed, from whichthe equation of state is derived. In equilibrium state, the deformation of cluster size, as well asits derivatives, are studied as a function of the temperature and the external magnetic and elec-tric fields. In particular, the thermal expansion, magnetostriction and electrostriction e ff ects areexamined for arbitrary temperature, in a wide range of Hamiltonian parameters. Keywords:
Hubbard model, pair cluster, dimer, exact diagonalization, grand canonical ensemble,elastic interaction, magnetostriction, electrostriction
1. Introduction
The Hubbard model [1–4] plays an important role in contemporary solid state physics. Sinceits formulation, numerous applications of the model have been developed and the model itselfhas been extensively investigated as a fundamental, prototypic model for description of corre-lated fermions [5–8]. Many of the studies were concentrated on the typical three-dimensional,two-dimensional or one-dimensional Hubbard model (e.g. [9–11]). However, with an increase ∗ Corresponding author
Email addresses: [email protected] (T. Balcerzak), [email protected] (K.Szałowski) a r X i v : . [ c ond - m a t . s t a t - m ec h ] J un ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740of interest in studies of nanoclusters, which is stimulated by the development of nanotechnology,intensive investigations of the Hubbard model applied to the systems containing small number ofatoms have been performed. Some of them involve the studies of small clusters with either exactor close to exact methods [12–29]. Among such systems, a two-atomic cluster (dimer) plays aspecial role, being the smallest nanosystem where the Hubbard model can be adopted, and theexact analytical solutions can be found. As a consequence, numerous aspects of physics of Hub-bard dimers focused considerable attention in the literature [30–49]. These facts serve as a soundmotivation for further comprehensive study of the two-atomic Hubbard system aimed at its fulland exact thermodynamic description. In line with this trend, recently the Hubbard pair clusterembedded simultaneously in the external magnetic and electric fields has been studied by the ex-act analytical methods in the framework of the grand canonical ensemble. In this approach thesystem is open, so that it is able to exchange the electrons with its environment [41–43]. Thestatistical and thermodynamic properties, both magnetic and electric ones, have been extensivelystudied there.Despite these investigations, there is still room for extending the model of interest. Supple-menting purely electronic models by including further degrees of freedom can cause new kinds ofbehaviour to emerge and new control parameters can be gained this way. Purely electronic modelscan be supplemented, for example, with subsystems composed of localized spins (e.g. [50–57]),which are still exactly solvable. Immersing of such system in the external fields can give accessto unique properties (like, for example, manifestations of the magnetoelectric e ff ect [56]). Enrich-ing the purely electronic models is also a step towards more complete characterization of the realphysical systems. A good example can be inclusion of the Hamiltonian terms, which are respon-sible for coupling of the electronic degrees of freedom with the lattice. Such a study, aimed at thedescription of magnetomechanical and electromechanical properties of the Hubbard dimer, hasnot been undertaken yet. Therefore, the aim of the present paper is extension of the recent studies,basing on the formalism developed in Ref. [41], by including elastic interactions in the descriptionof the Hubbard dimer. Such interactions are especially important since they are responsible forstability of the dimer structure and allow to account for the energy of thermal vibrations. It can beexpected that after combining with the Hubbard Hamiltonian, the elastic interactions will lead to2ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740new phenomena, connecting the magnetic, electric and mechanical properties.The elastic interaction between the atoms is assumed here in the form of the Morse poten-tial [58, 59]. The thermal vibrational energy is obtained in the quasi-harmonic approximation,where the anharmonicity is taken into account by the Gr¨uneisen parameter [60, 61]. In turn, theHubbard pair Hamiltonian is taken in its extended form, with the Coulomb interaction betweenthe electrons on di ff erent atoms and the hopping integral depending on dimer size.On the basis of the above assumptions, the total grand thermodynamic potential is constructed,from which all statistical and thermodynamic properties can be obtained in a self-consistent man-ner. In particular, the deformation of interatomic distance (dimer length), magnetostriction andelectrostriction coe ffi cients are calculated for various temperatures, and the e ff ect of the externalmagnetic and electric fields on the mentioned quantities is also investigated. In addition, the chem-ical potential is studied in the presence of elastic interactions, showing a new behaviour in com-parison with the conventional Hubbard model. A special attention is paid to the low-temperatureregion, where the discontinuous quantum changes of these properties can be demonstrated.The paper is organized as follows: In the next Section the theoretical model is presented andthe basic formulae, important for numerical calculations, are derived. In the successive Section theselected results of calculations are illustrated in figures and discussed. The last Section is devotedto a brief summary of the results and final conclusions.
2. Theoretical model
In the present section, a step-by-step development of the theoretical model for the Hubbarddimer including the elastic and vibrational properties is presented and the solution of the model isdescribed.
Extended Hubbard Hamiltonian for the pair of atoms ( a , b ) embedded in the external fields isof the form: H a , b = − t (cid:88) σ = ↑ , ↓ (cid:16) c + a ,σ c b ,σ + c + b ,σ c a ,σ (cid:17) + U (cid:0) n a , ↑ n a , ↓ + n b , ↑ n b , ↓ (cid:1) , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 − H (cid:16) S za + S zb (cid:17) − V ( n a − n b ) + Wn a n b , (1)where t > U ≥ W stands for theCoulomb interaction between electrons localized on a and b atoms. In 1, H = − g µ B H z denotesthe magnitude of an external uniform magnetic field H z , whereas the electric field E is determinedby the electrostatic potentials V a = − V b = V applied to both atoms. Thus, E = V / ( | e | d ), where d is the interatomic distance and e is the electron charge. The quantity W in 1 is assumed in theclassical form: W = / (4 π(cid:15) ) e / d .The creation ( c + γ,σ ) and annihilation ( c γ,σ ) operators for site γ = a , b and spin state σ can beused to construct the occupation number operators n γ,σ for the electrons with given spin: n γ,σ = c + γ,σ c γ,σ . (2)The total occupation number operators at the site γ = a , b are then given by: n γ = n γ, ↑ + n γ, ↓ . (3)Moreover, the spin operators S z γ are defined as follows: S z γ = (cid:16) n γ, ↑ − n γ, ↓ (cid:17) / . (4)Since we intend to study the elastic properties of the model, the interatomic distance d can bepresented as follows: d = d (1 + ε ) , (5)where d is the equilibrium distance defined in the absence of the external fields ( H = , E = T = ε defines the temperature and fields-induced deformation.As it can be noted, some of the coe ffi cients in Hamiltonian 1, like V , W and t , are interatomicdistance-dependent. For instance, for the electric potential we have: V = V (1 + ε ) , (6)where V = E | e | d /
2. In turn, for the interatomic Coulomb interaction we can write: W = W + ε , (7)4ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740where W = e π(cid:15) d . Moreover, for the hopping integral we assume the power-law distance depen-dence, namely: t = t (cid:32) dd (cid:33) − n = t (1 + ε ) − n , (8)with a constant power index n , which can be used if the deformation ε is small.We treat the Hubbard pair as an open system, with a variable number of electrons, which canbe exchanged with the environment. It corresponds to the physical situation when, for instance,the cluster is interacting with the metallic substrate, which forms the electronic reservoir, and theaverage number of electrons in the system results from the thermodynamic equilibrium. In suchsmall system the relative fluctuations of the number of electrons, when related to its statisticalmean value, can be significant. Therefore, in this case a formalism of the grand canonical ensembleshould be used, which is more appropriate here than the canonical ensemble, with a fixed numberof electrons.In the framework of the grand canonical ensemble, the Hamiltonian 1 should be extended byadding the term − µ ( n a + n b ), where µ is the chemical potential. For such extended Hamiltonian,the diagonalization procedure can be performed exactly using the method outlined in Ref. [41].It is worth noticing that by including the new term, Wn a n b , the complexity of diagonalizationprocedure will not increase in comparison to that presented in Ref. [41]. We note that n a n b is adiagonal 16 ×
16 matrix and corresponds to 16 possible states. Its diagonal elements are: a , = a , = a , = a , = a , = a , = a , = a , = a , = a , = a , = a , = a , = a , = a , = a , =
4. Such a matrix, multiplied by W / t , should be simply addedto that presented in Ref. [41] (Eq. B.5 therein).As a result of diagonalization, the grand thermodynamic potential Ω a , b for the Hubbard paircan be found from the general formula: Ω a , b = − k B T ln { Tr a , b exp[ − β (cid:0) H a , b − µ ( n a + n b ) (cid:1) ] } . (9)If we denote the mean number of the electrons per atom by x , then the following relation issatisfied: 2 x = (cid:104) n a (cid:105) + (cid:104) n b (cid:105) = − (cid:32) ∂ Ω a , b ∂µ (cid:33) T , H , E . (10)5ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740In Eq. 10, (cid:104)· · ·(cid:105) means the statistical average, which is carried out over all 16 states, and is depen-dent on the temperature and external fields. From Eq. 10 the chemical potential µ can be found,provided that x is established, and then the expression 9 is complete. Irrespective of the Hubbard contribution, the interatomic static interaction of the core electronscan be taken into account in the form of the Morse potential energy [58]: V M ( r ) = D (cid:104) − e − α ( r − r ) / r (cid:105) . (11)In Eq. 11, D is the potential depth, parameter α controls the potential width and its asymmetry,whereas r is the equilibrium distance (potential minimum position) for this interaction at T = ε , and its value can be normalized by therequirement that V M ( ε = =
0. In this way we obtain the elastic energy in the form of: V M ( ε ) = D (cid:34) − e − α (cid:18) d r (1 + ε ) − (cid:19) (cid:35) − D (cid:34) − e − α (cid:18) d r − (cid:19) (cid:35) . (12)The dimensionless parameter d / r can be found from the minimum condition of the total energy,taken at T = In order to take into account the thermal vibrations, the two-atomic system can be treated asa quantum oscillator having reduced mass and described in the centre of mass coordinate system.This situation corresponds, for instance, to interaction of the system with the solid substrate, whichforms a phononic bath. For the sake of simplicity we adopt the quasi-harmonic Einstein model, inwhich the anharmonicity is taken into account by the Gr¨uneisen parameter Γ [61]. The vibrationalfree energy of the Einstein oscillator is then given by: E ω = k B T ln (cid:34) (cid:32) β (cid:126) ω (cid:33)(cid:35) . (13)The frequency of oscillations, ω , is dependent on the interatomic distance and can be expressed asa function of the deformation ε as follows: ω = ω (1 + ε ) Γ , (14)6ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740where the Gr¨uneisen parameter Γ satisfies the relationship: Γ = − d ω (cid:32) ∂ω∂ d (cid:33) . (15)It is known that Γ can be determined theoretically for several interatomic potentials and di ff er-ent dimensionalities of the system. For instance, for the Morse potential and 1D system Γ = (3 / α [62]. This convenient relationship will be adopted in further calculations.Thus, the vibrational energy can be finally presented as: E ω = k B T ln (cid:34) (cid:32) t k B T (cid:126) ω t
12 (1 + ε ) Γ (cid:33)(cid:35) , (16)where k B T / t is the dimensionless temperature, and (cid:126) ω / t is the dimensionless characteristicenergy parameter depending on the atomic mass. The total grand potential, Ω , of the system in question can be constructed as a sum of theHubbard grand potential 9, elastic (static) energy 12 and vibrational Helmholtz free-energy 16: Ω = Ω a , b + V M ( ε ) + E ω . (17)This grand potential should be minimized first with respect to the deformation ε by imposing thenecessary equilibrium condition: (cid:32) ∂ Ω ∂ε (cid:33) T , H , E = . (18)Di ff erentiation of Ω will then lead to three forces arising out of the three terms in r.h.s. of Eq. 17.The Hubbard force, F a , b , is connected with the change of the Hubbard grand potential underthe deformation: F a , b = − d (cid:32) ∂ Ω a , b ∂ε (cid:33) . (19)With the help of Eq. 9 we get: − d F a , b = ∂∂ε (cid:10) H a , b (cid:11) − x ∂µ∂ε , (20)where (cid:10) H a , b (cid:11) is the statistical average of the Hubbard Hamiltonian. The energy (cid:10) H a , b (cid:11) can be foundexactly after the process of diagonalization. Then, the derivative ∂ (cid:10) H a , b (cid:11) /∂ε can be calculated7ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740analytically with the help of Eqs. 6–8. In turn, the derivative ∂µ/∂ε can be, in principle, calculatednumerically for the arbitrary concentration (0 ≤ x ≤
2) of Hubbard electrons. For such purposeEq. 10 for the chemical potential can be used. Interestingly, for x =
1, i.e., for half-filled energylevels, the analytical formula for the chemical potential µ can be found, namely: µ = W + ε + U , (21)which e ff ectively facilitates the calculations of ∂µ/∂ε . As a result, the Hubbard force for x = F a , b = − n (1 + ε ) n + t d (cid:88) σ = ↑ , ↓ (cid:16)(cid:68) c + a ,σ c b ,σ (cid:69) + (cid:68) c + b ,σ c a ,σ (cid:69)(cid:17) + V d ( (cid:104) n a (cid:105) − (cid:104) n b (cid:105) ) − W d (1 + ε ) ( (cid:104) n a (cid:105) + (cid:104) n b (cid:105) − (cid:104) n a n b (cid:105) ) . (22)The statistical averages in Eq. 22 can be calculated after performing the diagonalization of theHamiltonian. The examples of similar calculations (in the absence of elastic interactions) havealready been presented in Ref. [42].As a next step, the elastic force, F ε , is found from the formula: F ε = − d (cid:32) ∂ V M ( ε ) ∂ε (cid:33) , (23)where V M ( ε ) is given by Eq. 12. Hence, we obtain: F ε = − D α r (cid:34) − e − α (cid:18) d r (1 + ε ) − (cid:19) (cid:35) e − α (cid:18) d r (1 + ε ) − (cid:19) . (24)The last force, F ω , which is of vibrational origin, is defined as: F ω = − d (cid:32) ∂ E ω ∂ε (cid:33) , (25)where E ω is given by Eq. 16. From Eqs. ?? nosort]eq25,eq16 we obtain: F ω = Γ + ε ) Γ+ (cid:126) ω d tanh − (cid:32) t k B T (cid:126) ω t
12 (1 + ε ) Γ (cid:33) . (26)It follows from the numerical calculations that the Hubbard force, F a , b , is compressive (negative)around ε =
0. On the other hand, the elastic force, F ε , is compressive for d > r and expansive8ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740(positive) for d < r , whereas the vibrational force, F ω , is always of expansive type. The necessarycondition for thermodynamic equilibrium requires the balance of all these forces, namely: F a , b + F ε + F ω = . (27)Expression 27 is the EOS for the Hubbard dimer with the elastic interactions and thermal vibra-tions taken into account. It enables calculation of spontaneous deformation, ε , of the dimer lengthfor given values of the intensive parameters: T , H and E , whereas the system is in stable equilib-rium. At first, in the absence of deformation (i.e., when ε = T = H = E = d / r should be found from Eq. 27, giving the equilibrium distance d at theground state.The numerical calculations based on the EOS 27 will be presented in the next Section.
3. Numerical results and discussion
The numerical results will be presented for the half-filling case, i.e., when the average num-ber of electrons per atom is x =
1. In the Hubbard Hamiltonian we select the exponent n in thehopping integral (8) equal to n = d . It is in analogy with the possible dependence of the exchange integral vs. the distance [64]and refers to the attractive part of the Lennard-Jones potential. A more accurate description ofthe hopping integrals dependence on the bond length can be given by the Goodwin, Skinner andPettifor function [67]. The interatomic Coulomb interaction parameter W d r t = e π(cid:15) r t can be es-timated for the realistic values of r ≈ . t ≈ . ≈ t and the interatomic equilibrium distance d is related to the r constant. In this convention, we chose the value D / t = α = (cid:126) ω / t = .
5, whereas the Gr¨uneisen parameteris
Γ = (3 / α [62]. Such a set of parameter values, although not describing any specific mate-rial, seems physically acceptable and can be used to demonstrate how the formalism works whenthe numerical calculations are performed. We are aware that the choice of these values can be9ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740discussed and modified according to the physical situation. Unfortunately, because of the limitedscope of the paper we are not able to present the influence of all these parameters on the numericalresults. However, it can be mentioned that, for instance, some change of the asymmetry parameter α , or the potential depth D , may have remarkable influence on the thermodynamic properties, as ithas been demonstrated in Ref. [65].For the above set of parameters and for U / t =
0, from Eq. 27 we obtain the equilibriumdistance d / r = . T = H = E =
0, whereasdeformation is ε = U / t = , , d / r = . µ/ t , for U / t = , , , µ/ t = . d is not a monotonousfunction of the parameter U . Moreover, the chemical potential µ is not equal to U /
2, as it happensin the case when the interatomic Coulomb repulsion, as well as the elastic interactions, are ne-glected [5, 41]. The values of µ obtained here are in agreement with Eq. 21. The formula 21 alsopredicts that out of the ground state, i.e., for T >
0, and in the possible presence of the externalfields, the chemical potential is no longer constant, but it should be a function of T , H and E , viaits dependency on the deformation ε . This new behaviour is illustrated in Fig. 1-3. For instance,in 1 the di ff erence ( µ − U / / t is shown as a function of the dimensionless temperature k B T / t for several values of U / t , when the external fields are absent. In general, the functions decreasewith the increasing temperature and, at the same time, the dependency of ( µ − U /
2) on U is notmonotonous. In 2 also the decreasing tendency of ( µ − U / / t as a function of dimensionlessmagnetic field H / t is shown for the same values of U / t . In this case the electric field is equal to E =
0, and the moderate temperature, k B T / t = .
5, is assumed. In 3 the chemical potential is plot-ted for k B T / t = . H =
0, as a function of the dimensionless electric field E | e | r / t . Apartfrom the decreasing tendency of ( µ − U / / t vs. E , which is visible for small U , the increasingtendency of the chemical potential for U / t = U / t =
10 is worth noticing.4 presents relative deformation ε of the dimer size vs. magnetic field H / t at very low temper-ature, k B T / t = . E =
0. Several values of U / t are examined. As we know from10ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740our previous studies [41–43], in the low-temperature region the quantum phase transitions shouldoccur at the magnetic (or electric) critical fields. Such transitions are of the first order, as regardsthe behaviour of total dimer magnetization, which results from the change from the singlet state,where (cid:68) S za + S zb (cid:69) =
0, to a triplet one, where (cid:68) S za + S zb (cid:69) =
1. The rapid change of magnetizationis accompanied by instantaneous changes of the remaining thermodynamic quantities. In 4 thechanges of deformation ε are illustrated at the transition from the singlet state, where ε =
0, tothe triplet (ferromagnetically saturated) state. It should be mentioned that due to the presence ofinteratomic Coulomb interaction W , the critical fields H c found here have di ff erent values in com-parison with those from our previous papers [41–43]. Unfortunately, there is no analytical formulafor these critical fields, and they can be calculated numerically only by equating the total grandpotentials of neighbouring phases.In 5 the electrostriction coe ffi cient, ν T , H = d (cid:16) ∂ d ∂ E (cid:17) T , H = + ε (cid:16) ∂ε∂ E (cid:17) T , H , is presented in dimension-less units, ν T , H t | e | r , as a function of the magnetic field H / t . Several values of the parameter U areexamined at the temperature equal to k B T / t = . E | e | r / t = .
5. As it is seen, the characteristic steps corresponding to the quantum phasetransitions occur for the electrostriction coe ffi cient. The critical magnetic fields H c , where thejumps occur, decrease with an increase in the Coulombic repulsion energy U , which e ff ect is anal-ogous to that noted in the previous figure. Moreover, for higher U values, the existence of negativeelectrostriction is demonstrated. It can also be noted that the electrostriction coe ffi cient takes zerovalue in the triplet (ferromagnetically saturated) state, because then the electric field is not able tochange the charge distribution [43], whereas the electric potential is compensated by the chemicalpotential. As a result, the total energy is not influenced by the electric field and length of the dimerremains unchanged.In the next figure, 6, the deformation ε is presented at the temperature k B T / t = .
001 andthe magnetic field H / t = . E | e | r / t . Again, thecharacteristic jumps of ε are shown for various curves corresponding to di ff erent values of U .The critical electric fields E c , where the phase transition from triplet to singlet state takes place,increase with an increase in U . Below the critical electric field, i.e., in triplet state, the deformation ε is constant vs. E , which results in null electrostriction e ff ect. This fact supports the discussion11ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740presented in the preceding paragraph and is in agreement with the general behaviour of the criticalfields for di ff erent U calculated for the ground state in the absence of the elastic interactions [42](Fig. 1 therein). On the basis of Fig. 3-6, it may be interesting to notice that instead of speakingabout the critical fields E c (or H c ) for given U -parameter and the fields H (or E ), we can equallywell speak about the critical Coulomb potential U c for given H and E .7 is devoted to illustration of the magnetostriction behaviour in two temperatures: k B T / t = .
05 and k B T / t = . H / t = .
8. The magne-tostriction coe ffi cient is defined as λ T , E = d (cid:16) ∂ d ∂ H (cid:17) T , E = + ε (cid:16) ∂ε∂ H (cid:17) T , E and is presented in dimensionlessunits, λ T , E t , vs. dimensionless electric field E | e | r / t . It is seen that the magnetostriction coef-ficient in lower temperature shows the sharp peaks, whereas in higher temperature the maximaare lower and much more di ff used. These maxima appear approximately at the critical electricfields E c and are placed at di ff erent positions for the curves plotted for various parameters U . Inparticular, when U increases, the critical fields are shifted towards higher values of T . For T → ffi cient is zero everywhere aside from the singularities, because in thelow-temperature region the magnetic field is not able to change the magnetization, both in singletand in triplet state, except for the phase transition point [42]. As a result, the magnetic energy isconstant aside from the singularities and changes of the magnetic field do not influence the dimersize. This conclusion is also in agreement with 4 where, apart from the phase transition points,the plots of ε vs. H are flat, which corresponds to the null magnetostriction e ff ect. When thetemperature increases, the quantum phase transitions become di ff used and, as a result, the magne-tostriction coe ffi cient presents broadened maxima with reduced height. The positions of the peaksare quite stable and are only slightly shifted towards larger electric fields. At the same time, theheight of the peaks strongly depends on the Hubbard parameter U , showing a decrease of maximalmagnetostriction coe ffi cient when U increases.The coe ffi cient of the linear thermal expansion is presented in 8 as a function of the dimension-less temperature k B T / t . This coe ffi cient is defined as α H , E = + ε (cid:16) ∂ε∂ T (cid:17) H , E , where the derivative over T is taken at constant fields H and E . In 8 the dimensionless quantity α H , E t / k B is plotted, whereasthe fields H and E are absent. The most interesting e ff ect illustrated in this figure is the maximumof the thermal expansion and its evolution when the Coulombic repulsion energy U increases. This12ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740maximum is initially small for the curve corresponding to U = U / t = U / t =
2. Then, for U / t = U / t = U / t =
10. At the sametime, the maximum is shifted towards lower temperatures again. The occurrence of this maximumis not connected with the phase transition, since the system still remains in the singlet state. Itrather results from a complicated interplay and competition between di ff erent contributions to thetotal grand potential, which depend on the interatomic distance.In order to outline the origin of this maximum, it can be noticed that an increase in U -potentialbrings on the charge separation, and as a result the occupation correlation on di ff erent atoms, (cid:104) n a n b (cid:105) , increases. For instance, in 8 for the fixed temperature k B T / t = .
5, and U / t =
2, 6 and10, the occupation correlations amount to: (cid:104) n a n b (cid:105) = . W coe ffi cient (Eq. 7) mustdecrease, which means that deformation ε takes higher values there. One can also note that, if d / r <
1, some increase in ε helps to lower the Morse potential energy. The resulting deformationis obtained as the solution of EOS (Eq. 27). For the above example ( k B T / t = . U / t = ε = . U , the deformation is a non-linearly increasing function of temperature.In this case the increase in ε is also connected with the energy compensation mechanism, and it ismainly stimulated by the fact that in the singlet (paramagnetic) phase the increase in temperatureresults in a nonlinear increase of the mean hopping energy (see [42] and Fig. 5 therein). Then, thederivative ∂ε/∂ T can reveal a non-monotonous behaviour vs. temperature and the maximum of α H , E appears as a consequence. The non-linear changes of ε are more pronounced in the interme-diate range of temperatures and for larger U -parameters, where the competition between variouspressure components in EOS is more dynamic. In our opinion, the revealed maximum of α H , E is,to some extent, analogous to the paramagnetic maximum of the specific heat, observed in localizedspin systems and known as the Schottky anomaly. In that e ff ect the magnetic energy (or entropy)also shows a non-linear increase vs. temperature.It can also be mentioned that for T → ffi cient tends to zero, whichis a proper thermodynamic behaviour. On the other hand, in the high temperature region, the13ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740anharmonic vibrations play a dominant role and, as a result, an increase in α H , E occurs there in amore regular way.The next figures (Fig. 9-14) are prepared in form of contour graphs, presenting in the multiplecharts the deformation (in parts a), magnetostriction coe ffi cient (in parts b) and electrostrictioncoe ffi cient (in parts c). The isolines of the dimensionless quantities, ε , 1000 × λ T , E t , and 1000 × ν T , H t / ( | e | r ), are plotted as a function of the external fields, in the plane ( H / t , E | e | r / t ). Twotemperatures, which are far from the ground state: k B T / t = . k B T / t = . U , three values are selected: U / t = ?? , U / t = ?? , and U / t =
10 for ?? . We see that both T and U play a significant role in the behaviour of the quantities presentedin multiple contour plots. The shapes of the isolines, defining the areas with similar values ofthe presented quantities, change dynamically between di ff erent plots. In particular, analysing thedeformation, it can be noted that areas with the smallest values of ε , denoted by blue shades andsituated in the left parts of the charts, where rather small magnetic fields dominate, are shiftedquite irregularly along the vertical (electric field) axis for di ff erent T and U . At the same time,the areas with highest values of ε , denoted by red colour and other warm shades, are expanding asthe temperature increases for each value of U . For instance, analysing the isolines for ε = U , we see that the areas with ε ≥ T increases. Onthe other hand, all isolines in 14(a) correspond to larger values of ε than isolines in 13(a). The redareas are situated at high electric fields when U / t = U / t =
5, however, when U / t =
10, thearea with large ε extends over all values of the electric fields down to zero, provided the magneticfield is not very small. It should be stressed that deformation ε is positive for all these cases. Theincreasing values of ε when the temperature increases indicate that also the thermal expansioncoe ffi cient should be positive.Regarding the magnetostriction, its value is mostly positive, however, some areas with thenegative magnetostriction can be found, namely for high magnetic and electric fields. Interestingly,the area with negative magnetostriction expands with the increasing temperature. On the otherhand, the dependency on U is not monotonous: although the area with negative magnetostriction14ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740is larger for U / t = U / t =
2, however, for U / t =
10 this range becomes stronglyreduced and is only visible at higher temperature ( k B T / t = . U / t = U / t =
5, however, it is shifted towards large electric fields if U / t =
10. The magnetic fieldsare moderate in all cases of high magnetostriction ranges, being from the range of (2 (cid:46) H / t (cid:46) U / t =
2, and (0 . (cid:46) H / t (cid:46)
2) for higher U / t . The increase in temperature results also insome expansion of the area with strong positive magnetostriction. It can also be noted that isolinescorresponding to zero magnetostriction present an interesting behaviour, namely, the increase inelectric field causes the decrease in the magnetic field along the isoline, so that the slope of theline as a function of H is negative.As far as the electrostriction coe ffi cient is concerned, some areas with the negative values ofthis quantity are also observed in all the cases. For small values of U these areas lie at low electricfields. In particular, for U / t =
2, rather high magnetic fields ( H / t (cid:38)
3) are involved, whereasfor U / t = H / t =
0, as in the case of 12. On the other hand, for U / t =
10 the area with negative elec-trostriction occurs for considerably higher electric fields, for instance, in the order of E | e | r / t ≈ k B T / t = . k B T / t = .
0, makes the area of negative electrostriction larger, but, at the same time, its amplitudedecreases approximately by an order of magnitude. On the other hand, the areas corresponding tothe highest positive electrostriction are mainly situated at very large electric fields. An exception is9 where, for U / t = . k B T / t = .
5, the region of positive electrostriction coe ffi cient extendsdown to very small electric fields.In general, the isolines separating the areas of positive and negative electrostriction are moredensely packed for k B T / t = . k B T / t = .
0. The oblique arrangement of these lines at k B T / t = . ffi cient when crossingthem horizontally by an increase in the magnetic field. The isolines corresponding to zero elec-trostriction e ff ect can reveal either positive or negative slope as a function of the magnetic field H .The positive slope occurs, first of all, in the range of small magnetic fields, whereas for larger H values the negative slopes are dominating (see, for instance ?? for U / t = , 121740 (2019),DOI:10.1016 / j.physa.2019.121740the behaviour of the magnetostriction, where the slopes of the null magnetostriction coe ffi cientline were always negative.
4. Summary and conclusion
In the paper the total grand potential has been constructed for the dimer system, in which themagnetic, elastic and thermal vibration energies play important roles. The similar attempts atthe total Gibbs energy construction, i.e., consisting in summation of various energy contributions,have been undertaken, so far, for the bulk materials [61, 63–66]. According to our knowledge, forthe Hubbard dimer such method has not been presented previously.It should be mentioned that also another approach to the problem can be proposed, namely,including the term representing electron–phonon interaction into the common Hamiltonian, andapplication of canonical transformation of Lang–Firsov type [68–72]. Then, the procedure is fol-lowed by using the 2nd order perturbation calculus, or, equivalently, by the Schrie ff er–Wol ff trans-formation [69], which leads to the renormalization of Hamiltonian parameters. Such a method ismuch more complicated and, in our opinion, is therefore less convenient for the complete studiesof various thermodynamic properties. Moreover, the Lang–Firsov type transformation has beenproposed for the Hamiltonians describing the coupling of electrons with the harmonic oscilla-tors [70], whereas in our approach the anharmonic vibrations are taken into account via Gr¨uneisenparameter Γ . It is known that the anharmonicity plays an important role, for instance, in the studiesof thermal expansion e ff ect.As a result of the present approach, the EOS (Eq. 27) has been derived for the system inquestion. This equation yields stable solutions for the interatomic distance deformation, ε , cor-responding to the thermodynamic equilibrium state obtained at arbitrary temperature T and theexternal fields H and E . Investigation of this deformation leads to the determination of its deriva-tive quantities such as the thermal expansion or magneto- and electrostriction coe ffi cients.It has been shown in the paper that for the half-filling conditions ( x =
1) and for the elasticinteractions included, the chemical potential µ is not a constant but it depends, via the quantity ε , on the temperature and external fields (Fig. 1-3). For x = µ , 121740 (2019),DOI:10.1016 / j.physa.2019.121740has been presented (Eq. 21). It should be stressed that the method used here can be equally wellapplied to arbitrary average concentration of the electrons, i.e., when 0 ≤ x ≤ ε have been found at the phase transition points.These discontinuous changes are also accompanied by discontinuities of such quantities as themagnetostriction and electrostriction coe ffi cients (Fig. 4-7). It is known that with an increase intemperature the quantum phase transitions become di ff used, so then, the most spectacular e ff ectsof discontinuity occur at T → ff ect. In our opinion, an interesting finding concerns the broad maxima of the thermalexpansion coe ffi cient, found at some temperatures and illustrated in 8. It has been shown that thepositions of maxmima are strongly dependent on the values of Hubbard U parameter.In the high-temperature region, where the discontinuities of thermodynamic quantities do notoccur, the contour charts yield a comprehensive description of the interatomic distance deforma-tion and its changes vs. external fields. The contour charts (Fig. 9-14) have been used to illustratesimultaneously ε , as well as the magnetostriction and electrostriction coe ffi cients vs. H and E , forsome selection of temperatures T and parameters U . In particular, a distinctive result is an identi-fication of the ( H , E )-regions where the negative values of these coe ffi cients do appear. It has beenfound that both temperature T and U -parameter have important influence on the localization andarea of these regions, as well as on the values of ε , λ T , E , and ν T , H themselves.To summarize, by taking into consideration the elastic interactions and thermal vibrations,a new, more complete thermodynamic description of the Hubbard dimer has been achieved. Itshould be stressed that diagonalization of the Hubbard Hamiltonian has been performed exactlyby means of the analytical methods [41]. It can be expected that selection of other interatomicpotentials would change our results quantitatively; however, the case which we consider in thepresent paper captures the essence of the phenomena emerging due to coupling between purelyelectronic and elastic properties.It should also be pointed out that having constructed the total grand potential, all thermody-namic and statistical properties of the system can be studied. For instance, this allows the self-17ccepted manuscript. The final version was published in:Physica A: Statistical Mechanics and its Applications , 121740 (2019),DOI:10.1016 / j.physa.2019.121740consistent calculation of the magnetization, electric polarization, magnetic correlation functions,occupation correlations, susceptibility (both magnetic and electric one), entropy and specific heat,in the presence of the elastic interactions. References [1] P. W. Anderson, New Approach to the Theory of Superexchange Interactions, Phys. Rev. 115 (1) (1959) 2–13. doi:10.1103/PhysRev.115.2 .[2] J. Hubbard, Electron Correlations in Narrow Energy Bands, Proceedings of the Royal Society of London A:Mathematical, Physical and Engineering Sciences 276 (1365) (1963) 238–257. doi:10.1098/rspa.1963.0204 .[3] M. C. Gutzwiller, E ff ect of Correlation on the Ferromagnetism of Transition Metals, Phys. Rev. Lett. 10 (5)(1963) 159–162. doi:10.1103/PhysRevLett.10.159 .[4] J. Kanamori, Electron Correlation and Ferromagnetism of Transition Metals, Progress of Theoretical Physics30 (3) (1963) 275–289. doi:10.1143/PTP.30.275 .[5] W. Nolting, A. Ramakanth, Quantum Theory of Magnetism, Springer-Verlag, Berlin, 2009.[6] K. Yosida, Theory of Magnetism, Springer-Verlag, Berlin, 1998.[7] H. Tasaki, The Hubbard model - an introduction and selected rigorous results, Journal of Physics: CondensedMatter 10 (20) (1998) 4353. doi:10.1088/0953-8984/10/20/004 .[8] A. Mielke, The Hubbard Model and its Properties, in: E. Pavarini, E. Koch, P. Coleman (Eds.), Many-BodyPhysics: From Kondo to Hubbard. Modeling and Simulation, Vol. 5, Verlag des Forschungszentrum J¨ulich,J¨ulich, 2015. hdl.handle.net/2128/9255 .[9] B. S. Shastry, Exact Integrability of the One-Dimensional Hubbard Model, Phys. Rev. Lett. 56 (23) (1986)2453–2455. doi:10.1103/PhysRevLett.56.2453 .[10] J. E. Hirsch, S. Tang, Antiferromagnetism in the Two-Dimensional Hubbard Model, Phys. Rev. Lett. 62 (5)(1989) 591–594. doi:10.1103/PhysRevLett.62.591 .[11] S. Fuchs, E. Gull, L. Pollet, E. Burovski, E. Kozik, T. Pruschke, M. Troyer, Thermodynamics of the 3D Hub-bard Model on Approaching the N´eel Transition, Phys. Rev. Lett. 106 (3) (2011) 030401. doi:10.1103/PhysRevLett.106.030401 .[12] Y. Ohta, K. Tsutsui, W. Koshibae, S. Maekawa, Exact-diagonalization study of the Hubbard model with nearest-neighbor repulsion, Physical Review B 50 (18) (1994) 13594–13602. doi:10.1103/PhysRevB.50.13594 .[13] M. A. Ojeda-L´opez, J. Dorantes-D´avila, G. M. Pastor, Magnetic properties of Hubbard clusters: Noncollinearspins in N-atom clusters, Journal of Applied Physics 81 (8) (1997) 4170–4172. doi:10.1063/1.365128 .[14] C. Noce, M. Cuoco, A. Romano, Thermodynamical properties of the Hubbard model on finite-size clusters,Physica C: Superconductivity 282-287 (1997) 1705–1706. doi:10.1016/S0921-4534(97)00972-6 . , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 [15] M. A. Ojeda, J. Dorantes-D´avila, G. M. Pastor, Noncollinear cluster magnetism in the framework of the Hubbardmodel, Physical Review B 60 (8) (1999) 6121–6130. doi:10.1103/PhysRevB.60.6121 .[16] F. Becca, A. Parola, S. Sorella, Ground-state properties of the Hubbard model by Lanczos diagonalizations,Physical Review B 61 (24) (2000) R16287–R16290. doi:10.1103/PhysRevB.61.R16287 .[17] J. L. Ricardo-Ch´avez, F. Lopez-Ur´ıas, G. M. Pastor, Thermal Properties of Magnetic Clusters, in: J. L. Mor´an-L´opez (Ed.), Physics of Low Dimensional Systems, Kluwer Academic / Plenum Publishers, New York, 2001, pp.23–32. doi:10.1007/0-306-47111-6_3 .[18] Y. Hancock, A. E. Smith, Application of the Hubbard model with external field to quasi-zero dimensional, in-homogeneous systems, Physics Letters A 300 (4) (2002) 491–499. doi:10.1016/S0375-9601(02)00854-X .[19] Y. Hancock, Quasi-zero-dimensional quantum spin-switching system, Physical Review B 71 (22) (2005)224428. doi:10.1103/PhysRevB.71.224428 .[20] R. Schumann, Analytical solution of extended Hubbard models on three- and four-site clusters, Physica C:Superconductivity 460-462 (2007) 1015–1017. doi:10.1016/j.physc.2007.03.203 .[21] R. Schumann, Rigorous solution of a Hubbard model extended by nearest-neighbour Coulomb and exchangeinteraction on a triangle and tetrahedron, Annalen der Physik 17 (4) (2008) 221–259. doi:10.1002/andp.200710281 .[22] S. X. Yang, H. Fotso, J. Liu, T. A. Maier, K. Tomko, E. F. D’Azevedo, R. T. Scalettar, T. Pruschke, M. Jarrell,Parquet approximation for the 4x4 Hubbard cluster, Physical Review E 80 (4) (2009) 046706. doi:10.1103/PhysRevE.80.046706 .[23] R. Schumann, D. Zwicker, The Hubbard model extended by nearest-neighbor Coulomb and exchange interactionon a cubic cluster - rigorous and exact results, Annalen der Physik 522 (6) (2010) 419–439. doi:10.1002/andp.201010452 .[24] H. Feldner, Z. Y. Meng, A. Honecker, D. Cabra, S. Wessel, F. F. Assaad, Magnetism of finite graphene samples:Mean-field theory compared with exact diagonalization and quantum Monte Carlo simulations, Phys. Rev. B81 (11) (2010) 115416. doi:10.1103/PhysRevB.81.115416 .[25] M. Y. Ovchinnikova, Spin Structures of the Hubbard Clusters and the Pseudogap, Journal of Superconductivityand Novel Magnetism 26 (9) (2013) 2845–2846. doi:10.1007/s10948-013-2218-0 .[26] Y. Hancock, A family of spin-switching, inhomogeneous Hubbard chains, Physica E: Low-dimensional Systemsand Nanostructures 56 (2014) 141–150. doi:10.1016/j.physe.2013.08.021 .[27] T. X. R. Souza, C. A. Macedo, Ferromagnetic Ground States in Face-Centered Cubic Hubbard Clusters, PLOSONE 11 (9) (2016) e0161549. doi:10.1371/journal.pone.0161549 .[28] K. Szałowski, T. Balcerzak, M. Jaˇsˇcur, A. Bob´ak, M. ˇZukoviˇc, Exact Diagonalization Study of an ExtendedHubbard Model for a Cubic Cluster at Quarter Filling, Acta Physica Polonica A 131 (4) (2017) 1012–1014. doi:10.12693/APhysPolA.131.1012 . , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 [29] K. Szałowski, T. Balcerzak, Electrocaloric e ff ect in cubic Hubbard nanoclusters, Scientific Reports 8 (1) (2018)5116. doi:10.1038/s41598-018-23443-x .[30] B. Alvarez-Fern´andez, J. A. Blanco, The Hubbard model for the hydrogen molecule, European Journal ofPhysics 23 (1) (2002) 11. doi:10.1088/0143-0807/23/1/302 .[31] A. B. Harris, R. V. Lange, Single-Particle Excitations in Narrow Energy Bands, Phys. Rev. 157 (2) (1967)295–314. doi:10.1103/PhysRev.157.295 .[32] R. Wortis, M. P. Kennett, Local integrals of motion in the two-site Anderson-Hubbard model, Journal of Physics:Condensed Matter 29 (40) (2017) 405602. doi:10.1088/1361-648X/aa818e .[33] Y.-C. Cheng, K.-C. Liu, Thermodynamic Properties of a Two-Site Hubbard Hamiltonian: Half-Filled-BandCase, Chinese Journal of Physics 14 (3) (1976) 325–328.[34] S. Longhi, G. Della Valle, V. Foglietti, Classical realization of two-site Fermi-Hubbard systems, Phys. Rev. B84 (3) (2011) 033102. doi:10.1103/PhysRevB.84.033102 .[35] R. C. Juliano, A. S. de Arruda, L. Craco, Coexistence and competition of on-site and intersite Coulomb interac-tions in Mott-molecular-dimers, Solid State Communications 227 (2016) 51 – 55. doi:http://dx.doi.org/10.1016/j.ssc.2015.11.021 .[36] M. E. Kozlov, V. A. Ivanov, K. Yakushi, Development of a two-site Hubbard model for analysis of the electron-molecular vibration coupling in organic charge-transfer salts, Physics Letters A 214 (3) (1996) 167 – 174. doi:http://dx.doi.org/10.1016/0375-9601(96)00113-2 .[37] A. V. Silant’ev, A Dimer in the Extended Hubbard Model, Russian Physics Journal 57 (11) (2015) 1491–1502. doi:10.1007/s11182-015-0406-z .[38] H. Hasegawa, Nonextensive thermodynamics of the two-site Hubbard model, Physica A: Statistical Mechanicsand its Applications 351 (2-4) (2005) 273 – 293. doi:http://dx.doi.org/10.1016/j.physa.2005.01.025 .[39] H. Hasegawa, Thermal entanglement of Hubbard dimers in the nonextensive statistics, Physica A: StatisticalMechanics and its Applications 390 (8) (2011) 1486 – 1503. doi:http://dx.doi.org/10.1016/j.physa.2010.12.033 .[40] J. Spałek, A. M. Ole´s, K. A. Chao, Thermodynamic properties of a two-site Hubbard model with orbitaldegeneracy, Physica A: Statistical Mechanics and its Applications 97 (3) (1979) 552 – 564. doi:http://dx.doi.org/10.1016/0378-4371(79)90095-5 .[41] T. Balcerzak, K. Szałowski, Hubbard pair cluster in the external fields. Studies of the chemical potential, PhysicaA: Statistical Mechanics and its Applications 468 (2017) 252–266. doi:10.1016/j.physa.2016.11.004 .[42] T. Balcerzak, K. Szałowski, Hubbard pair cluster in the external fields. Studies of the magnetic properties,Physica A: Statistical Mechanics and its Applications 499 (2018) 395–406. doi:10.1016/j.physa.2018.02.017 . , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 [43] T. Balcerzak, K. Szałowski, Hubbard pair cluster in the external fields. Studies of the polarization and suscepti-bility, Physica A: Statistical Mechanics and its Applications 512 (2018) 1069–1084. doi:10.1016/j.physa.2018.08.152 .[44] C. A. Ullrich, Density-functional theory for systems with noncollinear spin: Orbital-dependent exchange-correlation functionals and their application to the Hubbard dimer, Physical Review B 98 (3) (2018) 035140. doi:10.1103/PhysRevB.98.035140 .[45] F. Sagredo, K. Burke, Accurate double excitations from ensemble density functional calculations, The Journalof Chemical Physics 149 (13) (2018) 134103. doi:10.1063/1.5043411 .[46] D. J. Carrascal, J. Ferrer, J. C. Smith, K. Burke, The Hubbard dimer: A density functional case study of a many-body problem, Journal of Physics: Condensed Matter 27 (39) (2015) 393001. doi:10.1088/0953-8984/27/39/393001 .[47] J. I. Fuks, N. T. Maitra, Charge transfer in time-dependent density-functional theory: Insights from the asym-metric Hubbard dimer, Physical Review A 89 (6) (2014) 062502. doi:10.1103/PhysRevA.89.062502 .[48] E. Kamil, R. Schade, T. Pruschke, P. E. Bl¨ochl, Reduced density-matrix functionals applied to the Hubbarddimer, Physical Review B 93 (8) (2016) 085141. doi:10.1103/PhysRevB.93.085141 .[49] S. Balasubramanian, J. K. Freericks, Exact Time Evolution of the Asymmetric Hubbard Dimer, Journal ofSuperconductivity and Novel Magnetism 30 (1) (2017) 97–102. doi:10.1007/s10948-016-3811-9 .[50] J. ˇCis´arov´a, J. Streˇcka, Exact solution of a coupled spin-electron linear chain composed of localized Ising spinsand mobile electrons, Physics Letters A 378 (38-39) (2014) 2801 – 2807. doi:http://dx.doi.org/10.1016/j.physleta.2014.07.049 .[51] M. Nalbandyan, H. Lazaryan, O. Rojas, S. Martins de Souza, N. Ananikian, Magnetic, Thermal, and Entangle-ment Properties of a Distorted Ising-Hubbard Diamond Chain, Journal of the Physical Society of Japan 83 (7)(2014) 074001. doi:10.7566/JPSJ.83.074001 .[52] L. G´alisov´a, J. Streˇcka, Magnetic Gr¨uneisen parameter and magnetocaloric properties of a coupled spin–electrondouble-tetrahedral chain, Physics Letters A 379 (39) (2015) 2474–2478. doi:10.1016/j.physleta.2015.07.007 .[53] L. G´alisov´a, J. Streˇcka, Vigorous thermal excitations in a double-tetrahedral chain of localized Ising spinsand mobile electrons mimic a temperature-driven first-order phase transition, Physical Review E 91 (2) (2015)022134. doi:10.1103/PhysRevE.91.022134 .[54] H. ˇCenˇcarikov´a, J. Streˇcka, M. L. Lyra, Reentrant phase transitions of a coupled spin-electron model on doublydecorated planar lattices with two or three consecutive critical points, Journal of Magnetism and MagneticMaterials 401 (2016) 1106 – 1122. doi:http://dx.doi.org/10.1016/j.jmmm.2015.11.018 .[55] N. Ananikian, ˇC. Burd´ık, L. Ananikyan, H. Poghosyan, Magnetization Plateaus and Thermal Entanglement ofSpin Systems, Journal of Physics: Conference Series 804 (2017) 012002. doi:10.1088/1742-6596/804/1/ , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 .[56] H. ˇCenˇcarikov´a, J. Streˇcka, Enhanced magnetoelectric e ff ect of the exactly solved spin-electron model on adoubly decorated square lattice in the vicinity of a continuous phase transition, Physical Review E 98 (6) (2018)062129. doi:10.1103/PhysRevE.98.062129 .[57] H. S. Sousa, M. S. S. Pereira, I. N. de Oliveira, J. Streˇcka, M. L. Lyra, Phase diagram and re-entrant fermionicentanglement in a hybrid Ising-Hubbard ladder, Physical Review E 97 (5) (2018) 052115. doi:10.1103/PhysRevE.97.052115 .[58] P. M. Morse, Diatomic Molecules According to the Wave Mechanics. II. Vibrational Levels, Physical Review34 (1) (1929) 57–64. doi:10.1103/PhysRev.34.57 .[59] L. A. Girifalco, V. G. Weizer, Application of the Morse Potential Function to Cubic Metals, Physical Review114 (3) (1959) 687–690. doi:10.1103/PhysRev.114.687 .[60] E. Gr¨uneisen, Theorie des festen Zustandes einatomiger Elemente, Annalen der Physik 344 (12) (1912) 257–306. doi:10.1002/andp.19123441202 .[61] T. Balcerzak, K. Szałowski, M. Jaˇsˇcur, A simple thermodynamic description of the combined Einstein andelastic models, Journal of Physics: Condensed Matter 22 (42) (2010) 425401. doi:10.1088/0953-8984/22/42/425401 .[62] A. M. Krivtsov, V. A. Kuz’kin, Derivation of equations of state for ideal crystals of simple structure, Mechanicsof Solids 46 (3) (2011) 387. doi:10.3103/S002565441103006X .[63] T. Balcerzak, K. Szałowski, M. Jaˇsˇcur, A self-consistent thermodynamic model of metallic systems. Applicationfor the description of gold, Journal of Applied Physics 116 (4) (2014) 043508. doi:10.1063/1.4891251 .[64] T. Balcerzak, K. Szałowski, M. Jaˇsˇcur, Self-consistent model of a solid for the description of lattice and magneticproperties, Journal of Magnetism and Magnetic Materials 426 (2017) 310–319. doi:10.1016/j.jmmm.2016.11.107 .[65] T. Balcerzak, K. Szałowski, M. Jaˇsˇcur, Thermodynamic model of a solid with RKKY interaction and magnetoe-lastic coupling, Journal of Magnetism and Magnetic Materials 452 (2018) 360–372. doi:10.1016/j.jmmm.2017.12.088 .[66] K. Szałowski, T. Balcerzak, M. Jaˇsˇcur, Thermodynamics of a model solid with magnetoelastic coupling, Journalof Magnetism and Magnetic Materials 445 (2018) 110–118. doi:10.1016/j.jmmm.2017.08.073 .[67] L. Goodwin, A. J. Skinner, D. G. Pettifor, Generating Transferable Tight-Binding Parameters: Application toSilicon, Europhys. Letters 9 (1989) 701–706. doi:10.1209/0295-5075/9/7/0151 .[68] I. G. Lang, Yu. A. Firsov, Kinetic Theory of Semiconductors with Low Mobility, Zh. Eksp. Teor. Fiz. 43 (1962)1843 (Sov. Phys. JETP 16 (1963) 1301). .[69] W. Stephan, M. Capone, M. Grilli, C. Castellani, Influence of electron-phonon interaction on superexchange,Phys. Lett. A 227 (1997) 120–126. doi:0.1016/S0375-9601(96)00950-4 . , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 [70] V. M. Stojanovi´c, P. A. Bobbert, M. A. J. Michels, Nonlocal electron-phonon coupling: Consequences for thenature of polaron states, Phys. Rev. B 69 (2004) 144302. doi:10.1103/PhysRevB.69.144302 .[71] J. Koch, Quantum transport through single-molecule devices. PhD thesis, Freie Universit¨at Berlin (2006). https://refubium.fu-berlin.de/handle/fub188/13330 .[72] M. Hohenadler, W. von der Linden, Lang-Firsov Approaches to Polaron Physics: From Variational Meth-ods to Unbiased Quantum Monte Carlo Simulations, In: Alexandrov A.S. (eds) Polarons in AdvancedMaterials. Springer Series in Materials Science, vol 103. Springer, Dordrecht (2007). doi:10.1007/978-1-4020-6348-0_11 . , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 (- U / ) / t k B T / t U / t : 0 1 2 5 8 10 x = 1 E |e| r / t = 0H / t = 0 Figure 1: The normalized deviation of the chemical potential from the value of U / U / t inthe absence of the external fields. (- U / ) / t H / t U / t : 0 1 2 5 8 10x = 1 E |e| r / t = 0k B T / t = 0.5 Figure 2: The normalized deviation of the chemical potential from the value of U / U / t ,in the absence of the electric field and for moderate temperature k B T / t = . , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 (- U / ) / t E |e| r / t U / t : 0 1 2 5 8 10x = 1 H / t = 0k B T / t = 0.5 Figure 3: The normalized deviation of the chemical potential from the value of U / U / t , inthe absence of the magnetic field and for moderate temperature k B T / t = . H / t U / t : 0 1 2 4
5 6 8 10x = 1 E |e| r / t = 0 k B T / t = 10 -3 Figure 4: The relative deformation of dimer length as a function of the normalized magnetic field, for a few selectedvalues of normalized on-site Hubbard energy U / t , in the absence of the electric field and for low temperature k B T / t = − . , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 T , H t / ( | e | r ) H / t U / t : 5 6 7 8 9 10 x = 1 E |e| r / t = 0.5k B T / t = 10 -3 Figure 5: The normalized electrostriction coe ffi cient as a function of the normalized magnetic field, for a few selectedvalues of normalized on-site Hubbard energy U / t , for the normalized electric field of E | e | d / t = . k B T / t = − . E |e| r / t U / t : 0 1 2 3 4 5x = 1 H / t = 4.8k B T / t = 10 -3 Figure 6: The relative deformation of dimer length as a function of the normalized electric field, for a few selected val-ues of normalized on-site Hubbard energy U / t , for the normalized magnetic field H / t = . k B T / t = − . , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 H / t = 4.8 U / t = T , E t E |e| r / t k B T / t : 0.05 0.1U / t = 1: U / t = 2: U / t = 3: U / t = 4: Figure 7: The normalized magnetostriction coe ffi cient as a function of the normalized electric field, for a few selectedvalues of normalized on-site Hubbard energy U / t , for the normalized magnetic field of H / t = . k B T / t = .
05 and 0 . x = 1 E |e| r / t = 0H / t = 0 H , E t / k B k B T / t U / t : 0 1 2 3 6 10 Figure 8: The normalized thermal expansion coe ffi cient as a function of the normalized temperature, for a few selectedvalues of normalized on-site Hubbard energy U / t , in the absence of the external fields. , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 Figure 9: Contour plots of (a) relative deformation of dimer length; (b) normalized magnetostriction coe ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t = k B T / t = . ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t = k B T / t = .0.
05 and 0 . x = 1 E |e| r / t = 0H / t = 0 H , E t / k B k B T / t U / t : 0 1 2 3 6 10 Figure 8: The normalized thermal expansion coe ffi cient as a function of the normalized temperature, for a few selectedvalues of normalized on-site Hubbard energy U / t , in the absence of the external fields. , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 Figure 9: Contour plots of (a) relative deformation of dimer length; (b) normalized magnetostriction coe ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t = k B T / t = . ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t = k B T / t = .0. , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 Figure 11: Contour plots of (a) relative deformation of dimer length; (b) normalized magnetostriction coe ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t = k B T / t = . ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t = k B T / t = .0.
05 and 0 . x = 1 E |e| r / t = 0H / t = 0 H , E t / k B k B T / t U / t : 0 1 2 3 6 10 Figure 8: The normalized thermal expansion coe ffi cient as a function of the normalized temperature, for a few selectedvalues of normalized on-site Hubbard energy U / t , in the absence of the external fields. , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 Figure 9: Contour plots of (a) relative deformation of dimer length; (b) normalized magnetostriction coe ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t = k B T / t = . ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t = k B T / t = .0. , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 Figure 11: Contour plots of (a) relative deformation of dimer length; (b) normalized magnetostriction coe ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t = k B T / t = . ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t = k B T / t = .0. , 121740 (2019),DOI:10.1016 / j.physa.2019.121740 Figure 13: Contour plots of (a) relative deformation of dimer length; (b) normalized magnetostriction coe ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t =
10 and k B T / t = . ffi cient; (c)normalized electrostriction coe ffi cient - as a function of normalized electric and magnetic field, for U / t =
10 and k B T / t = .0.