text:
# Hubbard Hamiltonian in the dimer representation. Large 𝑈 limit

## 1 Introduction

The single-band Hubbard model, Ref. , plays in the solid state physics a similar principal role as the hydrogen atom in the atomic physics. This explains a continuous interest in its properties. The Hubbard Hamiltonian reads

$$H=\underset{i,j,\sigma }{}t_{i,j}c_{i,\sigma }^+c_{j,\sigma }+U\underset{i}{}n_{i,}n_{i,}.$$
(1)
Here $`c_{i,\sigma }`$ $`(c_{i,\sigma }^+)`$ are annihilation (creation) operators of an electron with spin $`\sigma =`$, $``$ in the Wannier representation at the lattice site $`𝐑_i`$ and $`n_{i,\sigma }=c_{i,\sigma }^+c_{i,\sigma }`$. Moreover, $`t_{i,j}`$ is the hopping integral between different lattice sites $`i`$ and $`j`$ ($`t_{i,i}=0)`$ and $`U`$ is the intrasite Coulomb repulsion. The Bloch conduction band energy $`\epsilon _𝐤`$ is given by

$$\epsilon _𝐤=\underset{ij}{}t_{i,j}e^{i𝐤(𝐑_i𝐑_j)}.$$
(2)
In the following we restrict ourselves to the simple cubic (sc) lattice and assume that

$$t_{i,j}=\{\begin{array}{cc}t\hfill & i,j\text{-nearest neighbours}\hfill \\ 0\hfill & \text{otherwise}.\hfill \end{array}$$

Then

$$\epsilon _𝐤=2t(\mathrm{cos}k_xa+\mathrm{cos}k_ya+\mathrm{cos}k_za).$$
(3)
The hopping parameter $`t`$ is simply related to the bandwidth $`W`$ of the Bloch band (2), e. g. $`W=12t`$ for the sc lattice. The interplay between the two model parameters, $`W`$ and $`U,`$ is decisive for the properties of the model resulting in strong electron correlations, leading to band magnetism (see e.g. Refs ,  for a review), insulator-to-metal transition (see e.g. Refs ,  and papers cited therein) and high-$`T_c`$ superconductivity (negative $`U`$-model, see e.g. Ref. ). The Hubbard model very often plays also a role of a submodel for many other more complicated models (as e.g. Anderson model (see Ref. ), s-f model (see Refs ) and so on). Especially interesting, but difficult to handle are the properties of the Hubbard model in such a regime of the model parameters where the bandwidth $`W`$ is comparable to the Coulomb repulsion $`U`$. In a large number of papers \[9-30\] the authors tried to solve this model using many sophisticated methods. The exact solution, however, does not exist till now. Many authors tried to change the situation in this field by introducing the expansion parameter $`x=\frac{t}{U}`$ $`(x1)`$. This idea (cf Refs ) consists in replacing ”difficult physics” connected with the model by ”difficult mathematics” obtained by a laborious expansion with respect to $`x.`$ Different methods connected with this problem have been applied as e.g. the perturbation expansion (see Refs , ), canonical transformation (see Refs , , , , ) or ab initio derivations (see Refs , \[28-30\]). Most of the methods, leading to the $`tJ`$ model (or generalized $`tJ`$ model) are also summarized in Refs , .

The goal of the present paper is just to show that we can take another, nonperturbative way. In the first step we divide the crystal lattice into a set of interacting dimers. In other words, we can rewrite the Hubbard Hamiltonian (1) for the sc lattice (see Fig. 1) in the equivalent form 111The dimer Fourier transformation $`c_{I,1,\sigma }=\frac{1}{\sqrt{N}}\underset{𝐤}{}c_{𝐤,\sigma }e^{i𝐤𝐑_{I,1}},c_{I,2,\sigma }=\frac{1}{\sqrt{N}}\underset{𝐤}{}c_{𝐤,\sigma }e^{i𝐤𝐑_{I,2}}`$where $`N`$ is the number of lattice points, applied to (4), gives the well known result $`H=\underset{𝐤,\sigma }{}\epsilon _𝐤n_{𝐤,\sigma }+\frac{U}{N}\underset{𝐤,𝐤^{},𝐪}{}c_{𝐤+𝐪,}^+c_{𝐤,}c_{𝐤^{}𝐪,}^+c_{𝐤^{},}`$ with $`\epsilon _𝐤`$ given by (3).:

$$\begin{array}{cc}H=\hfill & \underset{I}{}H_I^dt\underset{I,\sigma }{}(c_{I,2,\sigma }^+c_{I+1,1,\sigma }+c_{I+1,1,\sigma }^+c_{I,2,\sigma })\hfill \\ & t\underset{IJ,\sigma }{}(c_{I,1,\sigma }^+c_{J,1,\sigma }+c_{I,2,\sigma }^+c_{J,2,\sigma })\hfill \end{array}$$
(4)
where

$$H_I^d=t\underset{\sigma }{}(c_{I,1,\sigma }^+c_{I,2,\sigma }+c_{I,2,\sigma }^+c_{I,1,\sigma })+U(n_{I,1,}n_{I,1,}+n_{I,2,}n_{I,2,}).$$
(5)
The indices $`I`$ and $`J`$ enumerate the dimers and $`H_I^d`$ is the dimer Hamiltonian. The second term in (4) describes the hopping between nearest dimers in the $`z`$-direction, the third one represents the hopping between nearest dimers (($`y,z`$)-plane) and between different dimer planes (see Fig. 1).

The equivalent form of the Hubbard Hamiltonian (4) is especially suitable because we can apply in the following the exact solution of the dimer problem (5). In the next step we express the construction operators $`a_\sigma `$ ($`a_\sigma ^+`$) as linear combinations of the transition operators between different dimer states. This, in turn, allows to find the exact dimer representation of the whole Hamiltonian (4). The space of the dimer eigenvectors consists, however, of two subspaces. One of them corresponds to the lowest lying energy levels, the second one contains the levels with energies which in the large $`U`$ limit take on large, positive values. These levels in a reasonable temperature range cannot be occupied by electrons and therefore we can exclude them from further considerations. It is interesting to note that this approch, without any additional assumptions, applied to the dimer Hamiltonian (5) produces the analogy of the well known $`tJ`$ model (cf Refs , ). A similar approach can be applied to the Hamiltonian (4), describing the crystal as a whole. By eliminating from the considerations the unoccupied dimer energy levels in the large $`U`$ limit (it is the only assumption) with the use of the projection technique we can find in a straightforward way the final form of the Hamiltonian in this limit without using perturbation expansion or canonical transformation. The resulting Hamiltonian, obtained in this way, is very complicated. It, however, explicitly shows all possible magnetic, nonmagnetic or more complex competitive interaction processes very deeply hidden in the original Hamiltonian written in the site representation (1). It is the aim of this paper just to reveal these important but normally invisible elementary interactions. One important advantage might be that a given approach to the unsolvable Hubbard problem can be tested with respect to the types of neglected interaction processes. Besides, the new and straightforward method, presented in this paper, can easily be adopted to clusters of any size and also to another types of model Hamiltonians.

The paper is organized as follows. In Sec. 2 we find the exact solution of the Hubbard dimer (5) and give the exact expressions for the annihilation operators $`c_{I,1(2),\sigma }`$ in the dimer representation. In Sec. 3 we derive the Hubbard dimer Hamilonian (5) in the large $`U`$ limit. With the use of the projection technique onto the lowest lying dimer states we derive in Sec. 4 the Hubbard Hamiltonian for the crystal in this limit (central formula of this paper). A simplified version of the derived Hamiltonian with the use of the Taylor expansion with respect to $`x=\frac{t}{U}`$ $`(x1)`$ and the case $`U\mathrm{}`$ is discussed in Sec. 5.

## 2 Exact solution of the Hubbard dimer

The eigenvalues and eigenvectors of the dimer Hamiltonian (5) can be found by using a standard procedure (cf Refs \[33-35\]). We start with the vectors $`|n_{1,},n_{1,};n_{2,},n_{2,}`$ $`(n_{i,\sigma }=0,1;`$ $`i=1,2`$; $`\sigma =,)`$ which constitute the Fock basis of the single-dimer space of states:

$$\begin{array}{ccc}\begin{array}{c}\begin{array}{c}|0=|0,0;0,0,\hfill \end{array}\hfill \\ \\ \begin{array}{c}|11=|1,0;0,0,\hfill \\ |12=|0,1;0,0,\hfill \\ |13=|0,0;1,0,\hfill \\ |14=|0,0;0,1,\hfill \end{array}\hfill \end{array}\hfill & \begin{array}{c}|21=|1,1;0,0,\hfill \\ |22=|1,0;1,0,\hfill \\ |23=|1,0;0,1,\hfill \\ |24=|0,1;1,0,\hfill \\ |25=|0,1;0,1,\hfill \\ |26=|0,0;1,1,\hfill \end{array}\hfill & \begin{array}{c}\begin{array}{c}|31=|0,1;1,1,\hfill \\ |32=|1,0;1,1,\hfill \\ |33=|1,1;0,1,\hfill \\ |34=|1,1;1,0,\hfill \end{array}\hfill \\ \\ \begin{array}{c}|4=|1,1;1,1.\hfill \end{array}\hfill \end{array}\hfill \end{array}$$
(6)
Starting with the vectors (6) we easily get the eigenvalues $`E_\alpha `$ and the eigenvectors $`|E_\alpha `$ of the Hubbard dimer (5). Here, we only mention that the space of 16 eigenvectors $`|E_\alpha `$ can be devided into some subspaces numbered by $`n=\underset{i,\sigma }{}n_{i,\sigma }`$. The subspace, belonging to $`n=0`$ and $`n=4`$ is 1-dimensional ($`E_0=0,`$ $`|E_0=|0;`$ $`E_4=2U,`$ $`|E_4=|4)`$. There are, however, two 2-dimensional subspaces corresponding to $`n=1`$ (as e.g. $`E_{11}=t,`$ $`|E_{11}=\frac{1}{\sqrt{2}}(|11+|13),`$ etc.) and two 2-dimensional subspaces corresponding to $`n=3`$ (as e.g. $`E_{31}=t+U,`$ $`|E_{31}=\frac{1}{\sqrt{2}}(|31+|33),`$ etc.). The subspace, belonging to $`n=2`$ consists of one 4-dimensional subspace (as e.g. $`E_{21}=0,`$ $`|E_{21}=\frac{1}{\sqrt{2}}(|23+|24),`$ etc.) and two 1-dimensional subspaces ($`E_{25}=E_{26}=0,`$ $`|E_{25}=|22,`$ $`|E_{26}=|25`$). The complete set of the eigenvalues $`E_\alpha `$ and eigenvectors $`|E_\alpha `$ of the Hubbard dimer is given in the Appendix A. It allows to express the dimer Hamiltonian (5) in the equivalent form

$$H^d=\underset{\alpha }{}E_\alpha |E_\alpha E_\alpha |.$$
(7)
The next important step in our calculations is the possibility to express the annihilation operators $`c_{1(2),\sigma }`$ as linear combinations of transition operators between dimer states $`P_{\alpha ,\beta }=|E_\alpha E_\beta |.`$ This procedure can easily be performed when acting with the annihilation operators on the basis vectors (6) and using the reciprocal relations to (A.1). In this way we can obtain the dimer representation of the annihilation (creation) operators, given in Appedix B. After this operation we can insert so prepared $`c_{I,1(2),\sigma }`$ $`(c_{I,1(2),\sigma }^+)`$ into the Hubbard Hamiltonian (4) to obtain the Hubbard model in the dimer representation for the sc lattice. This representation will be used later to derive the Hubbard Hamiltonian in the large $`U`$ limit (see Sec. 4).

## 3 Hubbard dimer for large $`U`$

Looking at the eigenvalues (A.1) of the Hubbard dimer it is easy to see that for large $`U`$ $`(Ut)`$ the energies $`E_{\alpha \text{ }}=E_{22},E_{23},E_{31},E_{32},E_{33},E_{34}`$ and $`E_4`$ take on large, positive values, producing in the partition function the terms which can practically be neglected. It means that the mentioned energies are not occupied in the reasonable temperature range ($`1`$ $`eV11604.5`$ $`K`$) and can be excluded from our considerations. Therefore the dimer Hamiltonian, given by (7), reduces to

$$\overline{H}^d=tP_{11,11}+tP_{12,12}tP_{13,13}+tP_{14,14}+\left(C+\frac{U}{2}\right)P_{24,24}.$$
(8)
To bring this expression into a compact (second quanization) form we introduce the Hubbard operators

$`a_{i,\sigma }`$ $`=`$ $`c_{i,\sigma }(1n_{i,\sigma }),`$ (9)
$`b_{i,\sigma }`$ $`=`$ $`c_{i,\sigma }n_{i,\sigma }`$ (10)
and spin operators

$`S_i^z`$ $`=`$ $`{\displaystyle \frac{1}{2}}(n_{i,}n_{i,})={\displaystyle \frac{1}{2}}(n_{i,}^an_{i,}^a),`$ (11)
$`S_i^+`$ $`=`$ $`c_{i,}^+c_{i,}=a_{i,}^+a_{i,},`$ (12)
$`S_i^{}`$ $`=`$ $`c_{i,}^+c_{i,}=a_{i,}^+a_{i,}`$ (13)
where $`n_{i,\sigma }^a=a_{i,\sigma }^+a_{i,\sigma }`$ $`(i=1,2;`$ $`\sigma =`$, $`).`$

With the use of (9)-(13) the Hamiltonian of the Hubbard dimer (8) for large $`U`$ can be presented in the form (we introduce the omitted earlier dimer index $`I`$)

$`\overline{H}_I^d`$ $`=`$ $`t{\displaystyle \underset{\sigma }{}}[a_{I,1,\sigma }^+a_{I,2,\sigma }+a_{I,2,\sigma }^+a_{I,1,\sigma }]`$
$`+{\displaystyle \frac{4t^2}{U\sqrt{1+(\frac{4t}{U})^2}}}[\stackrel{}{S}_{I,1}\stackrel{}{S}_{I,2}{\displaystyle \frac{n_{I,1}^an_{I,2}^a}{4}}]`$
$`+{\displaystyle \frac{t^2(1\sqrt{1+(\frac{4t}{U})^2})}{U(1+\sqrt{1+(\frac{4t}{U})^2})\sqrt{1+(\frac{4t}{U})^2}}}[2(b_{I,1,}^+a_{I,1,}^+a_{I,2,}b_{I,2,}`$
$`+b_{I,2,}^+a_{I,2,}^+a_{I,1,}b_{I,1,})`$
$`+(1n_{I,1}^a)n_{I,2}^b+(1n_{I,2}^a)n_{I,1}^bn_{I,1}^bn_{I,2}^b]`$
$`+{\displaystyle \frac{t(1\sqrt{1+(\frac{4t}{U}})^2)}{2\sqrt{1+(\frac{4t}{U})^2}}}{\displaystyle \underset{\sigma }{}}{\displaystyle \underset{\alpha =1,2}{}}[a_{I,\alpha ,\sigma }^+b_{I,\overline{\alpha },\sigma }+b_{I,\alpha ,\sigma }^+a_{I,\overline{\alpha },\sigma }]`$
where $`n_{I,\alpha }^{a,b}=\underset{\sigma }{}n_{I,\alpha ,\sigma }^{a,b},`$ $`n_{I,\alpha ,\sigma }^b=b_{I,\alpha ,\sigma }^+b_{I,\alpha ,\sigma }=n_{I,\alpha ,\sigma }n_{I,\alpha ,\sigma }`$ $`(\alpha =1,2)`$, $`\overline{\alpha }=1`$ when $`\alpha =2`$ and $`\overline{\alpha }=2`$ when $`\alpha =1`$. It is very important to stress that the formula (14) has been obtained in a nonperturbative way, starting from the exact form of the dimer Hamiltonian (see (5) or (7)), and excluding from the considerations unoccupied dimer energy levels in the large $`U`$ limit. Let us note that using the Taylor expansion in (14) with respect to $`x=\frac{t}{U}`$ and retaining the terms proportional to $`\frac{t^2}{U}`$ we obtain

$$\overline{H}_I^d=t\underset{\sigma }{}[a_{I,1,\sigma }^+a_{I,2,\sigma }+a_{I,2,\sigma }^+a_{I,1,\sigma }]+\frac{4t^2}{U}[\stackrel{}{S}_{I,1}\stackrel{}{S}_{I,2}\frac{n_{I,1}^an_{I,2}^a}{4}].$$
(15)
The formula (15) is the well-known $`tJ`$ model for the Hubbard dimer where the first part in (15), similarly to (14), represents the exact form of the dimer Hamiltonian (5) in the limit $`U\mathrm{}`$.

The same result (14) can also be obtained applying a more general approach which will be used later to derive the Hubbard Hamiltonian for large $`U`$ in the case of a crystal. Let us note that after elimination of the unoccupied levels the subspace of the eigenvectors of the Hubbard dimer (A.1) for large $`U`$ consists of the following eigenvectors: $`|E_0`$, $`|E_{11}`$, $`|E_{12}`$, $`|E_{13}`$, $`|E_{14}`$, $`|E_{21}`$, $`|E_{24}`$, $`|E_{25}`$, and $`|E_{26}`$. It means that we can define a projection operator onto this subspace

$$\begin{array}{cc}P_I=\hfill & P_{0,0}^{(I)}+P_{11,11}^{(I)}+P_{12,12}^{(I)}+P_{13,13}^{(I)}+P_{14,14}^{(I)}\hfill \\ & +P_{21,21}^{(I)}+P_{24,24}^{(I)}+P_{25,25}^{(I)}+P_{26,26}^{(I)}\hfill \end{array}$$
(16)
which, in the second quantization form, reads

$$\begin{array}{cc}P_I=\hfill & 1\frac{1}{2}(n_{I,1}^b+n_{I,2}^b)+\frac{1}{4}n_{I,1}^bn_{I,2}^b\hfill \\ & \\ & \begin{array}{cc}& +\frac{1}{8}(1\frac{1}{\sqrt{1+(\frac{4t}{U})^2}})[4(\stackrel{}{S}_{I,1}\stackrel{}{S}_{I,2}\frac{1}{4}n_{I,1}^an_{I,2}^a)\hfill \\ & \\ & +2(b_{I,1,}^+a_{I,1,}^+a_{I,2,}b_{I,2,}+b_{I,2,}^+a_{I,2,}^+a_{I,1,}b_{I,1,})\hfill \\ & \\ & +(1n_{I,1}^a)n_{I,2}^b+(1n_{I,2}^a)n_{I,1}^bn_{I,1}^bn_{I,2}^b]\hfill \end{array}\hfill \\ & \\ & +\frac{t}{U\sqrt{1+(\frac{4t}{U})^2}}\underset{\sigma }{}\underset{\alpha =1,2}{}(a_{I,\alpha ,\sigma }^+b_{I,\overline{\alpha },\sigma }+b_{I,\alpha ,\sigma }^+a_{I,\overline{\alpha },\sigma }).\hfill \end{array}$$
(17)
Now, the Hamiltonian (14) can also be obtained from (5) (or (7)) with the use of the projection operator (16) or (17) using the relation

$$\overline{H}_I^d=P_IH_I^dP_I$$
(18)
and applying a straightforward but laborious algebraic calculation.

## 4 Hubbard model for large $`U`$

The Hamiltonian of the whole crystal (4) can be expressed (similar to (7)) in the form

$$H=\underset{\gamma }{}\overline{E}_\gamma |\overline{E}_\gamma \overline{E}_\gamma |$$
(19)
with unknown energies $`\overline{E}_\gamma `$ and eigenvectors $`|\overline{E}_\gamma `$. We can, however, expand the eigenvectors $`|\overline{E}_\gamma `$ in the series of the dimer eigenvectors (see (A.1))

$$|\overline{E}_\gamma =\underset{\gamma _1,..,\gamma _M}{}c_{\gamma _{1,}.._,\gamma _M}^\gamma |E_{\gamma _1}..|E_{\gamma _M}$$
(20)
assuming that the crystal consists of $`M`$ dimers. Using (19) and (20) we obtain

$$H=\underset{\gamma }{}\overline{E}_\gamma \underset{\gamma _1,..,\gamma _M}{}\underset{\gamma _1^^,,..,\gamma _M^,}{}c_{\gamma _{1,}.._,\gamma _M}^\gamma c_{\gamma _{1,}^,.._,\gamma _M^,}^\gamma |E_{\gamma _1}..|E_{\gamma _M}E_{\gamma _1^,}|..E_{\gamma _M^,}|.$$
(21)
It is clear that to obtain the Hubbard Hamiltonian for large $`U`$ we have to project (21) onto the subspace of the lowest lying dimer states with the use of the projection operator

$$P=P_1P_2\mathrm{}P_M$$
(22)
where $`P_I`$ is given by (16) or (17). In analogy to (14) we denote the Hubbard Hamiltonian in the large $`U`$ limit by $`\overline{H}`$. Similar to (18) we write

$$\overline{H}=PHP$$
(23)
and instead of the form (21) for the Hamiltonian $`H`$ we can use (4). Taking into account that $`P^2=P`$ $`(P_I^2=P_I,`$ $`[P_I,P_J]=0)`$ we obtain

$$\begin{array}{cc}\overline{H}=\hfill & P[\underset{I}{}\overline{H}_I^dt\underset{I,\sigma }{}(\overline{c}_{I,2,\sigma }^+\overline{c}_{I+1,1,\sigma }+\overline{c}_{I+1,1,\sigma }^+\overline{c}_{I,2,\sigma })\hfill \\ & t\underset{IJ,\sigma }{}(\overline{c}_{I,1,\sigma }^+\overline{c}_{J,1\sigma }+\overline{c}_{I,2,\sigma }^+\overline{c}_{J,2,\sigma })]P\overline{\overline{H}}\hfill \end{array}$$
(24)
where $`\overline{H}_I^d`$ is given by (14) and $`(\alpha =1,2)`$

$$\overline{c}_{I,\alpha ,\sigma }=c_{I,\alpha ,\sigma }P_I,$$
(25)
$$\overline{c}_{I,\alpha ,\sigma }^+=P_Ic_{I,\alpha ,\sigma }^+.$$
(26)
Applying the projection operator (16) or (17) to (B.1) - (B.4) and introducing Hubbard- and spin operators (9) - (13) we obtain

$$\overline{c}_{I,1,}=\underset{¯}{\underset{¯}{a}}_{I,1,}+\beta \left[S_{I,2}^za_{I,1,}+S_{I,2}^{}a_{I,1,}\right]\delta \left[S_{I,1}^za_{I,2,}+S_{I,1}^{}a_{I,2,}\right],$$
(27)
$$\overline{c}_{I,1,}=\underset{¯}{\underset{¯}{a}}_{I,1,}\beta \left[S_{I,2}^za_{I,1,}S_{I,2}^+a_{I,1,}\right]+\delta \left[S_{I,1}^za_{I,2,}S_{I,1}^+a_{I,2,}\right]$$
(28)
where the corresponding expressions for $`\overline{c}_{I,2,\sigma }`$ ($`\sigma =,`$) can easily be obtained by changing the internal dimer index $`12`$ in (27) and (28).

The new operators $`\underset{¯}{\underset{¯}{a}}_{I,\alpha ,\sigma }\left(\alpha =1,2;\sigma =,\right)`$ are introduced to obtain a relatively compact form of (27) and (28). They are defined as follows

$$\begin{array}{cc}\underset{¯}{\underset{¯}{a}}_{I,\alpha ,\sigma }=\hfill & \underset{¯}{a}_{I,\alpha ,\sigma }+\beta \left(\underset{¯}{b}_{I,\alpha ,\sigma }+a_{I,\alpha ,\sigma }^+a_{I,\overline{\alpha },\sigma }b_{I,\overline{\alpha },\sigma }\right)\hfill \\ & +\delta \left[\underset{¯}{b}_{I,\overline{\alpha },\sigma }+a_{I,\overline{\alpha },\sigma }^+a_{I,\alpha ,\sigma }b_{I,\alpha ,\sigma }+a_{I,\overline{\alpha },\sigma }\frac{n_{I,\alpha }^a}{2}\right]\hfill \end{array}$$
(29)
where

$$\underset{¯}{a}_{I,\alpha ,\sigma }=a_{I,\alpha ,\sigma }\left(1\frac{\beta }{2}n_{I,\overline{\alpha }}^a\frac{1}{2}n_{I,\overline{\alpha }}^b\right),$$
(30)
$$\underset{¯}{b}_{I,\alpha ,\sigma }=b_{I,\alpha ,\sigma }\left(1n_{I,\overline{\alpha }}^a\frac{1}{2}n_{I,\overline{\alpha }}^b\right)$$
(31)
and

$$\beta =\frac{1}{4}\left(1\frac{1}{\sqrt{1+\left(\frac{4t}{U}\right)^2}}\right),$$
(32)
$$\delta =\frac{t}{U}\frac{1}{\sqrt{1+\left(\frac{4t}{U}\right)^2}}.$$
(33)
To write down the explicit form of the Hubbard Hamiltonian $`\overline{\overline{H}}`$ (see (24)) we have to insert (27) and (28) into $`\overline{\overline{H}}`$. This operation leads, however, to a very complicated form of $`\overline{\overline{H}},`$ given in the Appendix C (a simplified form of this Hamiltonian is discussed in the next Section). Here again the only assumption made to derive $`\overline{\overline{H}}`$ was the reduction of the whole dimer space (A.1) to the subspace of the dimer eigenvectors $`|E_0`$, $`|E_{11}`$, $`|E_{12}`$, $`|E_{13}`$, $`|E_{14}`$, $`|E_{21}`$, $`|E_{24}`$, $`|E_{25}`$ and $`|E_{26}`$, corresponding to the lowest lying dimer energy levels because in the large $`U`$ limit only these levels can be occupied. The Hamiltonian $`\overline{\overline{H}}`$, obtained in this way, contains many competing, magnetic, nonmagnetic and more complex interactions. Among them we can find a direct antiferromagnetic interaction generated by the term $`\underset{I,1}{\overset{}{S}}\underset{I,2}{\overset{}{S}}`$ (Heisenberg exchange interaction) multiplied by the positive coupling constant. Such a term appears in $`\overline{H}_I^d`$ (see (C.6) and (14)). Inside of $`\overline{\overline{H}}`$ (see (C.6) and (C.2) - (C.5)) a kind of ferromagnetic interactions between spins from different dimers, represented by (C.3), appears within the terms proportional to $`\beta ^2`$ and $`\delta ^2`$ (negative coupling constants). The antiferromagnetic interactions, however, appear again in terms proportional to $`\beta \delta `$. There are also many other magnetic, more exotic interactions, represented by (C.2), (C.4) and (C.5), entering into (C.6). The situation is, however, much more complicated when we consider the total Hamiltonian $`\overline{H}`$ (see (24)) in the large $`U`$ limit. $`\overline{H}`$ differs from $`\overline{\overline{H}}`$ by the multiplicative factor $`P`$ (a product of the projection operators $`P_I`$ (see (22) and (17)), standing on the left. Inside of each $`P_I`$ the mentioned antiferromagnetic interaction also appears (see (17)). In other words, the total Hamiltonian $`\overline{H}`$, we are interested in, is actually a sum of the products of many competitive ferromagnetic, antiferromagnetic and more exotic interactions. The thermodynamic properties of the system, described by the Hamiltoniam $`\overline{H}`$ (24), are then a result of the competition between all of them. Which interaction wins in such a competition it certainly depends on temperature, model parameters ($`t,U`$) and on the average number of electrons per lattice site, determining the chemical potential of the system.

The formalism presented in this paper is also applicable to a more complicated decomposition of the Hubbard Hamiltonian (1) into a set od interacting clusters consisting e.g. of one central atom and $`z`$ its nearest neighbours. We, however, know (see e.g. Refs - and papers cited therein) that, unfortunately, the mathematical problems in this case exponentially grows up with the size of the cluster.

## 5 Taylor expansion

The complicated form of the Hubbard Hamiltonian $`\overline{H}`$ in the large $`U`$ limit (see (24), (C.6)) where $`P`$ is given by (22) (see also (17)) can essentially be reduced when applying the Taylor expansion with respect to the parameter $`x=\frac{t}{U}1`$. To do it we have to expand all the coefficients in (14), (17) and (C.6) including also $`\beta `$ and $`\delta `$ (see (32), (33)). Such an expansion can be performed to any power of $`x`$, however, the most simple form we obtain when we restrict ourselves to the linear approximation, resulting in the terms proportional to $`t`$ and $`\frac{t^2}{U}.`$ The accuracy of this expansion can easily be verified assumming e.g. a typical value of the ratio $`\frac{W}{U}=\frac{1}{5}`$ (or less). Because the bandwidth of the conduction band for the sc lattice is $`W=12t,`$ it results in a small value of the expansion parameter $`x=\frac{t}{U}=\frac{1}{60}`$ in this case. It, however, means that the linear approximation is quite reasonable because all higher terms in the expansion, proportional to $`x^n`$ ($`n=2,3,\mathrm{}),`$ produce $`60`$ times smaller contribution.

To present the results of the Taylor expansion including all the terms proportional to $`t`$ and $`\frac{t^2}{U}`$ let us first define several auxiliary quantities

$$P^{(1)}=\underset{I=1}{\overset{M}{}}P_I^{(1)},$$
(34)
$$P^{(2)}=\frac{t}{U}\underset{I=1}{\overset{M}{}}P_1^{(1)}\mathrm{}P_{I1}^{(1)}P_I^{(2)}P_{I+1}^{(1)}\mathrm{}P_M^{(1)}$$
(35)
where

$$P_I^{(1)}=1\frac{1}{2}\left(n_{I,1}^b+n_{I,2}^b\right)+\frac{1}{4}n_{I,1}^bn_{I,2}^b,$$
(36)
$$P_I^{(2)}=\underset{\sigma }{}\underset{\alpha =1}{\overset{2}{}}\left(a_{I,\alpha ,\sigma }^+b_{I,\overline{\alpha },\sigma }+b_{I,\alpha ,\sigma }^+a_{I,\overline{\alpha },\sigma }\right)$$
(37)
and ($`\alpha =1,`$ $`2`$)

$$\stackrel{~}{a}_{I,\alpha ,\sigma }=a_{I,\alpha ,\sigma }\left(1\frac{1}{2}n_{I,\overline{\alpha }}^b\right),$$
(38)
$$\stackrel{~}{\stackrel{~}{a}}_{I,\alpha ,\sigma }=\underset{¯}{b}_{I,\overline{\alpha },\sigma }+a_{I,\overline{\alpha },\sigma }^+a_{I,\alpha ,\sigma }b_{I,\alpha ,\sigma }+a_{I,\overline{\alpha },\sigma }\frac{n_{I,\alpha }^a}{2}$$
(39)
where $`\underset{¯}{b}_{I,\alpha ,\sigma }`$ is given by (31).

The Hamiltonian $`\overline{\overline{H}}`$ (C.6), including all the terms proportional to $`t`$ and $`\frac{t^2}{U}`$, takes on the form

$$\overline{\overline{H}}=\overline{\overline{H}}^{(1)}+\overline{\overline{H}}^{(2)}$$
(40)
where

$`\overline{\overline{H}}^{(1)}`$ $`=`$ $`t{\displaystyle \underset{I,\sigma }{}}\left[a_{I,1,\sigma }^+a_{I,2,\sigma }+a_{I,2,\sigma }^+a_{I,1,\sigma }\right]`$
$`t{\displaystyle \underset{I,\sigma }{}}\left[\stackrel{~}{a}_{I,2,\sigma }^+\stackrel{~}{a}_{I+1,1,\sigma }+\stackrel{~}{a}_{I+1,1,\sigma }^+\stackrel{~}{a}_{I,2,\sigma }\right]`$
$`t{\displaystyle \underset{IJ,\sigma }{}}{\displaystyle \underset{\alpha =1}{\overset{2}{}}}\stackrel{~}{a}_{I,\alpha ,\sigma }^+\stackrel{~}{a}_{J,\alpha ,\sigma }`$
and

$`\overline{\overline{H}}^{(2)}`$ $`=`$ $`{\displaystyle \frac{4t^2}{U}}{\displaystyle \underset{I}{}}[\stackrel{}{S}_{I,1}\stackrel{}{S}_{I,2}{\displaystyle \frac{1}{4}}n_{I,1}^an_{I,2}^a]`$
$`{\displaystyle \frac{t^2}{U}}{\displaystyle \underset{I,\sigma }{}}[\stackrel{~}{a}_{I,2,\sigma }^+\stackrel{~}{\stackrel{~}{a}}_{I+1,1,\sigma }+\stackrel{~}{\stackrel{~}{a}}_{I,2,\sigma }^+\stackrel{~}{a}_{I+1,1,\sigma }`$
$`+\stackrel{~}{a}_{I+1,1,\sigma }^+\stackrel{~}{\stackrel{~}{a}}_{I,2,\sigma }+\stackrel{~}{\stackrel{~}{a}}_{I+1,1,\sigma }^+\stackrel{~}{a}_{I,2,\sigma }]`$
$`{\displaystyle \frac{t^2}{U}}{\displaystyle \underset{IJ,\sigma }{}}{\displaystyle \underset{\alpha =1}{\overset{2}{}}}[\stackrel{~}{a}_{I,\alpha ,\sigma }^+\stackrel{~}{\stackrel{~}{a}}_{J,\alpha ,\sigma }+\stackrel{~}{\stackrel{~}{a}}_{I,\alpha ,\sigma }^+\stackrel{~}{a}_{J,\alpha ,\sigma }]`$
$`+{\displaystyle \frac{2t^2}{U}}{\displaystyle \underset{I}{}}[\stackrel{}{S}_{I,2}(\underset{¯}{\overset{}{\underset{¯}{s}}}_{I,1;I+1,1}+\underset{¯}{\overset{}{s}}_{I+1,1;I,1})`$
$`+\stackrel{}{S}_{I+1,1}(\underset{¯}{\overset{}{\underset{¯}{s}}}_{I+1,2;I,2}+\underset{¯}{\overset{}{s}}_{I,2;I+1,2})]`$
$`+{\displaystyle \frac{2t^2}{U}}{\displaystyle \underset{IJ}{}}[\stackrel{}{S}_{I,1}(\underset{¯}{\overset{}{\underset{¯}{s}}}_{I,2;J,1}+\underset{¯}{\overset{}{s}}_{J,1;I,2})`$
$`+\stackrel{}{S}_{I,2}(\underset{¯}{\overset{}{\underset{¯}{s}}}_{I,1;J,2}+\underset{¯}{\overset{}{s}}_{J,2;I,1})].`$
The operators $`\underset{¯}{\underset{¯}{s}}_{I,\mu ;J,\nu }^{z,\pm }`$ and $`\underset{¯}{s}_{I,\mu ;J,\nu }^{z,\pm }`$ in (42) retain their forms introduced in (C.1) but $`\underset{¯}{\underset{¯}{a}}_{I,\mu ,\sigma }(\underset{¯}{\underset{¯}{a}}_{I,\mu ,\sigma }^+)`$ in (C.1) should actually be replaced by $`\stackrel{~}{a}_{I,\mu ,\sigma }(\stackrel{~}{a}_{I,\mu ,\sigma }^+)`$, defined by (38). The total Hamiltonian $`\overline{H}`$ in the large $`U`$ limit (24) including all terms proportional to $`t`$ and $`\frac{t^2}{U}`$ can thus be written in the form

$$\overline{H}=P^{(1)}(\overline{\overline{H}}^{(1)}+\overline{\overline{H}}^{(2)})+P^{(2)}\overline{\overline{H}}^{(1)}=P^{(1)}\overline{\overline{H}}^{(1)}+(P^{(1)}\overline{\overline{H}}^{(2)}+P^{(2)}\overline{\overline{H}}^{(1)}).$$
(43)
The first part, $`P^{(1)}\overline{\overline{H}}^{(1)}`$, contains the terms proportional to $`t`$ whereas the term in the parentheses is proportional to $`\frac{t^2}{U}`$ (see (35) and (42)). In the expression for $`\overline{\overline{H}}^{(2)}`$(see (42)) the first term describes the antiferromagnetic, Heisenberg intradimer interaction (see also the second term in (15)). The apearence of this term may suggest that such an interaction should also arise between different dimers. This is of course the case. However, because of applied procedure such terms do not appear explicitly. In the last analysis our method treats all interactions within and between dimers, on the same quality level, i.e. all interactions are taken into account. The magnetic interdimer interactions, represented by the fourth and fifth term in (42), have formally the same structure as the Heisenberg interactions but instead of the scalar products of spin operators there are the products of spin operators and ”hopping spin” operators (defined in C.1). All the terms presented in (42) are correct because they originate from the exact decomposition of the Hubbard Hamiltonian (1) into a set of interacting dimers (4) where each dimer problem has been exactly solved (exact dimer representation of the construction operators after applying the projection procedure, given by (27) and (28)). The total Hamiltonian $`\overline{H}`$ (43), we are interested in, is much more complicated than $`\overline{\overline{H}}^{(1)}`$ and $`\overline{\overline{H}}^{(2)}`$alone (see (41)-(43)) because of the presence of the projection operators $`P^{(1)}`$ and $`P^{(2)}`$ (see (34)-(37)) in the expression (43).

It is also interesting to see what happens in the special case when taking the limit $`U\mathrm{}`$. The second term in the parentheses of Eq.(43) vanishes in this limit (cf (35) and (42)). Besides, each lattice site cannot be at the same time occupied by two electrons what is equivalent to the assumption ($`\alpha =1,2)`$

$$n_{I,\alpha ,\sigma }^b=b_{I,\alpha ,\sigma }^+b_{I,\alpha ,\sigma }=n_{I,\alpha ,\sigma }n_{I,\alpha ,\sigma }=0,$$
(44)
$$n_{I,\alpha }^b=\underset{\sigma }{}n_{I,\alpha ,\sigma }^b=0$$
(45)
and (cf (38))

$$\stackrel{~}{a}_{I,\alpha ,\sigma }=a_{I,\alpha ,\sigma }.$$
(46)
The total Hamiltonian (43) is thus given by a simple formula

$$\begin{array}{cc}\overline{H}=P^{(1)}\overline{\overline{H}}=\hfill & t\underset{I,\sigma }{}\left[a_{I,1,\sigma }^+a_{I,2,\sigma }+a_{I,2,\sigma }^+a_{I,1,\sigma }\right]\hfill \\ & \\ & t\underset{I,\sigma }{}\left[a_{I,2,\sigma }^+a_{I+1,1,\sigma }+a_{I+1,1,\sigma }^+a_{I,2,\sigma }\right]\hfill \\ & \\ & t\underset{IJ,\sigma }{}\underset{\alpha =1}{\overset{2}{}}a_{I,\alpha ,\sigma }^+a_{J,\alpha ,\sigma }.\hfill \end{array}$$
(47)
Going back to the original lattice (cf (4),(5) and (1)) it is easy to see that the exact Hubbard Hamiltonian in the limit $`U\mathrm{}`$ takes on the form

$$\overline{H}=t\underset{i,j,\sigma }{}a_{i,\sigma }^+a_{j,\sigma }$$
(48)
where $`i`$ and $`j`$ (as before) number the lattice points and $`a_{i,\sigma }`$ $`(a_{i,\sigma }^+)`$ are the Hubbard operators, defined by (9).

## 6 Conclusions

Using a new, nonperturbative approach, basing on the equivalent form of the Hubbard Hamiltonian, represented by the collection of interacting dimers (4) where each dimer problem has been exactly solved, we have expressed the annihilation (creation) operators (B.1)-(B.4) as linear combinations of transition operators between different dimer states (dimer representation). This method made it possible to exclude from the considerations the unoccupied dimer energy levels in the large $`U`$ limit by means of the projection technique resulting in the final Hamiltonian for the dimer itself (14), as well as, for the crystal as a whole (see (23), (24) and (C.6)). It is important to stress that the elimination of the unoccupied dimer energy levels was the only assumption to derive the Hubbard Hamiltonian (24) in the large $`U`$ limit. Therefore we can be sure that expanding (24) with respect to $`x=\frac{t}{U}`$ (Taylor expansion) we obtain absolutely all terms proportional to $`x^n`$ $`(n=1,2,\mathrm{}).`$ In other words all the coefficients proportional to $`x^n`$ are easy to control what is not always the case when using another methods. The final form of the obtained Hamiltonian in the large $`U`$ limit (see (23), (24) and (C.6)) visualizes high degree of complexity of the model, deeply hidden in its original version, forming a mixture of many multiple ferromagnetic, antiferromagnetic and more complex interactions competing one with another. This fact seems to be decisive for our final conclusion. Because we are still dealing with the approximate solutions of the model (the exact solution does not exist till now) it may happen that we underestimate in this way some important interactions and overestimate the others. It is the reason why the resulting thermodynamic properties of the Hubbard model, obtained in an approximate way, so strongly depend on the quality of applied approximations.

We note in passing that the exact dimer solution may serve as a novel alloy analogy for the Hubbard model, which could be treated by coherent potential approximation (CPA). It is well-known that the standard alloy analogy, based on the atomic limit, does not allow for ferromagnetism in the Hubbard model. This may change by application of the dimer solution which already accounts for a restricted hopping of the band electrons. A corresponding study is in preparation.

Another example of a dimer approach to Hubbard-like models is the bond operator theory as an extension of the slave bosonic and fermionic operators (see Refs ,  and papers cited therein).

Appendix A

The exact solution of the dimer eigenvalue problem (5) reads:

$$\begin{array}{cc}E_0=0;\hfill & |E_0=|0,\hfill \end{array}$$

$$\begin{array}{cc}E_{11}=t;\hfill & |E_{11}=\frac{1}{\sqrt{2}}(|11+|13),\hfill \\ E_{12}=t;\hfill & |E_{12}=\frac{1}{\sqrt{2}}(|11|13),\hfill \\ E_{13}=t;\hfill & |E_{13}=\frac{1}{\sqrt{2}}(|12+|14),\hfill \\ E_{14}=t;\hfill & |E_{14}=\frac{1}{\sqrt{2}}(|12|14),\hfill \end{array}$$

$$\begin{array}{cc}E_{21}=0;\hfill & |E_{21}=\frac{1}{\sqrt{2}}(|23+|24),\hfill \\ E_{22}=U;\hfill & |E_{22}=\frac{1}{\sqrt{2}}(|21|26),\hfill \\ E_{23}=C+\frac{U}{2};\hfill & |E_{23}=a_1(|21+|26)a_2(|23|24),\hfill \\ E_{24}=C+\frac{U}{2};\hfill & |E_{24}=a_2(|21+|26)+a_1(|23|24),\hfill \\ E_{25}=0;\hfill & |E_{25}=|22,\hfill \\ E_{26}=0;\hfill & |E_{26}=|25,\hfill \end{array}$$
(A.1)
$$\begin{array}{cc}E_{31}=t+U;\hfill & |E_{31}=\frac{1}{\sqrt{2}}(|31+|33),\hfill \\ E_{32}=t+U;\hfill & |E_{32}=\frac{1}{\sqrt{2}}(|31|33),\hfill \\ E_{33}=t+U;\hfill & |E_{33}=\frac{1}{\sqrt{2}}(|32+|34),\hfill \\ E_{34}=t+U;\hfill & |E_{34}=\frac{1}{\sqrt{2}}(|32|34),\hfill \end{array}$$

$$\begin{array}{cc}E_4=2U;\hfill & |E_4=|4\hfill \end{array}$$

where

$$C=\sqrt{\left(\frac{U}{2}\right)^2+4t^2},$$
(A.2)
$$a_1=\frac{1}{2}\sqrt{1+\frac{U}{2C}},$$
(A.3)
$$a_2=\frac{1}{2}\sqrt{1\frac{U}{2C}}.$$
(A.4)
Appendix B

The exact dimer representation of the construction operators is given by the following expressions:

$$\begin{array}{cc}c_{1,}=\hfill & \frac{1}{\sqrt{2}}\left(P_{0,11}+P_{0,12}\right)+\frac{1}{\sqrt{2}}P_{11,25}\frac{1}{\sqrt{2}}P_{12,25}\hfill \\ & +\frac{1}{2}\left(P_{13,21}+P_{13,22}\right)+\frac{1}{\sqrt{2}}\left(bP_{13,23}+aP_{13,24}\right)\hfill \\ & \frac{1}{2}\left(P_{14,21}P_{14,22}\right)+\frac{1}{\sqrt{2}}\left(aP_{14,23}bP_{14,24}\right)\hfill \\ & +\frac{1}{2}\left(P_{21,33}P_{21,34}\right)\frac{1}{2}\left(P_{22,33}+P_{22,34}\right)\hfill \\ & +\frac{1}{\sqrt{2}}\left(aP_{23,33}+bP_{23,34}\right)\frac{1}{\sqrt{2}}\left(bP_{24,33}aP_{24,34}\right)\hfill \\ & +\frac{1}{\sqrt{2}}\left(P_{26,31}P_{26,32}\right)+\frac{1}{\sqrt{2}}P_{31,4}+\frac{1}{\sqrt{2}}P_{32,4},\hfill \end{array}$$
(B.1)
$$\begin{array}{cc}c_{1,}=\hfill & \frac{1}{\sqrt{2}}\left(P_{0,13}+P_{0,14}\right)+\frac{1}{\sqrt{2}}P_{13,26}\frac{1}{\sqrt{2}}P_{14,26}\hfill \\ & +\frac{1}{2}\left(P_{11,21}P_{11,22}\right)\frac{1}{\sqrt{2}}\left(bP_{11,23}+aP_{11,24}\right)\hfill \\ & \frac{1}{2}\left(P_{12,21}+P_{12,22}\right)\frac{1}{\sqrt{2}}\left(aP_{12,23}bP_{12,24}\right)\hfill \\ & \frac{1}{2}\left(P_{21,31}P_{21,32}\right)\frac{1}{2}\left(P_{22,31}+P_{22,32}\right)\hfill \\ & +\frac{1}{\sqrt{2}}\left(aP_{23,31}+bP_{23,32}\right)\frac{1}{\sqrt{2}}\left(bP_{24,31}aP_{24,32}\right)\hfill \\ & \frac{1}{\sqrt{2}}\left(P_{25,33}P_{25,34}\right)\frac{1}{\sqrt{2}}P_{33,4}\frac{1}{\sqrt{2}}P_{34,4},\hfill \end{array}$$
(B.2)
$$\begin{array}{cc}c_{2,}=\hfill & \frac{1}{\sqrt{2}}\left(P_{0,11}P_{0,12}\right)\frac{1}{\sqrt{2}}P_{11,25}\frac{1}{\sqrt{2}}P_{12,25}\hfill \\ & \frac{1}{2}\left(P_{13,21}+P_{13,22}\right)+\frac{1}{\sqrt{2}}\left(bP_{13,23}+aP_{13,24}\right)\hfill \\ & \frac{1}{2}\left(P_{14,21}P_{14,22}\right)\frac{1}{\sqrt{2}}\left(aP_{14,23}bP_{14,24}\right)\hfill \\ & \frac{1}{2}\left(P_{21,33}+P_{21,34}\right)+\frac{1}{2}\left(P_{22,33}P_{22,34}\right)\hfill \\ & +\frac{1}{\sqrt{2}}\left(aP_{23,33}bP_{23,34}\right)\frac{1}{\sqrt{2}}\left(bP_{24,33}+aP_{24,34}\right)\hfill \\ & \frac{1}{\sqrt{2}}\left(P_{26,31}+P_{26,32}\right)+\frac{1}{\sqrt{2}}P_{31,4}\frac{1}{\sqrt{2}}P_{32,4},\hfill \end{array}$$
(B.3)
$$\begin{array}{cc}c_{2,}=\hfill & \frac{1}{\sqrt{2}}\left(P_{0,13}P_{0,14}\right)\frac{1}{\sqrt{2}}P_{13,26}\frac{1}{\sqrt{2}}P_{14,26}\hfill \\ & \frac{1}{2}\left(P_{11,21}P_{11,22}\right)\frac{1}{\sqrt{2}}\left(bP_{11,23}+aP_{11,24}\right)\hfill \\ & \frac{1}{2}\left(P_{12,21}+P_{12,22}\right)+\frac{1}{\sqrt{2}}\left(aP_{12,23}bP_{12,24}\right)\hfill \\ & +\frac{1}{2}\left(P_{21,31}+P_{21,32}\right)+\frac{1}{2}\left(P_{22,31}P_{22,32}\right)\hfill \\ & +\frac{1}{\sqrt{2}}\left(aP_{23,31}bP_{23,32}\right)\frac{1}{\sqrt{2}}\left(bP_{24,31}+aP_{24,32}\right)\hfill \\ & +\frac{1}{\sqrt{2}}\left(P_{25,33}+P_{25,34}\right)\frac{1}{\sqrt{2}}P_{33,4}+\frac{1}{\sqrt{2}}P_{34,4}\hfill \end{array}$$
(B.4)
where

$$P_{\alpha ,\beta }=|E_\alpha E_\beta |$$
(B.5)
and

$$\begin{array}{c}a=a_1+a_2\hfill \\ b=a_1a_2.\hfill \end{array}$$
(B.6)
To obtain the annihilation operators of the $`I`$th dimer, the index $`I`$ should be added $`(c_{1(2),\sigma }c_{I,1(2),\sigma }`$, $`P_{\alpha ,\beta }P_{\alpha ,\beta }^{(I)})`$ in (B.1)-(B.4).

Appendix C

To present the Hamiltonian $`\overline{\overline{H}}`$ in a compact form we first introduce the following operators $`\left(\mu ,\nu =1,2\right)`$ :

$$\begin{array}{cc}s_{I,\mu ;J,\nu }^+=\hfill & a_{I,\mu ,}^+a_{J,\nu ,},\hfill \\ & \\ s_{I,\mu ;J,\nu }^{}=\hfill & a_{J,\nu ,}^+a_{I,\mu ,},\hfill \\ & \\ n_{I,\mu ;J,\nu ;\sigma }=\hfill & a_{I,\mu ,\sigma }^+a_{J,\nu ,\sigma },\hfill \\ & \\ n_{I,\mu ;J,\nu }=\hfill & \underset{\sigma }{}n_{I,\mu ;J,\nu ;\sigma },\hfill \\ & \\ s_{I,\mu ;J,\nu }^z=\hfill & \frac{1}{2}\left(n_{I,\mu ;J,\nu ;}n_{I,\mu ;J,\nu ;}\right);\hfill \end{array}$$

$$\begin{array}{cc}\underset{¯}{s}_{I,\mu ;J,\nu }^+=\hfill & \underset{¯}{\underset{¯}{a}}_{I,\mu ,}^+a_{J,\nu ,},\hfill \\ & \\ \underset{¯}{s}_{I,\mu ;J,\nu }^{}=\hfill & a_{J,\nu ,}^+\underset{¯}{\underset{¯}{a}}_{I,\mu ,},\hfill \\ & \\ \underset{¯}{n}_{I,\mu ;J,\nu ;\sigma }=\hfill & \underset{¯}{\underset{¯}{a}}_{I,\mu ,\sigma }^+a_{J,\nu ,\sigma },\hfill \\ & \\ \underset{¯}{n}_{I,\mu ;J,\nu }=\hfill & \underset{\sigma }{}\underset{¯}{n}_{I,\mu ;J,\nu ;\sigma },\hfill \\ & \\ \underset{¯}{s}_{I,\mu ;J,\nu }^z=\hfill & \frac{1}{2}\left(\underset{¯}{n}_{I,\mu ;J,\nu ;}\underset{¯}{n}_{I,\mu ;J,\nu ;}\right);\hfill \end{array}$$
(C.1)
$$\begin{array}{cc}\underset{¯}{\underset{¯}{s}}_{I,\mu ;J,\nu }^+=\hfill & a_{I,\mu ,}^+\underset{¯}{\underset{¯}{a}}_{J,\nu ,},\hfill \\ & \\ \underset{¯}{\underset{¯}{s}}_{I,\mu ;J,\nu }^{}=\hfill & \underset{¯}{\underset{¯}{a}}_{J,\nu ,}^+a_{I,\mu ,},\hfill \\ & \\ \underset{¯}{\underset{¯}{n}}_{I,\mu ;J,\nu ;\sigma }=\hfill & a_{I,\mu ,\sigma }^+\underset{¯}{\underset{¯}{a}}_{J,\nu ,\sigma },\hfill \\ & \\ \underset{¯}{\underset{¯}{n}}_{I,\mu ;J,\nu }=\hfill & \underset{\sigma }{}\underset{¯}{\underset{¯}{n}}_{I,\mu ;J,\nu ;\sigma },\hfill \\ & \\ \underset{¯}{\underset{¯}{s}}_{I,\mu ;J,\nu }^z=\hfill & \frac{1}{2}\left(\underset{¯}{\underset{¯}{n}}_{I,\mu ;J,\nu ;}\underset{¯}{\underset{¯}{n}}_{I,\mu ;J,\nu ;}\right).\hfill \end{array}$$

The operators $`s_{I,\mu ;J,\nu }^{z,\pm }`$, $`\underset{¯}{s}_{I,\mu ;J,\nu }^{z,\pm }`$, $`\underset{¯}{\underset{¯}{s}}_{I,\mu ;J,\nu }^{z,\pm }`$, are not strictly the spin operators, they, however, show some similarities to the true spin operators (as e.g. (11)-(13)) and therefore they can be called ”hopping spin operators”. Moreover, the following abreviations have to be used:

$$Q(I,\mu ;J,\nu )=\stackrel{}{S}_{I,\mu }(\underset{¯}{\overset{}{\underset{¯}{s}}}_{I,\overline{\mu };J,\nu }+\underset{J,\nu ;I,\overline{\mu }}{\overset{}{\underset{¯}{s}}}),$$
(C.2)
$$R(I,\mu ;J,\nu )=\underset{I,\mu }{\overset{}{S}}\underset{J,\nu }{\overset{}{S}}n_{I,\overline{\mu };J,\overline{\nu }},$$
(C.3)
$$R^{z,\pm }(I,\mu ;J,\nu )=\left(S_{I,\mu }^zS_{J,\nu }^\pm S_{I,\mu }^\pm S_{J,\nu }^z\right)s_{I,\overline{\mu };J,\overline{\nu }}^{},$$
(C.4)
$$R^{,+}(I,\mu ;J,\nu )=\left(S_{I,\mu }^{}S_{J,\nu }^+S_{I,\mu }^+S_{J,\nu }^{}\right)s_{I,\overline{\mu };J,\overline{\nu }}^z.$$
(C.5)
The Hamiltonian $`\overline{\overline{H}}`$ (cf (24)) with the use of (C.1) - (C.5) takes on the form

$$\begin{array}{cc}\overline{\overline{H}}=\hfill & \underset{I}{}\overline{H}_I^dt\underset{I,\sigma }{}\left[\underset{¯}{\underset{¯}{a}}_{I,2,\sigma }^+\underset{¯}{\underset{¯}{a}}_{I+1,1,\sigma }+\underset{¯}{\underset{¯}{a}}_{I+1,1,\sigma }^+\underset{¯}{\underset{¯}{a}}_{I+1,2,\sigma }\right]\hfill \\ & \\ & 2t\beta \underset{I}{}\left[Q(I,1;I+1,1)+Q(I+1,2;I,2)\right]\hfill \\ & \\ & +2t\delta \underset{I}{}\left[Q(I,2;I+1,1)+Q(I+1,1;I,2)\right]\hfill \end{array}$$

$$\begin{array}{cc}t\beta ^2\underset{I}{}\hfill & [R(I,1;I+1,2)+R^{z,}(I,1;I+1,2)\hfill \\ & \\ & +R^{z,+}(I,1;I+1,2)+R^{,+}(I,1;I+1,2)\hfill \\ & \\ & +R(I+1,2;I,1)+R^{z,}(I+1,2;I,1)\hfill \\ & \\ & +R^{z,+}(I+1,2;I,1)+R^{,+}(I+1,2;I,1)]\hfill \end{array}$$

$$\begin{array}{cc}t\delta ^2\underset{I}{}\hfill & [R(I,2;I+1,1)+R^{z,}(I,2;I+1,1)\hfill \\ & \\ & +R^{z,+}(I,2;I+1,1)+R^{,+}(I,2;I+1,1)\hfill \\ & \\ & +R(I+1,1;I,2)+R^{z,}(I+1,1;I,2)\hfill \\ & \\ & +R^{z,+}(I+1,1;I,2)+R^{,+}(I+1,1;I,2)]\hfill \end{array}$$

$$\begin{array}{cc}+t\beta \delta \underset{I}{}\underset{\mu =1}{\overset{2}{}}\hfill & [R(I,\mu ;I+1,\mu )+R^{z,}(I,\mu ;I+1,\mu )\hfill \\ & \\ & +R^{z,+}(I,\mu ;I+1,\mu )+R^{,+}(I,\mu ;I+1,\mu )\hfill \\ & \\ & +R(I+1,\mu ;I,\mu )+R^{z,}(I+1,\mu ;I,\mu )\hfill \\ & \\ & +R^{z,+}(I+1,\mu ;I,\mu )+R^{,+}(I+1,\mu ;I,\mu )]\hfill \end{array}$$
(C.6)
$$\begin{array}{cc}& t\underset{IJ,\sigma }{}\underset{\mu =1}{\overset{2}{}}\underset{¯}{\underset{¯}{a}}_{I,\mu ,\sigma }^+\underset{¯}{\underset{¯}{a}}_{J,\mu ,\sigma }\hfill \\ & \\ & 2t\beta \underset{IJ}{}\underset{\mu =1}{\overset{2}{}}Q(I,\mu ;J,\overline{\mu })+2t\delta \underset{IJ}{}\underset{\mu =1}{\overset{2}{}}Q(I,\mu ;J,\mu )\hfill \end{array}$$

$$\begin{array}{cc}& t(\beta ^2+\delta ^2)\underset{IJ}{}\underset{\mu =1}{\overset{2}{}}[R(I,\mu ;J,\mu )+R^{z,}(I,\mu ;J,\mu )\hfill \\ & \\ & +R^{z,+}(I,\mu ;J,\mu )+R^{,+}(I,\mu ;J,\mu )]\hfill \\ & \\ & +2t\beta \delta \underset{IJ}{}\underset{\mu =1}{\overset{2}{}}[R(I,\mu ;J,\overline{\mu })+R^{z,}(I,\mu ;J,\overline{\mu })\hfill \\ & \\ & +R^{z,+}(I,\mu ;J,\overline{\mu })+R^{,+}(I,\mu ;J,\overline{\mu })]\hfill \end{array}$$

where $`\overline{H}_I^d`$ is given by (14). It should also be noted that the operators $`\underset{¯}{\underset{¯}{a}}_{I,\alpha ,\sigma }`$ $`(\underset{¯}{\underset{¯}{a}}_{I,\alpha ,\sigma }^+),`$ entering in (C.1) and (C.2), are in fact, dependent on $`\beta `$ and $`\delta `$ (see (29)). The decomposition of the Hamiltonian $`\overline{\overline{H}}`$ (C.6) according to the terms proportional to $`\beta ,`$ $`\delta `$, $`\beta ^2`$,$`\delta ^2`$ and $`\beta \delta `$ has thus only a formal character in order to keep the presentation of the Hamiltonian $`\overline{\overline{H}}`$ in a compact form.

The formula (C.6) is the complete expression.


format: text

text:
# QUASI-SADDLES OF LIQUIDS: COMPUTATIONAL STUDY OF A BULK LENNARD-JONES SYSTEM

## 1 Introduction

A number of theoretical approaches have been developed to understand the relationship between the potential energy surface (PES) and the properties of condensed phases. Since the potential energy surface, $`U(𝐫)`$, is a complicated function of the $`3N`$ dimensional position coordinates of an $`N`$-atom system, most approaches have concentrated on identifying features of $`U(𝐫)`$ most relevant to understanding a particular phenomenon. The most widely used of such energy landscape approaches is inherent structure analysis which focuses on the sampling of local minima of the PES . To perform the inherent structure sampling, configurations are sampled from a well-defined ensemble and are referred to as instantaneous configurations in this work. A steepest-descent or equivalent local minimisation procedure is then initiated for each instantaneous configuration in order to obtain the nearest local minimum. All instantaneous configurations which quench to the same minimum or inherent structure are referred to as belonging to the basin of that minimum. In general, the phase with lowest entropy and internal energy will correspond to location of the system in the basin of the global minima and phase transitions will correspond to occupancy of higher energy metastable minima. Inherent structure analysis has proved to be an extremely useful tool in understanding melting, the glass transition and protein folding . It is, however, limited by the fact that it focuses entirely on minima whereas an understanding of dynamical phenomena requires a knowledge of factors which govern the passage of the system between basins of different minima. For example, the dynamical slowing down associated with the glass transition requires an understanding of why rate constants for inter-basin transitions become essentially zero on experimental time scales. In this context, it has recently been suggested that saddle configurations of the PES, which mark the border between adjoining basins of minima, must play a crucial role. Such saddle configurations will have $`U=0`$ and will show negative curvature in one or more orthogonal directions. Saddles of order 1 with just a single negative curvature direction will correspond to a simple transition state connecting two adjacent minima. These stationary points of $`U(𝐫)`$ will also correspond to absolute minima of the pseudo-potential surface, $`W(𝐫)=\frac{1}{2}|U|^2`$ . To identify the saddle points which play a significant role in the system dynamics, it has been suggested that inherent saddles, analogous to inherent structures, can be identified by a steepest descent minimisation from an instantaneous configuration to the nearest local minimum on the pseudo-potential surface, W(r)= $`\frac{1}{2}|U|^2`$, where $`U(𝐫)`$ is the potential energy as a function of the 3N-dimensional position vector $`𝐫`$. The statistics of inherent saddles has been recently analysed for several model glass-formers, including a binary Lennard-Jones mixture(BLJ), a modified monoatomic Lennard-Jones(MLJ) liquid and the binary soft-sphere (BSS) mixture. Simulations show that these systems undergo a transition from saddle-dominated to minima-dominated dynamical regimes at the mode-coupling transition temperature, $`T_{MC}`$. For $`T>T_{MC}`$, the system is localised largely in border regions and basin-hopping is facile. Below $`T_{MC}`$, the average saddle order goes to zero implying that the system occupies interiors of the basins of local minima and basin-hopping becomes an activated process. The saddle order is found to be an essentially linear function of the saddle configurational energy. By extrapolating the linear trend one can obtain a threshold energy, $`U_{thr}`$, at which the saddle order goes to zero. The average configurational energy at the mode-coupling temperature is found to be very close to $`U_{thr}`$. Since a glass may be viewed as a very high-viscosity liquid, it is expected that many of the statistical features of the stationary points of the PES of glass-formers will be shared by simple liquids. Recent studies of the stationary points of Lennard-Jones and Morse liquids verify this conjecture . The statistical properties of inherent saddles are shown to carry interesting signatures steming from variations in generic features of the interatomic potential, such as the range and curvature; for example, the configurational energy is shown to depend linearly on saddle order with a slope that is proportional to the range of the pair potential. These recent studies therefore provide strong evidence that stationary points of the potential energy surface constitute a physically significant set of topographic features of the energy landscape.

The computational strategy for identifying stationary points or inherent saddles in the above studies has been to locate the minima of the pseudo-potential surface, $`W(𝐫)=\frac{1}{2}|U|^2`$; within the limits of computational precision, the minima for which $`W(𝐫)=0`$ are taken to be stationary points . As pointed out in ref., this approach must, however, be applied very carefully because any minimisation algorithm will sample not only absolute minima of $`W(𝐫)`$ which correspond to the stationary points of the true potential $`U(𝐫)`$ but also low-lying minima of $`W(𝐫)`$ which correspond to inflexion points of $`U(𝐫)`$ . At an inflexion point, the Hessian matrix, $`H`$, must have one or more zero eigenvalues; the eigenvectors corresponding to these zero eigenvalues do not correspond to basin-crossing displacements. The question therefore arises as to whether the physical significance of inherent saddle analysis is compromised because, given finite numerical precision limits, it may be difficult to distinguish between true saddles and inflexion points. The available numerical studies of this issue are limited but they indicate that the number of inflexion directions is always very small and does not significantly affect the statistical estimates of saddle properties . For example, in the case of the 256 atom MLJ system studied in ref.14, no point with more that four inflexion directions was found. In the liquid regime, the ratio of inflexion directions to negative curvature directions was approximately 0.1. It rises to 0.8 very close to $`T_{MC}`$ but even then no major shift in the estimated value of $`T_{MC}`$ was seen.

Given the usefulness of stationary point analysis in understanding liquid-state dynamics, we feel it is worthwhile to perform a detailed computational study of a model system. In this work, we therefore investigate the sensitivity of statistical properties of inherent saddles to details of the minimisation procedure. For this purpose, we define a simple liquid bound by a quadratically shifted Lennard-Jones pair potential with continuous zeroth, first and second derivatives at the cut-off distance. This ensures that no loss of precision in the minimisation procedure is present due to discontinuities in the pair potential at a finite cut-off distance. Previous studies of stationary points of bulk liquids and glasses have used either unshifted pair potentials with long-range corrections  or quadratically shifted pair potentials with continuous zeroth and first derivatives but discontinuous second derivatives . The only exception in the literature known to us is the simulation of the binary soft-sphere system but this may represent a special case because of the absence of any attractive tail in the pair potential . Unlike in other studies, we consider a one-component bulk liquid, rather than a binary mixture; these two classes of systems differ only in the ease with which crystal nucleation takes place but the crucial statistical features of stationary points are common to both glass-formers and simple liquids. This monoatomic liquid bound by a quadratically shifted Lennard-Jones potential (QSLJ) defines a test system for which we can set reasonably demanding accuracy limits for the minimisation procedure. The low-lying minima of the pseudo-potential, $`W(𝐫)=\frac{1}{2}|U|^2`$, as located by the minimisation procedure are referred to as quasi-saddles; however, unless the context requires distinguishing clearly between true saddles and inflexion points, we also refer to quasi-saddles as saddle configurations. Changes in various properties of the saddle configurations are analysed as a function of the precision limits set for the minimisation algorithm. The possibility of setting up an unambiguous numerical criterion to distinguish between true saddles and inflexion points is also discussed This allows us to obtain a useful benchmark of the degree of error introduced by a given level of convergence in the minimisation procedure. The results should be useful in two contexts. Firstly, it allows one to estimate the computational effort to be expended in minimisation given the saddle properties of interest. Secondly, in a number of contexts, it is of interest to study liquids which do not have a finite range e.g. the Lennard-Jones and other well-known model systems. If the statistical analysis of saddle configurations is to be a generally applicable tool to study the disordered phase (liquids and glasses), it is essential to establish to what extent such an analysis will be limited by computational problems associated with the minimisation procedure. The paper is organised as follows. Section 2 contains details regarding the computational methods employed. Results are presented in Section 3. Section 4 contains the conclusions.

## 2 Computational Methods

### 2.1 Potential Energy Surface

The potential energy surface (PES), $`U(𝐫)`$, is defined as a scalar function of the 3N-dimensional position vector, r. The PES studied in this work is pair-additive i.e.

$$U(𝐫)=\underset{i=1}{\overset{N}{}}\underset{j<i}{}u_p(r_{ij})$$
(1)
where $`r_{ij}`$ is the distance between atoms $`i`$ and $`j`$ and the pair potential, $`u_p(r)`$, is a quadratically shifted Lennard-Jones function of the form

$`u_p(r)`$ $`=`$ $`4ϵ\left\{\left({\displaystyle \frac{\sigma }{r}}\right)^{12}\left({\displaystyle \frac{\sigma }{r}}\right)^6\right\}+ar^2+br+c\mathrm{for}rR_c`$ (2)
$`=`$ $`0\mathrm{for}r>R_c.`$
The $`a`$, $`b`$ and $`c`$ parameters are determined by requiring that $`u_p`$, $`du_p/dr`$ and $`d^2u_p/dr^2`$ be zero at $`r=R_c`$. We take $`ϵ`$, $`\sigma `$ and $`m`$ as reduced units for energy, length and mass respectively; all our results are given in reduced units. The pair potential defined in equation (2) is referred to as the QSLJ potential in this paper. The two-parameter quadratically shifted Lennard-Jones potential used in previous studies of saddles of binary LJ mixtures corresponds to a shifting function $`ar^2+c`$ and is referred to as the SLJ potential. The monoatomic SLJ system was originally discussed in ref.. We fix the spherical cut-off distance at 2.5$`\sigma `$ when computing the shifting parameters. Figure 1 compares the QSLJ, SLJ and unshifted Lennard-Jones (LJ) potentials.

### 2.2 NPT Ensemble Monte Carlo Simulations

All the simulations reported in this paper were performed using the isothermal-isobaric ensemble Monte Carlo method of McDonald . The partition function of an $`N`$-particle system in this ensemble may be written as:

$$Q_{NPT}=\frac{\beta P}{\mathrm{\Lambda }^{3N}N!}d(\mathrm{ln}v)v^{N+1}\mathrm{exp}(\beta Pv)𝑑𝐬\mathrm{exp}(\beta V(𝐬;L))$$
(3)
where $`P`$ is the external pressure, $`\beta =1/k_BT`$ is the inverse temperature, $`\mathrm{\Lambda }=\sqrt{h^2/2\pi mk_BT}`$ is the thermal de Broglie wavelength, $`v`$ and $`L`$ are the volume and length of the simulation cell respectively and $`𝐬=𝐫/L`$ is the 3N-component position vector scaled with respect to the length of the simulation cell. The volume, $`v`$, and the particle positions, $`𝐫`$, are the Metropolis variables. The random walk is performed using two types of moves: (i) volume moves which sample in $`\mathrm{ln}v`$ and (ii) particle moves which displace the spatial coordinates of a single, randomly chosen particle. In the case of the unshifted Lennard-Jones system, long-range corrections were employed and the spherical cut-off distance, $`R_c`$, was scaled whenever a volume move was made; the average value of $`R_c2.8\sigma `$. For the shifted potentials, the cut-off distance was fixed at 2.5$`\sigma `$.

A simulation cell with 125 particles and rhombic dodecahedral boundary conditions was used in the Monte Carlo simulations. The solid phase structure taken to be a face-centred cubic (fcc) lattice. Volume moves constituted 7.4% of trial moves. Acceptance ratios for volume and particle moves were kept at approximately 50% . Run lengths were between two and four million configurations with equilibriation periods of one million to two million. This was found to be sufficient to converge configurational averages, such as the average potential energy, density and bond orientational parameters, to better than 1%. In order to obtain statistical averages over inherent saddle configurations, 100 instantaneous configurations were sampled from the NPT-MC simulations at equal intervals and used to generate inherent saddles. The simulations were carried out at a reduced pressure of $`P^{}=0.67`$. At this pressure, the melting point of the unshifted Lennard-Jones system has been established from previous work as $`T^{}=0.75`$ .

### 2.3 Locating Inherent Saddles

Locating inherent saddles requires one to find the absolute minima for the pseudo-potential function, $`W(𝐫)=\frac{1}{2}|U|^2`$. Any local minimisation technique equivalent to a steepest descent procedure can be used to perform the minimisation. In a previous study , three minimisation procedures for locating inherent saddles were compared: (i) the conjugate gradient algorithm given in ref.. (ii) the variable metric LBFGS algorithm, as implemented by Liu and Nocedal  and (iii) Powell’s direction set method which requires no gradient information. The LBFGS algorithm was found to be computationally the most efficient and was therefore used in this study to locate local minima in $`W(𝐫)`$ using $`W=𝐇𝐅`$. The termination condition for the LBFGS algorithm is

$$g_i\tau \mathrm{max}(1,x_i)$$
(4)
where $`𝐠_i`$ and $`𝐱_i`$ are the gradient and solution vectors for the $`i`$-th iteration and $`||||`$ denotes the Euclidean norm. Clearly the more demanding the accuracy requirements i.e. the smaller the error tolerance $`\tau `$, the more computationally demanding the minimisation procedure will be. This is illustrated in Table 1 which shows the average number of iterations required to converge to the nearest saddle for a given value of $`\tau `$ for three different temperatures. The error tolerance $`\tau `$ is the key convergence parameter with respect to which we wish to test the sensitivity of properties of quasi-saddles.

While absolute minima of $`W(𝐫)`$ must correspond to true stationary points of $`U(𝐫)`$, local minima in $`W(𝐫)`$ need not correspond to true stationary points of $`U(𝐫)`$ with all elements of the gradient vector, $`U`$, equal to zero. Some local minima may be inflexion points where both the first and second derivatives in some directions are zero. In principle, such flex points can be distinguished from true saddles by looking at the absolute magnitude of $`W(𝐫)`$. However, given finite numerical precision it may be difficult to unequivocally distinguish between low-lying and absolute minima. Alternatively, such inflexion points can be located by examining whether any of the eigenvalues of the Hessian for a configuration corresponding to a minimum in $`W(𝐫)`$, other than the three translational modes, is unusually small. However, the lowest magnitude, non-translational, eigenvalue of the Hessian is very sensitive to the accuracy of the minimisation procedure. Relaxed tolerances can lead to dramatic undercounting of the number of flex points while rigorous convergence criteria lead to almost 99% of local minima of $`W(𝐫)`$ being identified as inflexion points (for example, one can contrast the results presented in refs. and .

## 3 Results and Discussion

### 3.1 Isobars of the Quadratically Shifted Lennard-Jones liquid

The behaviour of the qudratically shifted Lennard-Jones system (QSLJ), as defined in Section 2.1, is mapped out along the $`P^{}=0.67`$ isobar to facilitate comparison with the unshifted Lennard-Jones system. We have also performed NPT ensemble simulations at the same pressure for the two-parameter quadratically shifted Lennard-Jones (SLJ) potential with discontinuous second derivative, as defined by Stoddard and Ford. The binary Lennard-Jones mixtures for which inherent saddle analysis has been performed recently employ the SLJ type of cut-off conditions. Figures 2(a) and 2(b) compare the temperature dependence of the configurational energy per particle and the number density respectively along the $`P^{}=0.67`$ isobar for the QSLJ, SLJ and LJ systems. The thermodynamic melting temperature, $`T_m^{}`$, is known to be 0.75 for the LJ system at this pressure . The metastability limit of the solid from our NPT simulations is $`T_k^{}=0.88`$; such a 15-20% difference between $`T_m^{}`$ and $`T_k^{}`$ is fairly typical. Though the difference in the well depth of the pair potential between the LJ and SLJ systems is very small (see Figure 1), the $`T_k^{}`$ values for the SLJ system is about 13% lower than that of the unshifted LJ system because of the finite-range of the SLJ pair potential. The QSLJ system has an even lower value of $`T_k^{}=0.64`$ because of the reduced well-depth. At any given temperature, the QSLJ system has the lowest number density and highest configurational potential energy. Since we are interested in the stationary points sampled in the disordered liquid phase, the results of the inherent saddle analysis discussed below refer to liquid state simulations for $`T^{}0.75`$.

### 3.2 Sensitivity of Saddle Properties to Accuracy of Minimisation Procedure

A saddle point is characterised by its configurational energy, $`U_s`$, and its index density, $`n_s`$, which corresponds to the fraction of imaginary modes; the NPT averages are denoted by $`U_s`$ and $`n_s`$. In the work of Broderix et al, $`U_s`$ and $`n_s`$ are referred to as parametric averages with respect to $`T`$ as the parameter. For the binary LJ mixture, they showed that the points $`(U_s,n_s)`$ lie on a straight line. They also defined the geometric average, $`U_s(n_s)`$ as the average configurational energy of all saddles of index density $`n_s`$ and showed that the straight line fits to the geometric and parametric averages are identical, within the limits of statistical error. This linear relationship between the saddle order and configurational energy has been subsequently shown to be true for Lennard-Jones and Morse liquids; moreover, the slope, $`U_s/n_s`$ has been found to be proportional to the range of the potential . The linear form of the $`U_s(n_s)`$ curves in the disordered, liquid phase can be understood on the basis of a localised mechanism for generating imaginary frequency displacement modes which requires rearrangement of only a few neighbouring atoms.

Figure 3(a) shows the parameteric points $`(U_s,n_s)`$ obtained from minimisations corresponding to three different values of the error tolerance, $`\tau =10^{10}`$, $`10^{12}`$ and $`10^{15}`$. It is well know that error tolerances for systems with finite cut-offs and discontinuous seroth, first or seond derivatives cannot be made very small. The tolerance value of $`ϵ=10^{10}`$ corresponds to that used in a previous study of the unshifted LJ liquid . The tolerances used for the Morse potentials in a more recent study are similar. The slopes of the straight line fits are almost unchanged over five orders of magnitude in $`\tau `$. The results at $`\tau =10^{12}`$ and $`10^{15}`$ are indistinguishable. For $`\tau =10^{10}`$, for a given value of the saddle energy, the saddle order is somewhat overestimated. Figure 3(b) shows the temperature dependence of the ensemble averaged index density. The behaviour with $`\tau `$ is similar to that seen for the $`U_s(n_s)`$ curves. Figure 3(c) shows the ensemble-averaged saddle configurational energies as a function of temperature. The convergence with $`\tau `$ is much better for the saddle energies than for index densities. The reason is that the convergence of the near-zero eigenvalues of the Hessian is relatively strongly affected by the error tolerance compared to the overall saddle configurational energies. To substantiate this point, we have compared the change in the absolute magnitude of the non-translational eigenvalues, $`|\omega _i^2|`$, of the Hessian with $`\tau `$ for several configurations. Figure 4 shows the values of the lowest, $`|\omega _0^2|`$, and third lowest, $`|\omega _2^2|`$, magnitude eigenvalues of the Hessian for an arbitrary sequence of configurations for $`\tau =10^{14}`$ and $`10^{15}`$. The lowest magnitude eigenvalue varies by over a factor of 10 when $`\tau `$ is reduced by an order of magnitude. There is, on average, a decrease in the magnitude of $`\omega _0^2`$ as $`\tau `$ decreases. The magnitudes of the third-lowest eigenvalue (and all other higher magnitude eigenvalues) are practically unchanged with $`\tau `$. The sensitivity of the lowest magnitude eiegenvalues of the Hessian to the error tolerance of minimisation is responsible for the relatively large error in $`n_s`$ compared to $`U_s`$ discussed above.

### 3.3 Distinguishing between True Saddles and Inflexion Points

We have approached the possibility of setting up an unambiguous numerical criterion to separate the quasi-saddles into true saddles and inflexion points by comparing the location of each quasi-saddle on a two-dimensional plot with $`x`$ and $`y`$ axes corresponding to $`W(𝐫_s)`$ and $`|\omega _0^2|/\omega _E^2`$ respectively where $`\omega _0^2`$ is the lowest magnitude non-translational eigenvalue of the Hessian (see Figure 5) and $`\omega _E`$ is the Einstein frequency. True saddles should correspond to very small values of $`W`$ and relatively large values of $`|\omega _0^2|`$ and therefore should cluster in the upper left hand corner of the graph. Inflexion points should correspond to relatively larger values of $`W`$ and near-zero values of $`|\omega _0^2|`$ and therefore should cluster in the lower right hand corner of the graph. If two such well-defined clusters of points are seen, then an unambiguous numerical distinction between true saddles and inflexion points can be made. Figure 5 shows, however, that the situation is quite the opposite. Results for two different temperatures are shown and are representative of the behaviour in other simulation runs. Only a few points (less than 5 in 100) can be classed unambiguously as true saddlles; the others form an almost continuous distribution. This is presumably why in ref.16 the use of tight convergence criteria for minimisation led to an identification of 99% of the quasi-saddles as inflexion points. We suggest that the reason for this behaviour is that in the absence of symmetry constraints requiring some collective normal mode directions of quasi-saddles to have zero curvature, the distinction between saddles and inflexion points is not physically very meaningful because it attempts to distinguish between zero and near-zero curvature situations. Coupled with the fact that the number of near-zero curvature modes is always small, as can be seen even from the saddle normal mode spectra described in ref. as well as from the results in ref., it would appear that the effect of reducing $`\tau `$ will be to introduce small corrections in the saddle index density. Since the physical significance of such quasi-saddles stems from the fact that the majority of negative curvature directions will correspond to basin-crossing displacements, improving the error tolerance will not qualitatively alter any results.

## 4 Conclusions

The purpose of this study was to test the sensitivity of statistical properties of quasi-saddles of the potential energy surface, $`U(𝐫)`$, to the convergence criteria used for minimisation on the pseudo-potential surface, $`W(𝐫)=|U|^2`$. For this purpose a quadratically shifted Lennard-Jones (QSLJ) pair potential was defined which has continuous zeroth, first and second derivatives at the cut-off distance. Inherent saddles (or quasi-saddles) of the bulk QSLJ liquid were studied as a function of the error tolerance, $`\tau `$, of the LBFGS minimisation algorithm. The variation in statistical properties of saddles was studied over a range of $`\tau `$ values covering five orders of magnitude. The largest $`\tau `$ values correspond to those used for unshifted LJ and Morse liquids simulations previously. Based on our results, it is clear that there are no qualitative changes in statistical proeprties of saddles over this range of error tolerance and even the quantitative changes are small. The ensemble-averaged saddle index densities were found to be more sensitive to $`\tau `$ than the corresponding averages over saddle configurational energies. The lowest magnitude eigenvalue, $`|\omega _0|^2`$, of the Hessian was, however, found to be very sensitive to the tolerance; as $`\tau `$ was decreased, $`|\omega _0|^2`$, was found to show an overall decrease. This indicates that if convergence criteria are not strict, absolute or low-lying minima of $`W(𝐫)`$ willl be diagnosed as having no inflexion directions, as was found in ref.12 and 18. For the smallest value of $`\tau `$, we also considered the possibility of setting up a numerical criterion to clearly distinguish between true saddles and inflexion points by considering the value of both $`W(𝐫)`$ as well as $`|\omega _0^2|`$ at a quasi-saddle identified by minimisation on the pseudo-potential surface. We find however, that the quasi-saddles do not fall into two such clear and distinct categories. The reason appears to be that of the $`3N`$ normal mode directions associated with a low-lying or absolute minima of $`W(𝐫)`$, a few directions will correspond to near-zero curvatures. Deducing which of these near-zero curvatures is actually identically equal to zero and therefore corresponds to an inflexion direction is necessary in order to distinguish between true saddles and inflexion points Moreover, this appears not to be a physically very meaningful distinction since the physically important feature of quasi-saddles is the existence of negative curvature directions indicating their location at the border between adjacent minima.

Acknowledgements CC would like to thank the Council for Scientific and Industrial Research for financial support and the Computer Services Centre of I.I.T-Delhi for access to their computational resources. PS would like to thank the Council for Scientific and Industrial Research, New Delhi for the award of a Senior Research Fellowship.

Table 1

Average Number of iterations of the LBFGS minimisation algorithm for a given error tolerance, $`\tau `$, at a reduced temperature $`T^{}`$ for the quadratically shifted Lennard-Jones (QSLJ) liquid.

| $`T^{}`$ | $`\tau =10^{10}`$ | $`\tau =10^{12}`$ | $`\tau =10^{15}`$ |
| --- | --- | --- | --- |
| 0.67 | 66 | 1306 | 11469 |
| 0.88 | 54 | 1077 | 10895 |
| 1.29 | 35 | 898 | 17362 |

Figure Captions

1. Comparison of the quadratically shifted Lennard-Jones (QSLJ) pair potential with continuous zeroth, first and second derivatives with the two-parameter quadratically shifted Lennard-Jones (SLJ) potential defined in ref.22 as well with the unshifted Lennard-Jones (LJ) potential. The well-depth, $`ϵ`$, and size parameter, $`\sigma `$, of the LJ potential are taken as reduced units of energy and length respectively.
2. Dependence of configurational properties on reduced temperature $`T^{}`$ at a reduced pressure $`P^{}`$=0.67 for the bulk LJ, SLJ and QSLJ systems: (a) average potential energy per particle $`<U^{}>`$ and (b) reduced density $`\rho ^{}`$
3. Statistical properties of saddle configurations for three different values of the error tolerance, $`\tau `$, used in the LBFGS minimisation algorithm: (a) Parametric averages of saddle configurational energies, $`U_s`$, versus index densities, $`n_s`$ (b) saddle index densities as a function of the reduced temperature, $`T^{}`$ and (c) saddle configurational energies, $`U_s`$, as a function of the reduced temperature, $`T^{}`$. Configurational energies $`U`$ are in units of $`ϵ`$ per particle.
4. Sensitivity of the lowest magnitude, $`|\omega _0^2|`$, and third lowest magnitude, $`|\omega _2^2|`$ eigenvalues of the Hessian for a set of 20 saddle configurations sampled from the simulation at $`T^{}=1.05`$ for for $`\tau =10^{14}`$ and $`10^{15}`$. The y-axis has a logarithmic scale.
5. Scatter plots showing the correlation of the $`|\omega _0^2|/\omega _E^2`$ values with the value of the pseudo-potential, $`W(𝐫)=|U|^2`$, for saddle configurations sampled from simulations at two different temperatures where $`\omega _E`$ is the Einstein frequency.