Hubbard pair cluster with elastic interactions. Studies of thermal expansion, magnetostriction and electrostriction

TL;DR

This paper investigates the effects of elastic interactions on a Hubbard pair cluster, analyzing thermal expansion, magnetostriction, and electrostriction across various temperatures and external fields within an extended Hubbard model.

Contribution

It introduces a comprehensive study of elastic interactions in a Hubbard dimer, deriving an equation of state and analyzing deformation effects under different conditions.

Findings

Thermal expansion varies with temperature and external fields.

Magnetostriction and electrostriction effects are characterized across parameter ranges.

Deformation responses depend on the interplay of elastic and electronic interactions.

Abstract

The pair cluster (dimer) is studied within the framework of the extended Hubbard model and the grand canonical ensemble. The elastic interatomic interactions and thermal vibrational energy of the atoms are taken into account. The total grand potential is constructed, from which the equation of state is derived. In equilibrium state, the deformation of cluster size, as well as its derivatives, are studied as a function of the temperature and the external magnetic and electric fields. In particular, the thermal expansion, magnetostriction and electrostriction effects are examined for arbitrary temperature, in a wide range of Hamiltonian parameters.

Figures

Equations

$${\mathcal{H}}_{a,b}$$

$$n_{\gamma ,\sigma }=c_{\gamma ,\sigma }^{+}{c}_{\gamma ,\sigma }.$$

$$n_{\gamma }=n_{\gamma ,\uparrow }+n_{\gamma ,\downarrow }.$$

$$S_{\gamma }^z=\left(n_{\gamma ,\uparrow }-n_{\gamma ,\downarrow }\right)/2.$$

$$d=d_0\left(1+\epsilon \right),$$

$$V=V_0\left(1+\epsilon \right),$$

$$W=\frac{W_0}{1+\epsilon },$$

$$t=t_0{\left(\frac{d}{d_0}\right)}^{-n}=t_0{\left(1+\epsilon \right)}^{-n},$$

$${\Omega }_{a,b}=-k_{\mathrm{B}}T\ln \{{\mathrm{Tr}}_{a,b}\exp [-\beta \left({\mathcal{H}}_{a,b}-\mu \left(n_a+n_b\right)\right)]\}.$$

$$2x=⟨n_a⟩+⟨n_b⟩=-{\left(\frac{\partial {\Omega }_{a,b}}{\partial \mu }\right)}_{T,H,E}.$$

$$V_M\left(r\right)=D{\left[1-e^{-\alpha \left(r-r_0\right)/r_0}\right]}^2.$$

$$V_M\left(\epsilon \right)=D{\left[1-e^{-\alpha \left(\frac{d_0}{r_0}\left(1+\epsilon \right)-1\right)}\right]}^2-D{\left[1-e^{-\alpha \left(\frac{d_0}{r_0}-1\right)}\right]}^2.$$

$$E_{\omega }=k_{\mathrm{B}}T\ln \left[2\sinh \left(\frac{1}{2}\beta \hslash \omega \right)\right].$$

$$\omega =\frac{{\omega }_0}{{\left(1+\epsilon \right)}^{\Gamma }},$$

$$\Gamma =-\frac{d}{\omega }\left(\frac{\partial \omega }{\partial d}\right).$$

$$E_{\omega }=k_{\mathrm{B}}T\ln \left[2\sinh \left(\frac{t_0}{k_{\mathrm{B}}T}\frac{\hslash {\omega }_0}{t_0}\frac{1}{2{\left(1+\epsilon \right)}^{\Gamma }}\right)\right],$$

$$\Omega ={\Omega }_{a,b}+V_M\left(\epsilon \right)+E_{\omega }.$$

$${\left(\frac{\partial \Omega }{\partial \epsilon }\right)}_{T,H,E}=0.$$

$$F_{a,b}=-\frac{1}{d_0}\left(\frac{\partial {\Omega }_{a,b}}{\partial \epsilon }\right).$$

$$-d_0{F}_{a,b}=\frac{\partial }{\partial \epsilon }⟨{\mathcal{H}}_{a,b}⟩-2x\frac{\partial \mu }{\partial \epsilon },$$

$$\mu =\frac{W_0}{1+\epsilon }+\frac{U}{2},$$

$$F_{a,b}$$

$$F_{\epsilon }=-\frac{1}{d_0}\left(\frac{\partial V_M\left(\epsilon \right)}{\partial \epsilon }\right),$$

$$F_{\epsilon }=-\frac{2D\alpha }{r_0}\left[1-e^{-\alpha \left(\frac{d_0}{r_0}\left(1+\epsilon \right)-1\right)}\right]e^{-\alpha \left(\frac{d_0}{r_0}\left(1+\epsilon \right)-1\right)}.$$

$$F_{\omega }=-\frac{1}{d_0}\left(\frac{\partial E_{\omega }}{\partial \epsilon }\right),$$

$$F_{\omega }=\frac{\Gamma }{2{\left(1+\epsilon \right)}^{\Gamma +1}}\frac{\hslash {\omega }_0}{d_0}{\tanh }^{-1}\left(\frac{t_0}{k_{\mathrm{B}}T}\frac{\hslash {\omega }_0}{t_0}\frac{1}{2{\left(1+\epsilon \right)}^{\Gamma }}\right).$$

$$F_{a,b}+F_{\epsilon }+F_{\omega }=0.$$

Full text

††journal:

In Memory of Professor Leszek Wojtczak

1 Introduction

The Hubbard model [1, 2, 3, 4] plays an important role in contemporary solid state physics. Since its formulation, numerous applications of the model have been developed and the model itself has been extensively investigated as a fundamental, prototypic model for description of correlated fermions [5, 6, 7, 8]. Many of the studies were concentrated on the typical three-dimensional, two-dimensional or one-dimensional Hubbard model (e.g. [9, 10, 11]). However, with an increase of 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 of atoms have been performed. Some of them involve the studies of small clusters with either exact or close to exact methods [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29]. Among such systems, a two-atomic cluster (dimer) plays a special role, being the smallest nanosystem where the Hubbard model can be adopted, and the exact analytical solutions can be found. As a consequence, numerous aspects of physics of Hubbard dimers focused considerable attention in the literature [30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49]. These facts serve as a sound motivation for further comprehensive study of the two-atomic Hubbard system aimed at its full and exact thermodynamic description. In line with this trend, recently the Hubbard pair cluster embedded simultaneously in the external magnetic and electric fields has been studied by the exact analytical methods in the framework of the grand canonical ensemble. In this approach the system is open, so that it is able to exchange the electrons with its environment [41, 42, 43]. The statistical and thermodynamic properties, both magnetic and electric ones, have been extensively studied there.

Despite these investigations, there is still room for extending the model of interest. Supplementing purely electronic models by including further degrees of freedom can cause new kinds of behaviour to emerge and new control parameters can be gained this way. Purely electronic models can be supplemented, for example, with subsystems composed of localized spins (e.g. [50, 51, 52, 53, 54, 55, 56, 57]), which are still exactly solvable. Immersing of such system in the external fields can give access to unique properties (like, for example, manifestations of the magnetoelectric effect [56]). Enriching the purely electronic models is also a step towards more complete characterization of the real physical systems. A good example can be inclusion of the Hamiltonian terms, which are responsible for coupling of the electronic degrees of freedom with the lattice. Such a study, aimed at the description of magnetomechanical and electromechanical properties of the Hubbard dimer, has not 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 description of the Hubbard dimer. Such interactions are especially important since they are responsible for stability of the dimer structure and allow to account for the energy of thermal vibrations. It can be expected that after combining with the Hubbard Hamiltonian, the elastic interactions will lead to new phenomena, connecting the magnetic, electric and mechanical properties.

The elastic interaction between the atoms is assumed here in the form of the Morse potential [58, 59]. The thermal vibrational energy is obtained in the quasi-harmonic approximation, where the anharmonicity is taken into account by the Grüneisen parameter [60, 61]. In turn, the Hubbard pair Hamiltonian is taken in its extended form, with the Coulomb interaction between the electrons on different 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 manner. In particular, the deformation of interatomic distance (dimer length), magnetostriction and electrostriction coefficients are calculated for various temperatures, and the effect of the external magnetic and electric fields on the mentioned quantities is also investigated. In addition, the chemical potential is studied in the presence of elastic interactions, showing a new behaviour in comparison with the conventional Hubbard model. A special attention is paid to the low-temperature region, 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 and the basic formulae, important for numerical calculations, are derived. In the successive Section the selected results of calculations are illustrated in figures and discussed. The last Section is devoted to 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 Hubbard dimer including the elastic and vibrational properties is presented and the solution of the model is described.

2.1 Grand potential for the Hubbard dimer

Extended Hubbard Hamiltonian for the pair of atoms $(a,b)$ embedded in the external fields is of the form:

[TABLE]

where $t>0$ is the hopping integral, $U\geq 0$ is on-site Coulomb repulsion and $W$ stands for the Coulomb interaction between electrons localized on $a$ and $b$ atoms. In 1, $H=-g\mu_{\rm B}H^{z}$ denotes the magnitude of an external uniform magnetic field $H^{z}$, whereas the electric field $E$ is determined by the electrostatic potentials $V_{a}=-V_{b}=V$ applied to both atoms. Thus, $E=2V/\left(|e|d\right)$, where $d$ is the interatomic distance and $e$ is the electron charge. The quantity $W$ in 1 is assumed in the classical form: $W=1/(4\pi\epsilon_{0})\,e^{2}/d$.

The creation ($c_{\gamma,\sigma}^{+}$) and annihilation ($c_{\gamma,\sigma}$) operators for site $\gamma=a,b$ and spin state $\sigma$ can be used to construct the occupation number operators $n_{\gamma,\sigma}$ for the electrons with given spin:

[TABLE]

The total occupation number operators at the site $\gamma=a,b$ are then given by:

[TABLE]

Moreover, the spin operators $S_{\gamma}^{z}$ are defined as follows:

[TABLE]

Since we intend to study the elastic properties of the model, the interatomic distance $d$ can be presented as follows:

[TABLE]

where $d_{0}$ is the equilibrium distance defined in the absence of the external fields ($H=0,\,E=0$) and for zero temperature ($T=0$). In 5, $\varepsilon$ defines the temperature and fields-induced deformation. As it can be noted, some of the coefficients in Hamiltonian 1, like $V$, $W$ and $t$, are interatomic distance-dependent. For instance, for the electric potential we have:

[TABLE]

where $V_{0}=E|e|d_{0}/2$. In turn, for the interatomic Coulomb interaction we can write:

[TABLE]

where $W_{0}=\frac{e^{2}}{4\pi\epsilon_{0}d_{0}}$. Moreover, for the hopping integral we assume the power-law distance dependence, namely:

[TABLE]

with a constant power index $n$, which can be used if the deformation $\varepsilon$ is small.

We treat the Hubbard pair as an open system, with a variable number of electrons, which can be 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 the average number of electrons in the system results from the thermodynamic equilibrium. In such small system the relative fluctuations of the number of electrons, when related to its statistical mean value, can be significant. Therefore, in this case a formalism of the grand canonical ensemble should be used, which is more appropriate here than the canonical ensemble, with a fixed number of electrons.

In the framework of the grand canonical ensemble, the Hamiltonian 1 should be extended by adding the term $-\mu\left(n_{a}+n_{b}\right)$, where $\mu$ 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 diagonalization procedure will not increase in comparison to that presented in Ref. [41]. We note that $n_{a}n_{b}$ is a diagonal 16 $\times$ 16 matrix and corresponds to 16 possible states. Its diagonal elements are: $a_{1,1}=a_{2,2}=a_{3,3}=a_{4,4}=a_{5,5}=0$, $a_{6,6}=a_{7,7}=1$, $a_{8,8}=2$, $a_{9,9}=0$, $a_{10,10}=a_{11,11}=1$, $a_{12,12}=2$, $a_{13,13}=0$, $a_{14,14}=a_{15,15}=2$, $a_{16,16}=4$. Such a matrix, multiplied by $W/t$, should be simply added to that presented in Ref. [41] (Eq. B.5 therein).

As a result of diagonalization, the grand thermodynamic potential $\Omega_{a,b}$ for the Hubbard pair can be found from the general formula:

[TABLE]

If we denote the mean number of the electrons per atom by $x$, then the following relation is satisfied:

[TABLE]

In Eq. 10, $\left\langle\cdots\right\rangle$ means the statistical average, which is carried out over all 16 states, and is dependent on the temperature and external fields. From Eq. 10 the chemical potential $\mu$ can be found, provided that $x$ is established, and then the expression 9 is complete.

2.2 Elastic (static) interaction

Irrespective of the Hubbard contribution, the interatomic static interaction of the core electrons can be taken into account in the form of the Morse potential energy [58]:

[TABLE]

In Eq. 11, $D$ is the potential depth, parameter $\alpha$ controls the potential width and its asymmetry, whereas $r_{0}$ is the equilibrium distance (potential minimum position) for this interaction at $T=0$. The static energy can be written in terms of deformation $\varepsilon$, and its value can be normalized by the requirement that $V_{M}\left(\varepsilon=0\right)=0$. In this way we obtain the elastic energy in the form of:

[TABLE]

The dimensionless parameter $d_{0}/r_{0}$ can be found from the minimum condition of the total energy, taken at $T=0$ in the absence of the external fields.

2.3 Vibrational (thermal) energy

In order to take into account the thermal vibrations, the two-atomic system can be treated as a 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, which forms a phononic bath. For the sake of simplicity we adopt the quasi-harmonic Einstein model, in which the anharmonicity is taken into account by the Grüneisen parameter $\Gamma$ [61]. The vibrational free energy of the Einstein oscillator is then given by:

[TABLE]

The frequency of oscillations, $\omega$, is dependent on the interatomic distance and can be expressed as a function of the deformation $\varepsilon$ as follows:

[TABLE]

where the Grüneisen parameter $\Gamma$ satisfies the relationship:

[TABLE]

It is known that $\Gamma$ can be determined theoretically for several interatomic potentials and different dimensionalities of the system. For instance, for the Morse potential and 1D system $\Gamma=\left(3/2\right)\alpha$ [62]. This convenient relationship will be adopted in further calculations.

Thus, the vibrational energy can be finally presented as:

[TABLE]

where $k_{\mathrm{B}}T/t_{0}$ is the dimensionless temperature, and $\hbar\omega_{0}/t_{0}$ is the dimensionless characteristic energy parameter depending on the atomic mass.

2.4 The equation of state (EOS)

The total grand potential, $\Omega$, of the system in question can be constructed as a sum of the Hubbard grand potential 9, elastic (static) energy 12 and vibrational Helmholtz free-energy 16:

[TABLE]

This grand potential should be minimized first with respect to the deformation $\varepsilon$ by imposing the necessary equilibrium condition:

[TABLE]

Differentiation of $\Omega$ 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 under the deformation:

[TABLE]

With the help of Eq. 9 we get:

[TABLE]

where $\left\langle\mathcal{H}_{a,b}\right\rangle$ is the statistical average of the Hubbard Hamiltonian. The energy $\left\langle\mathcal{H}_{a,b}\right\rangle$ can be found exactly after the process of diagonalization. Then, the derivative $\partial\left\langle\mathcal{H}_{a,b}\right\rangle/\partial\varepsilon$ can be calculated analytically with the help of Eqs. 6–8. In turn, the derivative $\partial\mu/\partial\varepsilon$ can be, in principle, calculated numerically for the arbitrary concentration $\left(0\leq x\leq 2\right)$ of Hubbard electrons. For such purpose Eq. 10 for the chemical potential can be used. Interestingly, for $x=1$, i.e., for half-filled energy levels, the analytical formula for the chemical potential $\mu$ can be found, namely:

[TABLE]

which effectively facilitates the calculations of $\partial\mu/\partial\varepsilon$. As a result, the Hubbard force for $x=1$ can be presented in the following analytic form:

[TABLE]

The statistical averages in Eq. 22 can be calculated after performing the diagonalization of the Hamiltonian. The examples of similar calculations (in the absence of elastic interactions) have already been presented in Ref. [42].

As a next step, the elastic force, $F_{\varepsilon}$, is found from the formula:

[TABLE]

where $V_{M}\left(\varepsilon\right)$ is given by Eq. 12. Hence, we obtain:

[TABLE]

The last force, $F_{\omega}$, which is of vibrational origin, is defined as:

[TABLE]

where $E_{\omega}$ is given by Eq. 16. From Eqs. LABEL:[nosort]eq25,eq16 we obtain:

[TABLE]

It follows from the numerical calculations that the Hubbard force, $F_{a,b}$, is compressive (negative) around $\varepsilon=0$. On the other hand, the elastic force, $F_{\varepsilon}$, is compressive for $d>r_{0}$ and expansive (positive) for $d<r_{0}$, whereas the vibrational force, $F_{\omega}$, is always of expansive type. The necessary condition for thermodynamic equilibrium requires the balance of all these forces, namely:

[TABLE]

Expression 27 is the EOS for the Hubbard dimer with the elastic interactions and thermal vibrations taken into account. It enables calculation of spontaneous deformation, $\varepsilon$, of the dimer length for given values of the intensive parameters: $T$, $H$ and $E$, whereas the system is in stable equilibrium. At first, in the absence of deformation (i.e., when $\varepsilon=0$ for $T=0$, $H=0$ and $E=0$), the constant parameter $d_{0}/r_{0}$ should be found from Eq. 27, giving the equilibrium distance $d_{0}$ at the ground 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 number of electrons per atom is $x=1$. In the Hubbard Hamiltonian we select the exponent $n$ in the hopping integral (8) equal to $n=6$ in order to describe a rapid variability of the function around $d_{0}$. 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 of the hopping integrals dependence on the bond length can be given by the Goodwin, Skinner and Pettifor function [67]. The interatomic Coulomb interaction parameter $\frac{W_{0}d_{0}}{r_{0}t_{0}}=\frac{e^{2}}{4\pi\epsilon_{0}r_{0}t_{0}}$ can be estimated for the realistic values of $r_{0}\approx 1.5~{}\AA$ and $t_{0}\approx 2.5~{}{\rm eV}$ giving a figure of $\approx$ 3.83, and for further calculations a round value of $4$ is accepted. All dimensionless energies, temperature and fields magnitudes are normalized by $t_{0}$ and the interatomic equilibrium distance $d_{0}$ is related to the $r_{0}$ constant. In this convention, we chose the value $D/t_{0}=4$ for the Morse potential depth, while the asymmetry parameter is assumed with its typical value $\alpha=4$ [59]. For the thermal vibrations we assume the characteristic energy parameter $\hbar\omega_{0}/t_{0}=0.5$, whereas the Grüneisen parameter is $\Gamma=\left(3/2\right)\alpha$ [62]. Such a set of parameter values, although not describing any specific material, seems physically acceptable and can be used to demonstrate how the formalism works when the numerical calculations are performed. We are aware that the choice of these values can be discussed and modified according to the physical situation. Unfortunately, because of the limited scope of the paper we are not able to present the influence of all these parameters on the numerical results. However, it can be mentioned that, for instance, some change of the asymmetry parameter $\alpha$, or the potential depth $D$, may have remarkable influence on the thermodynamic properties, as it has been demonstrated in Ref. [65].

For the above set of parameters and for $U/t_{0}=0$, from Eq. 27 we obtain the equilibrium distance $d_{0}/r_{0}=0.922085$ in the ground state (i.e., for $T=0$, $H=0$ and $E=0$, whereas deformation is $\varepsilon=0$). Analogously, for $U/t_{0}=1,\,2,\,5$ and $10$, the corresponding equilibrium distances are: $d_{0}/r_{0}=0.918756$, 0.915598, 0.915319 and 0.941767, respectively. At the same time, in the ground state the chemical potentials, $\mu/t_{0}$, for $U/t_{0}=0,\,1,\,2,\,5$ and $10$ are as follows: $\mu/t_{0}=4.337994$, 4.853712, 5.368729, 6.870059 and 9.247335, respectively.

We see from the above results that the equilibrium interatomic distance $d_{0}$ is not a monotonous function of the parameter $U$. Moreover, the chemical potential $\mu$ is not equal to $U/2$, as it happens in the case when the interatomic Coulomb repulsion, as well as the elastic interactions, are neglected [5, 41]. The values of $\mu$ obtained here are in agreement with Eq. 21. The formula 21 also predicts that out of the ground state, i.e., for $T>0$, and in the possible presence of the external fields, the chemical potential is no longer constant, but it should be a function of $T$, $H$ and $E$, via its dependency on the deformation $\varepsilon$. This new behaviour is illustrated in Fig. 1-3. For instance, in 1 the difference $(\mu-U/2)/t_{0}$ is shown as a function of the dimensionless temperature $k_{\rm B}T/t_{0}$ for several values of $U/t_{0}$, when the external fields are absent. In general, the functions decrease with the increasing temperature and, at the same time, the dependency of $(\mu-U/2)$ on $U$ is not monotonous. In 2 also the decreasing tendency of $(\mu-U/2)/t_{0}$ as a function of dimensionless magnetic field $H/t_{0}$ is shown for the same values of $U/t_{0}$. In this case the electric field is equal to $E=0$, and the moderate temperature, $k_{\rm B}T/t_{0}=0.5$, is assumed. In 3 the chemical potential is plotted for $k_{\rm B}T/t_{0}=0.5$ and $H=0$, as a function of the dimensionless electric field $E|e|r_{0}/t_{0}$. Apart from the decreasing tendency of $(\mu-U/2)/t_{0}$ vs. $E$, which is visible for small $U$, the increasing tendency of the chemical potential for $U/t_{0}=8$ and $U/t_{0}=10$ is worth noticing.

4 presents relative deformation $\varepsilon$ of the dimer size vs. magnetic field $H/t_{0}$ at very low temperature, $k_{\rm B}T/t_{0}=0.001$, whereas $E=0$. Several values of $U/t_{0}$ are examined. As we know from our previous studies [41, 42, 43], in the low-temperature region the quantum phase transitions should occur at the magnetic (or electric) critical fields. Such transitions are of the first order, as regards the behaviour of total dimer magnetization, which results from the change from the singlet state, where $\left\langle S^{z}_{a}+S^{z}_{b}\right\rangle=0$, to a triplet one, where $\left\langle S^{z}_{a}+S^{z}_{b}\right\rangle=1$. The rapid change of magnetization is accompanied by instantaneous changes of the remaining thermodynamic quantities. In 4 the changes of deformation $\varepsilon$ are illustrated at the transition from the singlet state, where $\varepsilon=0$, to the triplet (ferromagnetically saturated) state. It should be mentioned that due to the presence of interatomic Coulomb interaction $W$, the critical fields $H_{c}$ found here have different values in comparison with those from our previous papers [41, 42, 43]. Unfortunately, there is no analytical formula for these critical fields, and they can be calculated numerically only by equating the total grand potentials of neighbouring phases.

In 5 the electrostriction coefficient, $\nu_{T,H}=\frac{1}{d}\left(\frac{\partial d}{\partial E}\right)_{T,H}=\frac{1}{1+\varepsilon}\left(\frac{\partial\varepsilon}{\partial E}\right)_{T,H}$, is presented in dimensionless units, $\frac{\nu_{T,H}t_{0}}{|e|r_{0}}$, as a function of the magnetic field $H/t_{0}$. Several values of the parameter $U$ are examined at the temperature equal to $k_{\rm B}T/t_{0}=0.001$, whereas the electric field is constant and equal to $E|e|r_{0}/t_{0}=0.5$. As it is seen, the characteristic steps corresponding to the quantum phase transitions occur for the electrostriction coefficient. The critical magnetic fields $H_{c}$, where the jumps occur, decrease with an increase in the Coulombic repulsion energy $U$, which effect is analogous to that noted in the previous figure. Moreover, for higher $U$ values, the existence of negative electrostriction is demonstrated. It can also be noted that the electrostriction coefficient takes zero value in the triplet (ferromagnetically saturated) state, because then the electric field is not able to change the charge distribution [43], whereas the electric potential is compensated by the chemical potential. As a result, the total energy is not influenced by the electric field and length of the dimer remains unchanged.

In the next figure, 6, the deformation $\varepsilon$ is presented at the temperature $k_{\rm B}T/t_{0}=0.001$ and the magnetic field $H/t_{0}=4.8$ as a function of dimensionless electric field $E|e|r_{0}/t_{0}$. Again, the characteristic jumps of $\varepsilon$ are shown for various curves corresponding to different 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 $\varepsilon$ is constant vs. $E$, which results in null electrostriction effect. This fact supports the discussion presented in the preceding paragraph and is in agreement with the general behaviour of the critical fields for different $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 speaking about the critical fields $E_{c}$ (or $H_{c}$) for given $U$-parameter and the fields $H$ (or $E$), we can equally well 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_{\rm B}T/t_{0}=0.05$ and $k_{\rm B}T/t_{0}=0.1$ for the same magnetic field as in previous figure, $H/t_{0}=4.8$. The magnetostriction coefficient is defined as $\lambda_{T,E}=\frac{1}{d}\left(\frac{\partial d}{\partial H}\right)_{T,E}=\frac{1}{1+\varepsilon}\left(\frac{\partial\varepsilon}{\partial H}\right)_{T,E}$ and is presented in dimensionless units, $\lambda_{T,E}t_{0}$, vs. dimensionless electric field $E|e|r_{0}/t_{0}$. It is seen that the magnetostriction coefficient in lower temperature shows the sharp peaks, whereas in higher temperature the maxima are lower and much more diffused. These maxima appear approximately at the critical electric fields $E_{c}$ and are placed at different positions for the curves plotted for various parameters $U$. In particular, when $U$ increases, the critical fields are shifted towards higher values of $T$. For $T\to 0$ the magnetostriction coefficient is zero everywhere aside from the singularities, because in the low-temperature region the magnetic field is not able to change the magnetization, both in singlet and in triplet state, except for the phase transition point [42]. As a result, the magnetic energy is constant aside from the singularities and changes of the magnetic field do not influence the dimer size. This conclusion is also in agreement with 4 where, apart from the phase transition points, the plots of $\varepsilon$ vs. $H$ are flat, which corresponds to the null magnetostriction effect. When the temperature increases, the quantum phase transitions become diffused and, as a result, the magnetostriction coefficient presents broadened maxima with reduced height. The positions of the peaks are quite stable and are only slightly shifted towards larger electric fields. At the same time, the height of the peaks strongly depends on the Hubbard parameter $U$, showing a decrease of maximal magnetostriction coefficient when $U$ increases.

The coefficient of the linear thermal expansion is presented in 8 as a function of the dimensionless temperature $k_{\rm B}T/t_{0}$. This coefficient is defined as $\alpha_{H,E}=\frac{1}{1+\varepsilon}\left(\frac{\partial\varepsilon}{\partial T}\right)_{H,E}$, where the derivative over $T$ is taken at constant fields $H$ and $E$. In 8 the dimensionless quantity $\alpha_{H,E}\,t_{0}/k_{\rm B}$ is plotted, whereas the fields $H$ and $E$ are absent. The most interesting effect illustrated in this figure is the maximum of the thermal expansion and its evolution when the Coulombic repulsion energy $U$ increases. This maximum is initially small for the curve corresponding to $U=0$ and gradually disappears for the curves with $U/t_{0}=1$ and $U/t_{0}=2$. Then, for $U/t_{0}=3$ it reappears at higher temperature and becomes very pronounced for the curves corresponding to $U/t_{0}=6$ and $U/t_{0}=10$. At the same time, the maximum is shifted towards lower temperatures again. The occurrence of this maximum is not connected with the phase transition, since the system still remains in the singlet state. It rather results from a complicated interplay and competition between different contributions to the total 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$-potential brings on the charge separation, and as a result the occupation correlation on different atoms, $\left\langle n_{a}n_{b}\right\rangle$, increases. For instance, in 8  for the fixed temperature $k_{\rm B}T/t_{0}=0.5$, and $U/t_{0}=2$, 6 and 10, the occupation correlations amount to: $\left\langle n_{a}n_{b}\right\rangle=0.241$, 0.822 and 0.9751, respectively. Then, in order to compensate the increase in interatomic Coulomb energy, the $W$ coefficient (Eq. 7) must decrease, which means that deformation $\varepsilon$ takes higher values there. One can also note that, if $d_{0}/r_{0}<1$, some increase in $\varepsilon$ helps to lower the Morse potential energy. The resulting deformation is obtained as the solution of EOS (Eq. 27). For the above example ($k_{\rm B}T/t_{0}=0.5$ and $U/t_{0}=2$, 6 and 10), the corresponding deformations are: $\varepsilon=0.0167$, 0.0358 and 0.0400, respectively. On the other hand, for given $U$, the deformation is a non-linearly increasing function of temperature. In this case the increase in $\varepsilon$ is also connected with the energy compensation mechanism, and it is mainly stimulated by the fact that in the singlet (paramagnetic) phase the increase in temperature results in a nonlinear increase of the mean hopping energy (see [42] and Fig. 5 therein). Then, the derivative $\partial\varepsilon/\partial T$ can reveal a non-monotonous behaviour vs. temperature and the maximum of $\alpha_{H,E}$ appears as a consequence. The non-linear changes of $\varepsilon$ are more pronounced in the intermediate range of temperatures and for larger $U$-parameters, where the competition between various pressure components in EOS is more dynamic. In our opinion, the revealed maximum of $\alpha_{H,E}$ is, to some extent, analogous to the paramagnetic maximum of the specific heat, observed in localized spin systems and known as the Schottky anomaly. In that effect the magnetic energy (or entropy) also shows a non-linear increase vs. temperature.

It can also be mentioned that for $T\to 0$ the thermal expansion coefficient tends to zero, which is a proper thermodynamic behaviour. On the other hand, in the high temperature region, the anharmonic vibrations play a dominant role and, as a result, an increase in $\alpha_{H,E}$ occurs there in a more regular way.

The next figures (Fig. 9-14) are prepared in form of contour graphs, presenting in the multiple charts the deformation (in parts a), magnetostriction coefficient (in parts b) and electrostriction coefficient (in parts c). The isolines of the dimensionless quantities, $\varepsilon$, $1000\times\lambda_{T,E}t_{0}$, and $1000\times\nu_{T,H}t_{0}/(|e|r_{0})$, are plotted as a function of the external fields, in the plane ($H/t_{0},\,E|e|r_{0}/t_{0}$). Two temperatures, which are far from the ground state: $k_{\rm B}T/t_{0}=0.5$ (for Fig. 9, 11, and 13), and $k_{\rm B}T/t_{0}=1.0$ (for Fig. 10, 12, and 14) are chosen. In order to illustrate the role of Coulombic on-site repulsion $U$, three values are selected: $U/t_{0}=2$ for LABEL:fig9,fig10, $U/t_{0}=5$ for LABEL:fig11,fig12, and $U/t_{0}=10$ for LABEL:fig13,fig14.

We see that both $T$ and $U$ play a significant role in the behaviour of the quantities presented in multiple contour plots. The shapes of the isolines, defining the areas with similar values of the presented quantities, change dynamically between different plots. In particular, analysing the deformation, it can be noted that areas with the smallest values of $\varepsilon$, denoted by blue shades and situated in the left parts of the charts, where rather small magnetic fields dominate, are shifted quite irregularly along the vertical (electric field) axis for different $T$ and $U$. At the same time, the areas with highest values of $\varepsilon$, denoted by red colour and other warm shades, are expanding as the temperature increases for each value of $U$. For instance, analysing the isolines for $\varepsilon$ $=$ 0.1 in Fig. 9(a)-12(a) for given $U$, we see that the areas with $\varepsilon\geq$0.1 are increasing when $T$ increases. On the other hand, all isolines in 14(a) correspond to larger values of $\varepsilon$ than isolines in 13(a). The red areas are situated at high electric fields when $U/t_{0}=2$ and $U/t_{0}=5$, however, when $U/t_{0}=10$, the area with large $\varepsilon$ extends over all values of the electric fields down to zero, provided the magnetic field is not very small. It should be stressed that deformation $\varepsilon$ is positive for all these cases. The increasing values of $\varepsilon$ when the temperature increases indicate that also the thermal expansion coefficient should be positive.

Regarding the magnetostriction, its value is mostly positive, however, some areas with the negative magnetostriction can be found, namely for high magnetic and electric fields. Interestingly, the area with negative magnetostriction expands with the increasing temperature. On the other hand, the dependency on $U$ is not monotonous: although the area with negative magnetostriction is larger for $U/t_{0}=5$ than for $U/t_{0}=2$, however, for $U/t_{0}=10$ this range becomes strongly reduced and is only visible at higher temperature ($k_{\rm B}T/t_{0}=1.0$). The area where the highest magnetostriction occurs lies in the range belonging to small electric fields, provided that $U/t_{0}=2$ or $U/t_{0}=5$, however, it is shifted towards large electric fields if $U/t_{0}=10$. The magnetic fields are moderate in all cases of high magnetostriction ranges, being from the range of ($2\lesssim H/t_{0}\lesssim 4$) for $U/t_{0}=2$, and ($0.5\lesssim H/t_{0}\lesssim 2$) for higher $U/t_{0}$. The increase in temperature results also in some expansion of the area with strong positive magnetostriction. It can also be noted that isolines corresponding to zero magnetostriction present an interesting behaviour, namely, the increase in electric field causes the decrease in the magnetic field along the isoline, so that the slope of the line as a function of $H$ is negative.

As far as the electrostriction coefficient is concerned, some areas with the negative values of this quantity are also observed in all the cases. For small values of $U$ these areas lie at low electric fields. In particular, for $U/t_{0}=2$, rather high magnetic fields ($H/t_{0}\gtrsim 3$) are involved, whereas for $U/t_{0}=5$ the corresponding magnetic fields are much lower. They may extend down even until $H/t_{0}=0$, as in the case of 12. On the other hand, for $U/t_{0}=10$ the area with negative electrostriction occurs for considerably higher electric fields, for instance, in the order of $E|e|r_{0}/t_{0}\approx 4$, whereas the magnetic fields are still very small. The increase in temperature, from $k_{\rm B}T/t_{0}=0.5$ to $k_{\rm B}T/t_{0}=1.0$, makes the area of negative electrostriction larger, but, at the same time, its amplitude decreases approximately by an order of magnitude. On the other hand, the areas corresponding to the highest positive electrostriction are mainly situated at very large electric fields. An exception is 9 where, for $U/t_{0}=2.0$ and $k_{\rm B}T/t_{0}=0.5$, the region of positive electrostriction coefficient extends down to very small electric fields.

In general, the isolines separating the areas of positive and negative electrostriction are more densely packed for $k_{\rm B}T/t_{0}=0.5$ than for $k_{\rm B}T/t_{0}=1.0$. The oblique arrangement of these lines at $k_{\rm B}T/t_{0}=0.5$ may lead to considerable dynamics of the electrostriction coefficient when crossing them horizontally by an increase in the magnetic field. The isolines corresponding to zero electrostriction effect 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 LABEL:fig11,fig12 for $U/t_{0}=5$). This contrasts with the behaviour of the magnetostriction, where the slopes of the null magnetostriction coefficient line were always negative.

4 Summary and conclusion

In the paper the total grand potential has been constructed for the dimer system, in which the magnetic, elastic and thermal vibration energies play important roles. The similar attempts at the total Gibbs energy construction, i.e., consisting in summation of various energy contributions, have been undertaken, so far, for the bulk materials [61, 63, 64, 65, 66]. According to our knowledge, for the 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, and application of canonical transformation of Lang–Firsov type [68, 69, 70, 71, 72]. Then, the procedure is followed by using the 2nd order perturbation calculus, or, equivalently, by the Schrieffer–Wolff transformation [69], which leads to the renormalization of Hamiltonian parameters. Such a method is much more complicated and, in our opinion, is therefore less convenient for the complete studies of various thermodynamic properties. Moreover, the Lang–Firsov type transformation has been proposed for the Hamiltonians describing the coupling of electrons with the harmonic oscillators [70], whereas in our approach the anharmonic vibrations are taken into account via Grüneisen parameter $\Gamma$. It is known that the anharmonicity plays an important role, for instance, in the studies of thermal expansion effect.

As a result of the present approach, the EOS (Eq. 27) has been derived for the system in question. This equation yields stable solutions for the interatomic distance deformation, $\varepsilon$, corresponding to the thermodynamic equilibrium state obtained at arbitrary temperature $T$ and the external fields $H$ and $E$. Investigation of this deformation leads to the determination of its derivative quantities such as the thermal expansion or magneto- and electrostriction coefficients.

It has been shown in the paper that for the half-filling conditions ($x=1$) and for the elastic interactions included, the chemical potential $\mu$ is not a constant but it depends, via the quantity $\varepsilon$, on the temperature and external fields (Fig. 1-3). For $x=1$ the analytic formula describing $\mu$ has been presented (Eq. 21). It should be stressed that the method used here can be equally well applied to arbitrary average concentration of the electrons, i.e., when $0\leq\;x\leq 2$.

In the low-temperature region, where the quantum phase transitions are induced by the external fields, the discontinuous changes of deformation $\varepsilon$ have been found at the phase transition points. These discontinuous changes are also accompanied by discontinuities of such quantities as the magnetostriction and electrostriction coefficients (Fig. 4-7). It is known that with an increase in temperature the quantum phase transitions become diffused, so then, the most spectacular effects of discontinuity occur at $T\to 0$.

For the asymmetric interatomic potential, an increase in temperature results in the linear expansivity effect. In our opinion, an interesting finding concerns the broad maxima of the thermal expansion coefficient, found at some temperatures and illustrated in 8. It has been shown that the positions of maxmima are strongly dependent on the values of Hubbard $U$ parameter.

In the high-temperature region, where the discontinuities of thermodynamic quantities do not occur, the contour charts yield a comprehensive description of the interatomic distance deformation and its changes vs. external fields. The contour charts (Fig. 9-14) have been used to illustrate simultaneously $\varepsilon$, as well as the magnetostriction and electrostriction coefficients vs. $H$ and $E$, for some selection of temperatures $T$ and parameters $U$. In particular, a distinctive result is an identification of the ($H,E$)-regions where the negative values of these coefficients do appear. It has been found that both temperature $T$ and $U$-parameter have important influence on the localization and area of these regions, as well as on the values of $\varepsilon$, $\lambda_{T,E}$, and $\nu_{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. It should be stressed that diagonalization of the Hubbard Hamiltonian has been performed exactly by means of the analytical methods [41]. It can be expected that selection of other interatomic potentials would change our results quantitatively; however, the case which we consider in the present paper captures the essence of the phenomena emerging due to coupling between purely electronic and elastic properties.

It should also be pointed out that having constructed the total grand potential, all thermodynamic and statistical properties of the system can be studied. For instance, this allows the self-consistent 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.