# Classical Analogue of the Ionic Hubbard Model

###### Abstract

In our earlier work [M. Hafez, et al., Phys. Lett. A 373 (2009) 4479] we employed the flow equation method to obtain a classical effective model from a quantum mechanical parent Hamiltonian called, the ionic Hubbard model (IHM). The classical ionic Hubbard model (CIHM) obtained in this way contains solely Fermionic occupation numbers of two species corresponding to particles with and spin, respectively. In this paper, we employ the transfer matrix method to analytically solve the CIHM at finite temperature in one dimension. In the limit of zero temperature, we find two insulating phases at large and small Coulomb interaction strength, , mediated with a gap-less metallic phase, resulting in two continuous metal-insulator transitions. Our results are further supported with Monte Carlo simulations.

###### pacs:

71.30.+h, 68.35.Rh## I Introduction

To understand the magnetism and magnetic phenomena, one of the basic interactions is the exchange mechanism which is deeply rooted in Coulomb interactions and quantum mechanical indistinguishability. Therefore a fair understanding of the magnetic behavior of materials is not possible without investigating appropriate quantum spin models. Introducing uni-axial anisotropy to the Heisenberg model amounts to suppression of transverse quantum fluctuations (), leading to the so called Ising model Newell (). The resulting Ising Hamiltonian turns out to contain the basic magnetic phases of the original Heisenberg model, namely ferromagnetism, and antiferromagnetism Feynman (). Although many of the interesting possible aspects such as spin liquid phases, spin-wave excitations, etc. Fradkin () can not be captured by the Ising model. Classical Ising model has the merit of being much simpler to solve, and admits analytical Newell (); Onsager () and graphical solutions Feynman () in various geometries JafariIJMP () in one and two dimensions which is lacking in the original quantum Heisenberg model.

Dielectric properties are among the most important properties characterizing materials. The question can be asked here, is there any Ising-like model that can provide basic informations about their phase diagram, and at the same time being simple enough to allow for analytical solutions? We have taken the example of the (IHM) Nagaosa (); Egami (). This model was introduced to study the neutral-to-ionic transition in organic compounds, as well as, understanding the role of strong electronic correlations in the dielectric properties of metal oxides Egami (); HafezCUT (). This model is as follows:

(1) | |||||

where () is the usual annihilation (creation) operator at site with spin . is the on-site Coulomb interaction parameter, and denotes a one-body staggered ionic potential. The kinetic energy scale is given by the real hopping amplitude which prefers to gain kinetic energy by spreading the wave-function over the whole system, leading to quantum fluctuation of the charge density. The zero temperature phase diagram of this model contains Mott and band insulating states when the energy scales corresponding to or dominate, respectively HafezCUT (). When these two scales become of the same order of magnitude, the nature of intermediate phase still remains controversial. Some authors argue that the intermediate state is a spontaneously broken symmetry phase SBSP (), while some others argue that the phase in between is metallic Nagaosa (); Bouadim2007 (); HafezCUT (); Metals (). In our previous investigation HafezCUT () we employed the method of flow equations for the quantum Hamiltonian (1) to obtain an effective Hamiltonian in which the hopping term is renormalized to zero, producing a new nearest neighbor Coulomb attraction. The resulting classical Hamiltonian is of the following lattice gas form:

(2) |

where the renormalized parameters have closed form expressions in terms of the bare parameters HafezCUT (). We take the kinetic energy scale as the unit of energy.

By calculating the spin and charge gaps, we showed that this simple classical model is capable of capturing the physics of a metallic state sandwiched between two distinct insulating phases at zero temperature as one increases for a fixed value of HafezCUT (). Here contains two species (or colors) corresponding to , respectively. Since for each ”color” is either or , it is an Ising-like variable. Therefore our effective CIHM Hamiltonian can be thought of, as two inter-penetrating Ising models on a lattice. Being a natural extension of two species lattice gas model, allows for analytical solution in one spatial dimensions (1D). Here we can employ the transfer matrix method, constructed in terms of matrices to calculate the thermodynamic properties of this model at finite temperatures. At high temperatures where the thermal fluctuations wash out the quantum effects, we expect the results of our calculations to be a good description of the class of materials modeled in terms of IHM. Also a novel graphical solution in 2D similar to Feynman’s construction can be worked out JafariUnpublished ().

The paper is organized as follows: In the next section we calculate the grand canonical potential for this model in one dimension and discuss the particle density and ionicity in various fillings. It is followed by a section focused on half-filling situation, and calculate the specific heat, compressibility and ionicity to assess the nature of the phases in the parameter space of . Throughout the paper we have fixed the value of HafezCUT (). The final section is devoted to summary and discussion.

## Ii Grand Potential

In this section we calculate the grand canonical potential that can be derived from the grand partition function (GPF) which is defined as follows:

(3) |

where is Temperature, is lattice size, and is chemical potential. The denotes all possible configurations of the occupation numbers, which must be summed over. is the number of particle and is the energy that is defined by the model Hamiltonian Eq. (2).

The summation in Eq. (3) can be calculated analytically as follows. The values of are only zero and one, which is the Fermionic memory left in the commuting variables . Hence , and the second term in Eq. (2) becomes:

(4) |

The GPF becomes,

(5) | |||||

As can be seen at the Hamiltonian level, only the spin-summed occupations appear and therefore this classical model does not contain spin-polarized (magnetic) solutions. So we change the summation variable from and to that has three possible values 0,1,2. Therefore we have:

(6) | |||||

where the coefficient takes into account the two-fold degeneracy for in Eq. (5), which corresponds to either , , or , . Eq. (6) can be written in a more symmetric form,

(7) | |||||

where is assumed to be even and the periodic boundary conditions, is implied. Defining the matrix elements of the transfer matrix by,

(8) | |||||

Eq. (7) can be written as:

(9) | |||||

where is the transpose of , . Then the partition becomes:

(10) |

where , , and are the eigenvalues of . In the thermodynamic limit where is very large, the grand potential per site becomes:

(11) |

Here is the maximum eigenvalue. Hence to obtain the grand potential one needs to calculate the eigenvalues of . Later on, when we discuss the ionicity, we will in addition need the eigenvectors of . The eigenvalues of are solutions to the third order equation,

(12) |

with , , and given in the appendix. Through the paper we will report the numerical plots for in units in which . For this value of , there would be no level-crossing among the eigenvalues of the matrix when one varies , , and , as expected from Perron’s theorem Goldenfeld ().

From the partition function, one can in principle calculate various averages of the form . The simplest of these averages are the average particle density (symmetric combination), and the ionicity (antisymmetric combination) which contain important information about the nature of the thermodynamic phases of the model. In the following we calculate these averages as a function of .

### ii.1 Particle Density

Once the grand potential is known, one can derive many other thermodynamic quantities. The particle density, , can be calculated as,

(13) |

In Fig. 1, the particle density versus has been plotted for different values of . In this figure the value of is fixed to be . This value of corresponds to a band insulating regime at zero temperature HafezCUT (). As can be seen in the figure, apart from trivial cases of the empty () and the fully-filled () lattice, there are three plateaus. The first one corresponds to the half-filling (i.e. one particle per lattice site), and two others correspond to commensurate fillings, (quarter-filling), (three-quarter-filling) which is quite similar to the plateaus of the parent ionic Hubbard model Bouadim2007 (). Indeed the and are related by a particle-hole transformation. Therefore in what follows, we focus on . With increasing the temperature, the plateaus get more and more rounded. The value of chemical potential corresponding to the half-filling is independent of temperature, as can be seen in Fig. 1. The inset plot indicates this point more transparently. This ”isosbestic” behavior is observed for all values of values of .

In Fig. 2, we have plotted the particle density as a function of chemical potential at a constant (low) temperature , for different values of . As can be seen in the figure, increasing , causes the plateaus which are particle-hole counterpart of each other, get wider. However the plateau at gets narrower and finally vanishes around , resulting in a gap-less phase at half-filling. Upon further increasing , the half-filling plateau is recovered, and gets wider, indicating the emergence of a growing new gap in the system, which is reminiscent of the Mott insulating behavior. To identify the nature of gap at one needs to calculate the ionicity, which provides information about how the unit cell is being filled. This will be done in the sequel.

For later reference, we report the value of chemical potential corresponding to half-filling, which by examining the numerical plots turns out to be:

(14) |

The temperature range at which the above relation is valid, is roughly . For temperatures outside this range, there will be deviations from the above rule of thumb relation.

### ii.2 Ionicity

If denote the sublattice of odd and even sites, respectively, then and are defined as:

(15) | |||||

(16) |

the ionicity per site becomes , where , and is the total number of lattice sites. This quantity can also be calculated analytically with the aid of transfer matrices. In the grand canonical ensemble is given by:

(17) |

Depending on whether is odd or even, similar to calculations for grand partition function, we obtain:

(18) | |||||

(19) |

where is the following matrix:

(20) |

In Fig. 3, we have plotted the ionicity per site as a function of for various temperatures at a fixed value of . For this value of , at zero temperature one would expect a band insulator at half-filling, where A-sublattice (with on-site energy ) are doubly occupied. Let us focus around which corresponds to quarter filling (c.f. Fig. 1). The value of for in Fig. 3 shows that the only electron of each unit cell belongs to A-sublattice, which represents a charge density wave insulator.

As we increase from , the ionicity increases, which indicates that sublattice A continues to be filled. When , one approaches the half-filling (Fig. 1), where at lower temperatures the A-sublattice is totally occupied, hence . As can be seen in Fig. 3, this saturation value is decreased as the temperature is increased. This reduction in the ionicity, indicates that the particles are bing thermally excited across the band gap. Further increasing of the chemical potential, one reaches the plateau around of Fig. 3 (c.f. Fig. 1). The decrease in the ionicity indicates that the added particles essentially start to occupy the B-sublattice. The symmetry of Fig. 3 around , is due to the apparent particle-hole symmetry of the Hamiltonian.

In Fig. 4 we plot the ionicity per site at the constant temperature and different values of indicated in the legend. Again there are three plateaus corresponding to 1/4,1/2 and 3/4-filling, respectively. This figure indicates that the center of the half-filling plateau is at . This can be understood as a Hartree like energy for the Hubbard model. Now concentrating around half-filling in this figure, for small values of , the ionicity reaches , which indicates a complete electric polarization of the unit cell, and hence we have a band insulator. For intermediate values of , the maximum value of ionicity does not reach , which shows some of the added particles start to occupy the B-sublattice. Also note that for the intermediate values of the half-filling plateau starts to disappear, i.e. the emergence of a metallic behavior. For large values of , the half-filling plateau is restored, but at zero ionicity. Therefore we have an insulator with unpolarized unit cell. Such an insulating state can be thought of as classical counterpart of the Mott insulating state.

## Iii Half-filling

As we saw in Fig. 2, for the energy scales cooperate with each other, to give rise to a charge density wave insulating ground state. Rather at , these two energy scales compete against each other to destroy the insulating behavior for intermediate (), giving rise to a richer phase diagram. Therefore in this section we focus at half-filling and calculate the specific heat, compressibility, and ionicity. Before doing so, we compare some physical quantities evaluated by a fully numerical Monte Carlo simulation, with our exact transfer matrix results. In Fig. 5, ionicity and specific heat per site are plotted at half-filling and show a good agreement between analytical and numerical results. This ensures that both Monte Carlo and analytic results are reliable.

Now let us proceed with the calculation of various thermodynamic quantities. For the fixed value of , we have two ways to plot thermodynamic quantities. First way is to plot them as a function of temperature , at some selected values of . These results indicate that in the present one dimensional model, there will be no finite temperature phase transition. This is obviously due to the analytic behavior of the partition function as a function of . The second way is to plot them as a function of , for some selected temperatures. This second way of presenting the data, reveals that as one lowers the temperature, there will be sharper features as a function of , indicating the zero temperature phase transition. Our calculations are based on Eqs. (11,14). For very low temperatures, where the validity of Eq. (14), might be questioned, we employ Monte Carlo simulation data.

### iii.1 Specific heat

Specific heat per site can be calculated as:

(21) |

where is entropy per site at half-filling that can be derived from the grand potential. Figure 6 shows the specific heat versus for various values of . As can be seen in Fig. 6 for values of , there is a single peak in the . For a peak-dip-hump structure can be observed. For the hump is quite clear, while for , it can be interpreted as a precursor to the hump. For , the peak merges into the hump structure. In terms of the parent quantum Hamiltonian, such hump structure indicates incoherent excitations. For , Eq. (14) is not reliable at very low temperatures. Therefore we report Monte Carlo simulation results which indicates a very sharp peak at . As moves towards from both sides (lower and upper panels), the peak gets sharper and moves towards lower temperatures. This indicates that in limit one expects a transition between gaped and gap-less states.

This behavior can be simply understood in terms of a two-state model with level spacing , whose specific heat is given by

(22) |

where . Behavior of Eq. (22) for is like , while for it vanishes as . For the intermediate region a Schottky peak around () arises in the specific heat. Fitting the specific data to Eq. (22), in Fig. 7 we have plotted the estimated gap versus , which indicates two gaped phases. According to the above two-state formula, this peak corresponds to , from which a gap of can be estimated. If we extrapolate the estimated gap for , zero gap region is expected to occur for which is compatible with our previous work HafezCUT ().

Now let us look at the specific heat data from a different angle. In Fig. 8 we plot the specific heat versus for different values of temperatures. As can be seen, there are two peak structures at all temperatures, which get sharper and by lowering the temperature, they tend to accumulate around . Extrapolating the trend of this double-peak structure to the limit of , suggests two phase transitions at and HafezCUT (), compatible with the behavior of the vanishing gap region in Fig. 7. The characteristic quadratic behavior around seen in Fig. 8, which according to the two-state model is expected to be like , indicates two continuous metal-insulator transitions, with , .

### iii.2 Compressibility

Another quantity that can be treated in our consideration CIHM is the compressibility which is given by:

(23) |

where is the density of particles per site.

In Fig. 9 we plot the compressibility as a function of for different values of . Zero compressibility is a characteristic of gaped states. As can be seen in this figure for small value of the range of temperatures at which the compressibility is close to zero is substantial, which means that the gap is so large that up to such temperature the insulating behavior is still manifest. By increasing , this temperature range shrinks and becomes smaller and smaller, until around , it extrapolates to zero. Increasing beyond , again recovers a finite temperature range in which the compressibility is zero. This behavior confirms that a gap-less state is sandwiched between two gaped states.

To see the above statement more clearly, in Fig. 10 we plot compressibility as a function of for selected temperatures.

As can be seen in this figure, there is a region with non-zero compressibility, which characterizes a gap-less phase. Outside this region, the compressibility decays to zero. By decreasing , the width of the compressibility peak becomes smaller and while the height of the peak diverges as ; a typical characteristic of a continuous phase transition. This confirms the existence of a metallic phase at zero temperature HafezCUT (). The effect of thermal fluctuations is to smear the edges of metallic region. This is quite intuitive, as for values of near the zero temperature boundary of metallic phase with neighboring insulating phases, the gaps are small, and hence the thermal energy can overcome the gap.

### iii.3 Ionicity

In Fig. 11 we plot the ionicity per site for and various values of at half-filling. As can be seen in the figure, by lowering the temperature, for , the ionicity tends to ; for it reaches a zero temperature value of ; while for , it reaches a value between these two limits; characterizing a phase with charge fluctuations.

This behavior can be understood as follows: At small regime the unit cell is fully polarized at low temperatures, with both and particles occupying the A-sublattice. As the temperature is increased, some of the particles get excited to B-sublattice by absorbing the thermal energy, . Similarly for large values of , at lower temperatures the unit cell is not polarized, due to the Coulomb term . By increasing the temperature, thermal excitations with doubly occupied sites will be created, thereby increasing the ionicity. For intermediate , the weights of polarized and unpolarized configurations in the unit cell become comparable; hence giving the ionicity .

In Fig. 12 we plot the ionicity at half-filling as a function of for and various values of . As can be seen at all temperatures, the ionicity smoothly varies between for small values , and for large values of . The width of the transition region decreases by lowering the temperature, and is expected to approach the figure 2 of Ref. HafezCUT, .

## Iv Conclusion

In this work we studied a classical model consisting of two Ising like variables on a one dimensional chain. Despite the simplicity which essentially results form the lack of Fermionic minus sign problem, our model captures some of the interesting properties of the ionic Hubbard model. Various thermodynamic properties, such as, specific heat, ionicity, particle density and compressibility when viewed as a function of in a given temperature, indicate the presence of two gaped states at small and large values of , with a gap-less state sandwiched between them (around ). When the same quantities viewed as a function of temperature, there is no sign of phase transition down to zero temperature. The three phase scenario of the zero temperature with clear zero temperature phase transition boundaries at and extrapolates to higher temperatures. However, the boundaries get smeared due to thermal fluctuations, giving rise to a cross-over behavior.

Mapping of dimensional quantum models to dimensional classical Hamiltonians is a well known paradigm in statistical physics. Our flow equation approach HafezCUT () suggests an alternative approach to construct ”” dimensional classical models which might be useful in capturing basic aspects of the original quantum Hamiltonian.

## V Acknowledgements

We would like to thank Professor N. Nafari for useful comments and critique of the manuscript. This research was supported by the Vice Chancellor for Research Affairs of the Isfahan University of Technology (IUT). S.A.J was supported by the National Elite Foundation (NEF) of Iran.

## Vi Appendix

### vi.1 Coefficients of Eq. (12)

## References

- (1) G.F. Newell and E.W. Montroll, Rev. Mod. Phys. 25, (1953) 353.
- (2) R.P. Feynman, Statistical Mechanics, Westview Press, 1998.
- (3) E. Fradkin, Field Theories of Condensed Matter Systems, Westview Press, 1998.
- (4) L. Onsager, Phys. Rev. 65, (1944) 117.
- (5) S. A. Jafari Int. J. Mod. Phys. B 23, (2009) 395.
- (6) N. Nagaosa, J. Takimoto, J. Phys. Soc. Jpn. 55 (1986) 2735.
- (7) T. Egami S. Ishihara, M. Tachiki Science 261 (1993) 1307.
- (8) M. Hafez, S. A. Jafari M. R. Abolhassani Phys. Lett. A, 373 (2009) 4479.
- (9) M. Fabrizio, A. O. Gogolin and A. A. Nersesyan, Phys. Rev. Lett. 83 (1999) 2014; C. D. Batista and A. A. Aligia, Phys. Rev. Lett. 92 (2004) 246405; S. R. Manmana, V. Meden, R. M. Noack, and K. Schönhammer, Phys. Rev. B 70 (2004) 155115; Ö. Legeza, K. Buchta, and J. Sólyom, Phys. Rev. B 73 (2006) 165124; T. Wilkens and R.M. Martin, Phys. Rev. B 63 (2001) 235108; S. S. Kancharla and E. Dagotto, Phys. Rev. Lett. 98 (2007) 016402.
- (10) K. Bouadim, N. Paris, F. Hébert, G. G. Batrouni and R.T. Scalettar, Phys. Rev. B 76 (2007) 085112;
- (11) A. Garg, H. R. Krishnamurthy, and M. Randeria, phys. Rev. Lett. 97 (2006) 046403; N. Gidopoulos, S. Sorella, E. Tosatti E. Phys. J. B 14 (2000) 217; N. Paris, K. Bouadim, F. Hébert, G. G. Batrouni and R. T. Scalettar, Phys. Rev. Lett. 98 (2007) 046403;
- (12) S. A. Jafari et al., unpublished.
- (13) N. Goldenfeld, Lecture on Phase Transitions and The Renormalization Group, Westview Press, 1992.