Electrocaloric effect in cubic Hubbard nanoclusters

In the paper a computational study of the electrocaloric effect is presented for a cubic nanocluster consisting of 8 sites. The system of interest is described by means of an extended Hubbard model in external electric field at half filling of the energy levels. The thermodynamic description is obtained within canonical ensemble formalism on the basis of exact numerical diagonalization of the system Hamiltonian. In particular, the entropy and the specific heat are determined as a function of temperature and external electric field. The electrocaloric effect is described quantitatively by isothermal entropy change. The behaviour of this quantity is thoroughly analysed as a function of extended Hubbard model parameters, temperature and electric field variation magnitude. The existence of direct and inverse electrocaloric effect is predicted for some range of model parameters. A high sensitivity to Hubbard model parameters is shown, what paves the way towards controlling and tuning the effect. A non-linear, quadratic dependence of isothermal entropy change on electric field variation magnitude is demonstrated. The potential for applications of electrocaloric effect in strongly correlated nanoclusters is shown.

external field, whereas the screening effects are negligible. The model cube, with the atoms situated in the corners, is an ideal candidate for such studies, since in this case all the atoms are located on the surface of the cluster, and the strongest interaction with the external field can be expected.
As a consequence, the aim of our study is theoretical characterization of the electrocaloric effect in a model nanostructure, being a zero-dimensional nanocluster of cubic shape. It consists of N = 8 atoms, as shown schematically in Fig. 1. The bond length defining the distance between nearest neighbour atoms in the cube is a. The behaviour of the electrons in the system is described by means of the extended single-orbital Hubbard model. The presence of n = N = 8 electrons is assumed, so the half-filling condition for the electronic levels is satisfied. The electrons occupy only a single band. The cube is placed in external electric field E directed along one of the bonds -parallel to z direction in the scheme (Fig. 1). The values of z = 0 and z = a correspond to the lower and upper side of the cube, respectively. The extended Hubbard model for the system in electric field has the following Hamiltonian: ) The fundamental Hubbard model is parametrized only by the energy of electronic hopping which takes place between nearest neighbours (hopping integral t 1 ) and by the energy of on-site coulombic interactions between electrons U. However, the extended Hubbard model used here contains additional parameters. They are: the hopping integral for electronic hopping between the second nearest neighbours t 2 , as well as the energy of coulombic repulsion for nearest-neighbour electrons V. All the mentioned quantities are schematically marked in Fig. 1.
The uniform constant electric field E is introduced by its potential term eEz t / i 1 , in which z i denotes the coordinate of the site i along z direction (see Fig. 1), while e is the elementary charge. The dimensionless potential is expressed in t 1 units. For the dimensionless unit, eEa/t 1 = 1, where a is the lattice constant, the approximate value of the electric field E can be estimated as 1 V/Å. In that context, it can be mentioned that the studies of nanoclusters for this range of fields are present in the literature 58 .

Results
Entropy and specific heat. The most fundamental thermodynamic quantity from the point of view of description of any caloric effect is the entropy and its dependence on the temperature. The entropy is defined as, = − ∂ ∂ S G T ( / ) E where G is the Gibbs free energy. Therefore, we commence the presentation of the numerical data with the plot of normalized entropy as a function of the temperature and the on-site coulombic interaction parameter, which is shown in Fig. 2. This contour plot presents the results obtained for the case of the Hubbard model without extensions and in the absence of the electric field. Please note the logarithmic temperature scale used to emphasize the low-temperature behaviour. The limiting value of reduced entropy at T = 0 is equal to 0, while in the high temperature limit it is equal to = . N ln ln12870 9 463 s . It is visible that for low on-site coulombic parameter the entropy increase between the low value (close to zero) and high value (close to the limiting one) takes place in a rather limited interval of moderate temperatures. The effect of increasing U is seen as the widening of this interval of temperatures, so that the isolines of constant entropy shift towards lower temperatures for low entropy and show the opposite behaviour for high entropy as the on-site coulombic parameter increases. In order to emphasize some more subtle features of the temperature dependence of entropy, it is instructive to refer to the specific heat C E defined by Eq. (7), i.e., 2 , as the derivative is more sensitive in this case. Moreover, C E is an interesting quantity itself as a thermal response function. The effect of the on-site coulombic parameter on the temperature dependence of the specific heat can be followed in Fig. 3, again in the form of the contour plot. The case of Fig. 3(a) corresponds to the absence of the electric field. For low U a pronounced single peak of specific heat is observed at moderate temperatures, which is correlated with the interval of the fastest entropy change. However, if U is increased, the peak tends to lose its magnitude and becomes shifted towards lower temperatures. For sufficiently strong on-site coulombic energy a second maximum emerges, which becomes shifted to higher temperatures under the influence of increasing U. This structure with two peaks exhibiting the opposite tendencies to shift when the on-site coulombic parameter changes, explains the behaviour of entropy discussed in the context of Fig. 2, when the temperature interval in which entropy varies most significantly expands with increase of U.
A similar behaviour of the specific heat can be followed in Fig. 3(b), which illustrates the case of the non-zero electric field, = eEa t / 2 1 . The emerging double-peaked structure of the temperature dependence of C E is conserved in the presence of E. However, the range of low U is mostly affected, because the magnitude of specific heat maximum is reduced there. Moreover, the maximum of C E as a function of U, occurring at moderate temperatures, is shifted to non-zero value of U-parameter. The high-U behaviour of the specific heat is rather weakly influenced by the presence of the electric field.
It is instructive to follow the temperature variability of specific heat as a function of the electric field, for fixed coulombic on-site parameter. Such an analysis can be performed on the basis of Fig. 4, which presents the case of unmodified Hubbard model with weak on-site parameter, U/t 1 = 2. For these parameters and low fields E, a single peak of specific heat occurs at intermediate temperatures. Application of the electric field in general tends to reduce the peak magnitude, but this effect is not very pronounced for weak fields. For stronger E the maximum of C E is shifted somehow towards higher temperatures and an additional branch being a weak second maximum  , the influence of the electric field is even less pronounced, with a general tendency to reduce the magnitude of specific heat peaks. The second, weak maximum, which is seen in the low-temperature region in the vicinity of the field eEa t / 45 1 ≈ . does not exist in the limit k T t / 0 3 for the set of Hamiltonian parameters in Fig. 4. It has been checked that in the limit k T t / 0 B 1 → the specific heat C E always tends to zero value, in accordance with the third law of thermodynamics. This particular behaviour of the specific heat can be interrelated with the ground-state properties of the system. In Electronic Supplementary Material we provide Fig. S1, which presents a dependence of the energy gap between the ground state and the first excited state on the electric field (with an inset showing an analogous dependence of the individual energy states of interest). At the point eEa t / 45 1 ≈ . the gap closes and for stronger electric fields it remains rather small, moreover, the ground state acquires three-fold degeneracy.
In order to verify the influence of the energy of coulombic interactions between electrons localized at nearest-neighbour sites V on the specific heat, we include Fig. S2 in Electronic Supplementary Material, which is a contour plot of C E as a function of the temperature and normalized repulsion V between nearest neighbours, in the absence of electric field, for U/t 1 = 2.
Isothermal entropy change. Let us commence the characterization of electrocaloric properties from the systematic analysis of the effect of coulombic interaction U on the isothermal entropy change, , when the electric field is varied between fixed non-zero value E and zero value. The isothermal entropy change is usually considered as a quantity characterizing the vertical size of Ericsson cycle 9 , when this cycle is presented in (T,S) coordinates, and therefore is very important for the description of electrocaloric effect, being also used to define refrigerant capacity 7,9 . This quantity is also interesting from the point of view of its dependence on the electric field strength E. Let us only mention that alternatively, the derivative ∂ ∂ S E ( / ) T could also be studied. Both quantities are related by . It can be mentioned that the entropy change for finite electric field variation magnitude can be directly measured for example with differential scanning calorimetry technique 7,12,59 . First, an unmodified Hubbard model will be considered, i.e. with t 2 = 0 and V = 0. In order to visualise the effect we present a contour plot in Fig. 5, showing the isothermal entropy change in the temperature-coulombic interaction energy plane (please note the logarithmic temperature scale). In general, it can be noticed that both ranges of direct ECE (with ΔS T > 0) and inverse ECE (corresponding to ΔS T < 0) are present in our system of interest. One inverse ECE range emerges at low temperatures for rather weak coulombic interactions. Another, much more extensive range is present at intermediate temperatures and it appears only when U/t 1 exceeds a certain critical value. If the coulombic interaction rises further, the range shifts toward higher temperatures, but the magnitude of inverse ECE achieves a maximum around U t / 6 1 and then falls down. The rest of the phase diagram is filled with the range of direct ECE, which is always present at high enough or low enough temperature (excluding the range of low values of U). However, at a certain range of coulombic interaction energy, only the direct ECE is found at arbitrary temperature. The direct ECE magnitude is the largest at intermediate temperature and for weak or zero U/t 1 . Another maximum (smaller) is present at low temperatures and U/t 1 close to 4. What is important, the extremal magnitude of direct ECE is about 3 times larger than the corresponding magnitude of inverse ECE.
As it is seen from the above analysis, the electrocaloric effect is rather sensitive to the energy of the coulombic interactions, which are capable of modifying severely its magnitude as well as of switching between the direct and inverse ECE at given temperature, and also of moving the position of its extrema. First the effect of inclusion of t 2 will be discussed. In Fig. 6(a) the evolution of the thermal dependence of entropy change is shown for low value of on-site coulombic parameter U/t 1 = 2. Under such conditions, only direct ECE is present for t 2 = 0 (see Fig. 5). It is visible that by introducing t 2 the extremal magnitude of direct ECE tends to increase. Moreover, in the low-temperature range, a region of inverse ECE builds up, with formation of a minimum on ΔS T curve. This minimum shifts towards lower temperatures and increases in magnitude when t 2 rises. Therefore, for U/t 1 = 2 a profound influence of second nearest neighbours hopping on the character of ECE can be expected in the whole range of temperatures.
The situation is rather different for the case of stronger coulombic interactions, for instance, for U/t 1 = 10, depicted in Fig. 6(b). In this case, increase in t 2 barely affects the high-temperature minimum of ECE, just reducing slightly its magnitude without shifting the position. On the contrary, the low-temperature maximum corresponding to direct ECE gradually flattens and even a range of inverse ECE emerges for the lowest temperatures, forming a minimum, which further develops. This proves that for stronger coulombic interactions only the low-temperature ECE can be effectively controlled and tuned with t 2 .
Another parameter of the extended Hubbard model is the energy of coulombic interactions between electrons localized at nearest-neighbour sites V. The influence of this parameter on the form of ECE in the studied system, for low energy of coulombic on-site interactions U/t 1 = 2 and t 2 = 0 is visualized as a contour plot in Fig. 7(a). For Hubbard model without extensions we deal with only direct ECE. Under the influence of increasing V the maximum at intermediate temperatures gradually flattens and shifts somehow to higher temperatures. Then, for  1 the different behaviour emerges -the range of direct ECE is divided into two ranges separated with inverse ECE with a minimum of ΔS T . When V grows further, this minimum develops, whereas the low temperature maximum of direct ECE moves to lower temperatures and the high-temperature maximum of ECE flattens. As a consequence, a profound effect of including V on the isothermal entropy change in ECE is found.
The case of stronger coulombic on-site interactions, for U/t 1 = 10, appears more complex, as it is seen in Fig. 7(b). When weak V > 0 is present, both the low-temperature maximum of direct ECE and the high-temperature minimum of inverse ECE are not greatly influenced by V. However, above V/t 1 = 3, the low-temperature range of direct ECE disappears completely, just after passing through a very pronounced peak. At the same time the high-temperature minimum of inverse ECE becomes shifted towards higher temperatures, what is accompanied with decrease in magnitude.
A crucial feature of the electrocaloric effect is the behaviour of its quantitative measures as a function of the magnitude of electric field change. Such an analysis can be performed on the basis of Fig. 8, for three cases of on-site coulombic interactions, U/t 1 = 2 ( Fig. 8(a)), 5 ( Fig. 8(b)) and 10 ( Fig. 8(c)). In the case of low coulombic interactions, U/t 1 = 2.0 ( Fig. 8(a)), it is visible that for low electric field only a broad range of direct ECE is present. When the field increases, this maximum develops and the temperature at which entropy change reaches the extremal magnitude rises slowly. For even stronger field amplitude, the low-temperature range of inverse ECE builds up, being very sensitive to E. It can be, therefore, seen, that for low coulombic on-site interactions energy the amplitude of electric field not only influences the magnitude of ECE but also changes qualitatively the shape of its temperature dependence. On the contrary, for stronger coulombic on-site interactions (shown in Fig. 8(b) and (c)), only the magnitude of ΔS T is varied when the electric field amplitude changes, while the general shape of the dependence of ΔS T on T is conserved, and no switching between direct and inverse ECE takes place. In that cases a pronounced nonlinearity in the dependence of the entropy change on the electric field can be observed, as ΔS T rises faster than linearly with E.
In order to supplement the physical picture, also the thermal variability of the differential ECE measure, S E ( / ) T ∂ ∂ for various electric fields can be traced, as in Fig. S3 in Electronic Supplementary Material. In general, the curves bear significant similarity to the corresponding ones for the quantity ΔS T (Fig. 8). Moreover, the contour plot of ∂ ∂ S E ( / ) T as a function of temperature and electric field is presented as Fig. S4. In order to get more detailed insight into the variability of ΔS T , let us plot this quantity as a function of the electric field amplitude for a fixed temperature. For this purpose we select U t / 10 0 1 = . and the temperatures corresponding approximately to the extrema of ECE (both direct one, at lower temperature, and inverse one at higher temperature). The dependence is presented in Fig. 8 (in main plot the linear scale is used) and the non-linear increase of the entropy change magnitude with the increase of the electric field can be clearly confirmed. The non-linear variability calls for plotting the data in double logarithmic scale, what is done in the inset to Fig. 9 (where the absolute value of ΔS T is presented). In double logarithmic scale both dependences linearize and it can be extracted from the inset that the dependence is quadratic, that is ∆ ∝ S E T 2 , for the studied ranges of parameters (close to the extrema). This is an important finding, since, for example, in the case of the usually considered magnetocaloric effect, the entropy change is linearly dependent on the magnetic field amplitude.
The above quadratic dependence means that the derivative ∂ ∂ S E / is linear vs. E in the considered range of parameters. In turn, taking into account the Maxwell relation: ∂ ∂ = ∂ ∂ S E P T / / , where P is the electric polarization of the cluster, we see that P should be also a linear function of E. This fact proves that our system is in a range of linear response to the external force. For the discussion of the quadratic dependence of the ECE measure on electric field see also ref. 33 .

Discussion
In the paper a theoretical study of the electrocaloric effect in a cubic nanocluster embedded in the electric field is presented. The system is modelled with the help of extended Hubbard model at half-filling, taking into account coulombic interactions between on-site and nearest-neighbour electrons, as well as electronic hopping between nearest-and next nearest neighbours. For the resulting Hamiltonian an exact numerical diagonalization was performed within canonical ensemble formalism, for fixed number of electrons. Such situation corresponds to chemically isolated system, where the chemical potential is constant and no exchange of electrons with the environment takes place. It should be mentioned here that for nanoscopic systems interacting chemically with their environment, for instance, when the cluster interacts with the conducting substrate, the grand canonical ensemble should be more appropriate. Both ensembles are applied in the studies of nanoscopic Hubbard systems and they are inequivalent 60,61 . When the statistical fluctuations of the number of electrons are taken into account within grand canonical ensemble, for small clusters, the discussed properties are additionally influenced. We might state that both ensembles correspond to manifestly different physical situations and their comparison would inspire further works. The thermodynamic description included in particular such quantities as the specific heat and entropy, which are the functions of the external electric field and temperature. In order to characterize the ECE, the isothermal entropy change was thoroughly studied as a function of temperature, Hubbard model parameters and electric field amplitude. Significant ranges of both direct and inverse ECE were identified. Moreover, significant sensitivity of the isothermal entropy change to Hubbard model parameters was shown. This enables the design and control of the electrocaloric properties of the nanoclusters.
It should be mentioned that having obtained the entropy as a function of T and E, the change of temperature (ΔT) S can be studied, when the field E is switched on (off) adiabatically, i.e., for S const = . The parameter (ΔT) S has a meaning for the prediction of cooling (heating) effect in the system, when rapid changes of the external field take place. However, in this paper we concentrated rather on the comprehensive investigations of the (ΔS) T and C E quantities. We would like to emphasize that the knowledge of the isothermal entropy change allows the determination of the heat transfer to or from the system, which is the vital process during the electrocaloric refrigeration.
It should be particularly emphasized that the obtained results are exact, since they originate from the exact numerical diagonalization of the many-body extended Hubbard Hamiltonian for the system in question. Such  an approach is computationally demanding, what limits its applicability to small systems, however, it provides a rigorous description of the thermodynamic properties, capturing all the many-body effects.
In this paper we concentrate on the rigorous solutions of the extended Hubbard model for small clusters when the electric field E is applied. In this model approach the influence of the electric field on the Hubbard Hamiltonian parameters was neglected, and these parameters are considered to be constant. This is the usual approach (e.g. 62,63 ). However, it should be noted here that for the electric field with the potential changing markedly on the scale of lattice constant, the atomic states which give rise to the Hubbard model could be deformed in comparison with the states without field. The example of the study when the atomic 1s states were modified to follow the varying interatomic distance parameter was given, for example, in ref. 64 . As a result, this deformation under the action of the electric field will influence both the hopping integrals as well as the Coulomb interaction potentials. Therefore, one should be aware that for very strong electric fields, the Hubbard model with constant interaction parameters can serve only as an approximate description of the physical reality, being, on the other hand, a starting point for more sophisticated treatments of the problem.
The future developments may include, for example, the consideration of various cluster sizes and shapes, as well as larger structures (with the help of modified numerical approach to the Hubbard model). Another potentially fruitful direction involves the study of charge doping on the ECE in nanostructures. Moreover, another external field such as magnetic field could be incorporated to investigate multicaloric properties 65 . In addition, lattice-related degrees of freedom might be also taken into account to exploit possible coupling to lattice entropy.

Methods
In order to obtain the full thermodynamic description of the system, an exact diagonalization method is used, allowing the treatment of the Hubbard terms in the Hamiltonian without any approximations 66  and the corresponding eigenvectors. The exact thermodynamic description of our system is constructed within the canonical ensemble, with the external electric field taken into account, assuming fixed number of electrons n and the fixed temperature equal to T. Under such assumptions, the statistical sum is given by: The knowledge of the statistical sum allows the direct calculation of the Gibbs free energy: On the other hand, the enthalpy of the system is given as the thermal average of the Hamiltonian: Let us put emphasis on the fact that the Hamiltonian given by Eq. (1) contains the potential energy term for the electrons in the electric field, therefore, its average H is not a sole internal energy but enthalpy instead. For the same reason the free energy determined directly from the statistical sum is the Gibbs free energy and not the Helmholtz one.
The knowledge of both quantities i.e., G and H, yields the entropy value, which is the most interesting quantity for us in the present study. Namely: The entropy S T E ( , ) is then given exactly for the system of interest, as a function of the temperature and external electric field, for the whole range of the extended Hubbard model parameters.
The knowledge of entropy allows, among other things, the study of the electrocaloric effect (ECE). Its crucial measure is the isothermal entropy change when the external electric field varies at constant temperature from the value E = 0 to E > 0, i.e.

S S T E S T E
T ∆ = = − .
If the value of ΔS T is positive, so the electric field decreases the entropy of the system, the direct ECE occurs. If the opposite case takes place, it is a manifestation of inverse ECE.
Another quantity of interest may be the specific heat of the system for constant electric field = ∂ ∂ C T S T ( / ) E E , which can be conveniently calculated on the basis of fluctuation-dissipation theorem, namely: