Mean-field solution of the Hubbard model: the magnetic phase diagram

The Hubbard Hamiltonian in the simplest case of a non-degenerate band [1], i.e., with one orbital per site, can be expressed in the second-quantization formalism [4] as:

$$\begin{equation} {\hat{H}}_H=-\sum _{ij\sigma }t_{ij}{\hat{c}}^{\dagger }_{i\sigma }{\hat{c}}_{j\sigma } + U\sum _i {\hat{n}}_{i\uparrow }{\hat{n}}_{i\downarrow }\equiv {\hat{H}}_t + {\hat{H}}_U. \end{equation} \tag{ 1 }$$

Here, as usual, ${\hat{c}}^{\dagger }_{i\sigma }$ is an operator representing the creation of an electron of spin σ (=↑, ↓) at site i, ${\hat{c}}_{j\sigma }$ is the annihilation of an electron of spin σ at site j and t_ij is the amplitude of the process, the so-called hopping amplitude from site j, where the electron is destroyed, to site i, where the electron is created. We finally notice that the sum over i, j is unrestricted and t_ij = (t_ji)*, see equation (2): therefore the Hamiltonian is Hermitian, as it should be. The term t_ij is the translation in the second-quantization language of both the kinetic energy and the crystal-potential energy associated with an electron at site i:

$$\begin{equation} t_{ij} \equiv t_{ij,\sigma } = \int \mathrm{d}^3\vec{r} \ \varphi ^*_{\vec{R}_i,\sigma }(\vec{r}) \left( -\frac{\hbar ^2\nabla ^2}{2m} +V(\vec{r}) \right) \varphi _{\vec{R}_j,\sigma }(\vec{r}). \end{equation} \tag{ 2 }$$

Wannier wave functions [5] $\varphi _{\vec{R}_i,\sigma }(\vec{r})$ are centred at site $\vec{R}_i$. We suppose that the hopping amplitude does not depend on the spin variable and therefore we drop the σ label. The term $V(\vec{r})$ represents the periodic crystal-potential energy. As stated above, one approximation that is usually employed for the hopping term t_ij is to consider it as being different from zero only when i and j are nearest-neighbour sites (t_ij = t), as their overlap is usually the largest. In this case the conventional minus sign associated to the hopping term in equation (1) allows having a minimum at the Γ-point in the reciprocal space ($\vec{k}=0$) for positive t.

Operators ${\hat{n}}_{i\sigma } \equiv {\hat{c}}^{\dagger }_{i\sigma }{\hat{c}}_{i\sigma }$ count the number of particles at site i with spin σ. They are projection operators, i.e., ${\hat{n}}^2_{i\sigma } = {\hat{n}}_{i\sigma }$ (either there is one electron with spin σ at site i or there are no electrons). The expectation value $\langle {\hat{n}}_{i\sigma }\rangle \equiv \langle \Psi _0|{\hat{n}}_{i\sigma }|\Psi _0\rangle$ in the many-body ground-state |Ψ₀〉 represents the electron density n_iσ at site i with spin σ (mean occupation number). The physical origin of ${\hat{H}}_U$ is the Coulomb repulsion of the electrons: when at site i both spin-up and spin-down electrons are present, from equation (1) they contribute to the total energy with a term +U, as both ${\hat{n}}_{i\uparrow }$ and ${\hat{n}}_{i\downarrow }$ in ${\hat{H}}_U$ give one. If, on the other side, the two electrons belong to two separate atoms, they do not feel any Coulomb repulsion (this is of course a strong constraint). Formally, the on-site Coulomb repulsion can be written as:

$$\begin{eqnarray} &&U \equiv \frac{e^2}{4\pi \varepsilon _0}\int \mathrm{d}^3\vec{r} \ \mathrm{d}^3\vec{r} \ ^{\prime } \ |\varphi _{\vec{R}_i,\sigma }(\vec{r}) |^2 \frac{1}{ |\vec{r}-\vec{r}\,^{\prime } |} |\varphi _{\vec{R}_i,\bar{\sigma }}(\vec{r} \ ^{\prime }) |^2. \end{eqnarray} \tag{ 3 }$$

The Coulomb term does not depend on the site label, i, if we suppose the system homogeneous. Notice that in the extreme limit U/t = 0, we recover a purely band-like (tight-binding) picture, with just kinetic energy and crystal periodic potential, and in the opposite limit t/U = 0, we find a purely atomic picture.

The bi-dimensional square lattice is drawn in figure 1. The unit cell is spanned by the two vectors $\vec{a}_1$ and $\vec{a}_2$, of common length a. However, for our calculations, we consider the cell of double area spanned by the two vectors $\vec{b}_1$ and $\vec{b}_2$, with the idea of looking for possible antiferromagnetic ground-states. Such a cell encloses two atomic sites that can in principle be inequivalent (e.g., spin ↑ and spin ↓). In what follows, we adopt the site label $\vec{R}_i$ for the double unit cell and, within each cell, the two atoms are labelled by an extra index α = 1, 2. For example, as shown in figure 1, nearest neighbours of α = 1 sites are necessarily α = 2 and vice-versa (bipartite lattice). The general creation (annihilation) operator for an electron at site $\vec{R}_i$, position α, spin σ can be written as: ${\hat{c}}^{\dagger }_{i\alpha \sigma }$ (${\hat{c}}_{i\alpha \sigma }$). Equation (1) should be changed accordingly:

$$\begin{equation} {\hat{H}}_{H}=-\sum _{ij\alpha \alpha ^{\prime }\sigma }t_{ij}^{\alpha \alpha ^{\prime }}{\hat{c}}^{\dagger }_{i\alpha \sigma }{\hat{c}}_{j\alpha ^{\prime }\sigma } + U\sum _{i\alpha } {\hat{n}}_{i\alpha \uparrow }{\hat{n}}_{i\alpha \downarrow }\equiv {\hat{H}}_{\tilde{t}} + {\hat{H}}_{\tilde{U}}. \end{equation} \tag{ 4 }$$

Zoom In Zoom Out Reset image size

The structure constants $t_{ij}^{\alpha \alpha ^{\prime }}$ for our problem (nearest-neighbour hopping) can be written as follows:

$$\begin{eqnarray} &&t_{ij}^{\alpha \alpha ^{\prime }} = - t\delta _{\alpha ,1}\delta _{\alpha ^{\prime },2} [ \delta _{\vec{R}_i,\vec{R}_j} \!\!\!+ \delta _{\vec{R}_i+\vec{b}_1,\vec{R}_j} \!\!\!+ \delta _{\vec{R}_i+\vec{b}_2,\vec{R}_j} \!\!\!+ \delta _{\vec{R}_i+\vec{b}_1+\vec{b}_2,\vec{R}_j} ] \nonumber \\ &&-\, t\delta _{\alpha ^{\prime },1} \delta _{\alpha ,2} [ \delta _{\vec{R}_i,\vec{R}_j}\!\!\! + \delta _{\vec{R}_i-\vec{b}_1,\vec{R}_j} \!\!\!+ \delta _{\vec{R}_i-\vec{b}_2,\vec{R}_j} \!\!\!+ \delta _{\vec{R}_i-\vec{b}_1-\vec{b}_2,\vec{R}_j} ]. \end{eqnarray} \tag{ 5 }$$

In spite of the cumbersome form, the meaning of terms appearing in (5) is quite straightforward: the first line represents the four nearest-neighbour hopping energies (represented by dashed green lines in figure 1) when site $\vec{R}_i$ is of α = 1-type and site $\vec{R}_j$ is of α = 2-type, and the second line the opposite case. As we employ the reciprocal lattice of the ($\vec{b}_1$, $\vec{b}_2$) direct lattice, we must express the nearest neighbours in terms of linear combinations of $\vec{b}_1$, $\vec{b}_2$ vectors.

Consider first the hopping part of the Hamiltonian, ${\hat{H}}_{\tilde{t}}$. By inserting the $\vec{k}$-transformed operators, we get the usual tight-binding band structure [5]:

$$\begin{eqnarray} {\hat{H}}_{\tilde{t}} & = -\frac{1}{N} \sum _{ij\alpha \alpha ^{\prime }\sigma } t_{ij}^{\alpha \alpha ^{\prime }}\sum _{\vec{k}\vec{k}^{\prime }} \mathrm{e}^{-\mathrm{i}\vec{k}\cdot \vec{R}_i}{\hat{c}}^{\dagger }_{\vec{k}\alpha \sigma } \mathrm{e}^{\mathrm{i}\vec{k}^{\prime }\cdot \vec{R}_j}{\hat{c}}_{\vec{k}^{\prime }\alpha ^{\prime }\sigma } \nonumber \\ & = -\frac{1}{N}\sum _{\vec{k}\vec{k}^{\prime }\alpha \alpha ^{\prime }\sigma }{\hat{c}}^{\dagger }_{\vec{k}\alpha \sigma }{\hat{c}}_{\vec{k}^{\prime }\alpha ^{\prime }\sigma } \sum _{ij} t_{ij}^{\alpha \alpha ^{\prime }}\mathrm{e}^{-\mathrm{i}\vec{k}\cdot \vec{R}_i}\mathrm{e}^{\mathrm{i}\vec{k}^{\prime }\cdot \vec{R}_j} \nonumber \\ & = -\sum _{\vec{k}\vec{k}^{\prime }\alpha \alpha ^{\prime }\sigma }{\hat{c}}^{\dagger }_{\vec{k}\alpha \sigma }{\hat{c}}_{\vec{k}^{\prime }\alpha ^{\prime }\sigma } \sum _{j}t_{\vec{\eta }_j}^{\alpha \alpha ^{\prime }}\mathrm{e}^{\mathrm{i}\vec{k}^{\prime }\cdot \vec{\eta }_j} \frac{1}{N} \sum _i \mathrm{e}^{-\mathrm{i}(\vec{k}-\vec{k}^{\prime })\cdot \vec{R}_i} \nonumber \\ & = \sum _{\vec{k}\alpha \alpha ^{\prime }\sigma }{\hat{c}}^{\dagger }_{\vec{k}\alpha \sigma }{\hat{c}}_{\vec{k}\alpha ^{\prime }\sigma } \varepsilon _{\vec{k}}^{\alpha \alpha ^{\prime }} \end{eqnarray} \tag{ 6 }$$

where we defined the matrix Hamiltonian in the α-α' basis as: $\varepsilon _{\vec{k}}^{\alpha \alpha ^{\prime }}= -\sum _{j}t_{\vec{\eta }_j}^{\alpha \alpha ^{\prime }}\mathrm{e}^{\mathrm{i}\vec{k}\cdot \vec{\eta }_j}$. In passing from the second to the third line of (6) we used the translational invariance of the structure factor $t_{ij}^{\alpha \alpha ^{\prime }}$. This means that if we write $\vec{R}_j=\vec{R}_i+\vec{\eta }_{ij}$, the vectors $\vec{\eta }_{ij}$, and therefore the way of counting nearest neighbours, are independent of the starting point $\vec{R}_i$. So, we can write: $\vec{\eta }_{ij}\rightarrow \vec{\eta }_{j}$. Moreover we used the fact that $\frac{1}{N} \sum _i \mathrm{e}^{-\mathrm{i}(\vec{k}-\vec{k}^{\prime })\cdot \vec{R}_i}=\delta _{\vec{k}\vec{k}^{\prime }}$.

As we have two atoms per unit cell, the matrix Hamiltonian is a 2 × 2 matrix (α, α' = 1, 2) that should be diagonalized in order to have the band structure $\varepsilon ^{\pm }_{0\vec{k}}$. From the explicit expression of the structure constants (5), we get the matrix:

$$\begin{eqnarray} H_{\tilde{t}}&=\left[\begin{array}{@{}cc@{}}0 & t\gamma _{\vec{k}} \\ t\gamma ^*_{\vec{k}} & 0 \end{array} \right] \end{eqnarray} \tag{ 7 }$$

where $\gamma _{\vec{k}}=- \lbrace 1+{\rm e}^{-{\rm i}\vec{k}\cdot (\vec{b}_1+\vec{b}_2)} +{\rm e}^{-{\rm i}\vec{k}\cdot \vec{b}_1}+{\rm e}^{-{\rm i}\vec{k}\cdot \vec{b}_2}\rbrace$. Its diagonalization leads to two bands given by:

$$\begin{eqnarray} &&\varepsilon ^{\pm }_{0\vec{k}} = \pm t |\gamma _{\vec{k}}| = \pm t \left|2\cos \left(\vec{k}\cdot \frac{\vec{b}_1+\vec{b}_2}{2}\right)+2\cos \left(\vec{k}\cdot \frac{\vec{b}_2-\vec{b}_1}{2}\right)\right|. \end{eqnarray} \tag{ 8 }$$

The band structure and density of states (DOS) per unit cell are depicted in figure 2, together with those for the single cell ($\vec{a}_1$, $\vec{a}_2$), for comparison. Bandwidth is W = 8t. At half filling, we have a nested Fermi surface, leading to a van Hove logarithmic singularity in the DOS [5]. We remind readers that nesting is the property by which any point of the Fermi surface is related to another point of the Fermi surface by a fixed vector, in this case the vector (π/a, π/a), because the Fermi surface is a square at half filling, as shown in figure 2(b). Clearly, at this level, changing the choice of the unit cell (a pure convention) does not change our results. In fact, the second band is just the folding of the first band of the single cell, as can be also seen by considering that $\frac{\vec{b}_1+\vec{b}_2}{2}=\vec{a}_1$ and $\frac{\vec{b}_2-\vec{b}_1}{2}=\vec{a}_2$. This comes from the fact that the point M_a ≡ (π/a, π/a) in the reciprocal cell of the direct cell ($\vec{a}_1,\vec{a}_2$) becomes equivalent to Γ_b ≡ (0, 0) in the reciprocal cell of the direct cell ($\vec{b}_1,\vec{b}_2$). The two DOS of figure 2 are the double of each other just because they are normalized per unit cell and there are two atoms per unit cell in the case of the double cell.

Zoom In Zoom Out Reset image size

The presence of ${\hat{H}}_{\tilde{U}}$, equation (4), changes these simple results even at a mean-field level. We are reminded that the mean-field approximation corresponds to neglecting the fluctuations around the mean density. Such fluctuations are defined as $\Delta {\hat{n}}_{i\alpha \sigma } \equiv {\hat{n}}_{i\alpha \sigma } - \langle {\hat{n}}_{i\alpha \sigma }\rangle$, i.e., the difference between the exact number operator ${\hat{n}}_{i\alpha \sigma }$ and the mean occupation number $\langle {\hat{n}}_{i\alpha \sigma }\rangle$. From this relation, we get: ${\hat{n}}_{i\alpha \sigma } \equiv \langle {\hat{n}}_{i\alpha \sigma }\rangle + \Delta {\hat{n}}_{i\alpha \sigma }$. If we suppose homogeneity of the system, then $\langle {\hat{n}}_{i\alpha \sigma }\rangle$ is independent of the cell position $\vec{R}_i$, and we can drop the label i and write:

$$\begin{eqnarray} {\hat{n}}_{i\alpha \uparrow } {\hat{n}}_{i\alpha \downarrow } & = [ \Delta {\hat{n}}_{i\alpha \uparrow } + \langle {\hat{n}}_{\alpha \uparrow }\rangle ] \cdot [ \Delta {\hat{n}}_{i\alpha \downarrow } + \langle {\hat{n}}_{\alpha \downarrow }\rangle ] \nonumber \\ & = \Delta {\hat{n}}_{i\alpha \uparrow } \Delta {\hat{n}}_{i\alpha \downarrow } + \sum _{\sigma } {\hat{n}}_{i\alpha \sigma } \langle {\hat{n}}_{\alpha \bar{\sigma }} \rangle - \langle {\hat{n}}_{\alpha \uparrow }\rangle \langle {\hat{n}}_{\alpha \downarrow }\rangle . \end{eqnarray} \tag{ 9 }$$

In a mean-field approximation we can neglect the first term of the previous equation, quadratic in the fluctuations. This implies that ${\hat{H}}_{\tilde{U}}$ becomes:

$$\begin{eqnarray} &&{\hat{H}}^{\mathrm{MF}}_{\tilde{U}} = \underbrace{U \sum _{i\alpha \sigma } {\hat{n}}_{i\alpha \sigma } \langle n_{\alpha \bar{\sigma }} \rangle }_{\tilde{H}_U} - \underbrace{U N \sum _{\alpha } \langle n_{\alpha \uparrow }\rangle \langle n_{\alpha \downarrow }\rangle }_{E_U}. \end{eqnarray} \tag{ 10 }$$

The second term E_U is a constant for a given magnetic configuration and number of particles, as it does not depend on creation or annihilation operators but only on their average values. However, it must be integrated in the calculation of the magnetic phase diagram, as such a term is advantageous for paramagnetic (PM) configurations with respect to ferromagnetism and antiferromagnetism (E_U is negative and the product 〈n_α↑〉〈n_α↓〉, for fixed number of particles per site, is maximum when 〈n_α↑〉 = 〈n_α↓〉).

By using the Fourier transform as for equation (6), the first term of (10) becomes: $\tilde{H}_U = U \sum _{i\alpha \sigma } \langle n_{\alpha \bar{\sigma }}\rangle {\hat{c}}^{\dagger }_{i\alpha \sigma } {\hat{c}}_{i\alpha \sigma }=U \sum _{\vec{k}\alpha \sigma } \langle n_{\alpha \bar{\sigma }}\rangle {\hat{c}}^{\dagger }_{\vec{k}\alpha \sigma } {\hat{c}}_{\vec{k}\alpha \sigma }$, diagonal in α. Therefore the energy per $\vec{k}$-point for a spin-σ electron is obtained by diagonalizing the 2 × 2 matrix:

$$\begin{eqnarray} &\left[\begin{array}{@{}cc@{}}U\langle n_{1\bar{\sigma }}\rangle & t\gamma _{\vec{k}} \\ t\gamma ^*_{\vec{k}} & U\langle n_{2\bar{\sigma }} \rangle \end{array} \right], \end{eqnarray} \tag{ 11 }$$

where $\gamma _{\vec{k}}$ has been defined after (7). This Hamiltonian can be easily diagonalized analytically, as the associated eigenvalue problem leads to the second-order algebraic equation for the eigenenergies $\varepsilon ^{\pm }_{\vec{k}}$:

$$\begin{eqnarray} &&\varepsilon ^{2}_{\vec{k}} - \varepsilon _{\vec{k}} \left(U\langle n_{1\bar{\sigma }}\rangle +U\langle n_{2\bar{\sigma }}\rangle \right) + U^2 \langle n_{1\bar{\sigma }}\rangle \langle n_{2\bar{\sigma }}\rangle -t^2|\gamma _{\vec{k}}|^2 =0. \end{eqnarray} \tag{ 12 }$$

The solutions of this second-order equation provide the mean-field band energies per $\vec{k}$-point for a spin-σ electron, once we add again the term E_U (divided by N, because $\varepsilon _{\vec{k}}$ is the energy per spin-σ electron):

$$\begin{eqnarray} &&\varepsilon ^{\pm }_{\vec{k}\sigma } = U\left(\frac{\langle n_{1\bar{\sigma }}\rangle +\langle n_{2\bar{\sigma }}\rangle }{2}\right) -U (\langle n_{1\uparrow }\rangle \langle n_{1\downarrow }\rangle + \langle n_{2\uparrow }\rangle \langle n_{2\downarrow }\rangle ) \nonumber \\ &&\pm \frac{1}{2}\textstyle {\sqrt{U^2(\langle n_{1\bar{\sigma }}\rangle -\langle n_{2\bar{\sigma }}\rangle )^2+16t^2\left(\cos \left[\vec{k}\cdot \frac{\vec{b}_1+\vec{b}_2}{2}\right]+\cos \left[\vec{k}\!\cdot \frac{\vec{b}_1-\vec{b}_2}{2}\right]\right)^2}} \nonumber \\ &&= U\left(\frac{\langle n_{1\bar{\sigma }}\rangle +\langle n_{2\bar{\sigma }}\rangle }{2}\right) -U (\langle n_{1\uparrow }\rangle \langle n_{1\downarrow }\rangle + \langle n_{2\uparrow }\rangle \langle n_{2\downarrow }\rangle ) \nonumber \\ &&\pm \frac{1}{2}\textstyle {\sqrt{U^2 (\langle n_{1\bar{\sigma }}\rangle -\langle n_{2\bar{\sigma }}\rangle)^2+4 \big(\varepsilon ^{\pm }_{0\vec{k}}\big)^2}} . \end{eqnarray} \tag{ 13 }$$

From the last line we see that ${\hat{H}}^{{\rm MF}}_{\tilde{U}}$ contributes to the total energy as a density-dependent and spin-dependent (because of $\bar{\sigma }$) renormalization of the tight-binding energy $\varepsilon ^{\pm }_{0\vec{k}}$. This double dependence on the density and on the spin is at the basis of the richness in the phase diagram shown in figure 6, allowing to get PM, ferromagnetic (FM) and antiferromagnetic (AFM) phases.

In order to find the ground-state energy for a given magnetic configuration, we need to employ a self-consistent scheme as, for a given spin σ, the energy depends on the mean occupation number of opposite spin, $\langle {\hat{n}}_{\alpha \bar{\sigma }}\rangle$, which in turns depends on the eigenvectors of the energy matrix itself. As shown in figure 3, we therefore start from a given configuration of input parameters (frames 1 and 2), we follow steps 3–6, up to the self-consistency condition expressed in frames 7 and 8, and finally move back to frame 2 or end to frame 9, depending on whether the condition is not satisfied or satisfied, respectively. The condition is satisfied when the output occupation number is the same as the input occupation number within a threshold that we fixed to 10⁻⁵.

Zoom In Zoom Out Reset image size

In more detail, we proceeded as follows: for practical reasons we found it simpler to numerically diagonalize the Hamiltonian of (11) through the Lapack libraries [6] (frame 3) in order to determine eigenvalues $\varepsilon _{j\sigma }(\vec{k})$ and eigenvectors $|\Psi ^{\,j}_{\vec{k}\sigma }\rangle =\sum _{\alpha }A^{\,j}_{\alpha ,\sigma }(\vec{k}) |\alpha \vec{k}\rangle$ (with j = 1, 2). We actually work with the chemical potential in order to fix the particle density a posteriori: ${\hat{H}}_H^{\mu }={\hat{H}}_H-\mu \sum _{i\alpha \sigma }{\hat{n}}_{i\alpha \sigma }$. The calculations performed in frames 4–6 are based on equations (14)–(16) given below. The first of these equations expresses the DOS per unit cell, ρ(ε):

$$\begin{equation} \rho (\varepsilon )=\frac{1}{N}\sum _{{\vec{k}}, j, \sigma } \delta [\varepsilon -\varepsilon _{j \sigma }({\vec{k}})]. \end{equation} \tag{ 14 }$$

The Dirac delta function in (14) is numerically calculated through a Gaussian function whose broadening is optimized after a convergence study. If the width of the Gaussian is too wide compared to the step between two energy points in ε, the DOS is too smooth, whereas if the width is too narrow, artificial oscillations appear. A more elegant approach might be to use the Methfessel–Paxton method [7], usually implemented in more advanced packages for electronic structure calculations.

As a second step we determine the chemical potential μ by the implicit equation:

$$\begin{equation} \int _{-\infty }^{+\infty } \mathrm{d}\varepsilon \,\rho (\varepsilon ) \frac{1}{\exp [\beta (\varepsilon -\mu )]+1}=n_e. \end{equation} \tag{ 15 }$$

Here β = 1/k_BT is the Boltzmann factor and we have defined the total number of electrons per unit cell (as used in figures 3, 5, 6): n_e ≡ ∑_{α, σ}〈n_ασ〉. Finally, the average number of particles per site α and per spin σ at a given temperature T is given by:

$$\begin{equation} \langle n_{\alpha \sigma } \rangle =\frac{1}{N} \sum _{{\vec{k}}, j=1,2} |A_{\alpha , \sigma }^{\,j}({\vec{k}})|^2 \frac{1}{\exp [\beta (\varepsilon _{j\sigma }({\vec{k}})-\mu )]+1}. \end{equation} \tag{ 16 }$$

At finite temperature, the thermodynamic variable to minimize is of course not the total energy, E, but the free energy, F = E − TS, what implies a calculation of the entropy of the system, through the equation:

$$\begin{equation} S(T) = - k_B \int _{-\infty }^{\infty } \mathrm{d}\varepsilon \rho (\varepsilon ) \lbrace f( \varepsilon ) \ln [f(\varepsilon )] \, + \, [1-f(\varepsilon )] \ln [1-f(\varepsilon )] \rbrace . \end{equation} \tag{ 17 }$$

Here f(ε) is the Fermi–Dirac distribution. However, in what follows, we are interested in the ground-state phase diagram of the model, i.e., at T = 0. Therefore, we minimized the total energy and used the parameter T of equations (15) and (16) for convergence purposes only. It is in fact common practice to use a Fermi–Dirac distribution, characterized by a parameter T different from zero, even in the T → 0 limit, in order to smooth the step function that appears in equations (15) and (16) for T = 0.

In general, we computed the total energy for each magnetic phase as a function of t/U. A phase transition is then characterized by the crossing of two (or more) free-energy curves: for example, we have represented the magnetic free energy as a function of t/U in figure 4 for a specific occupation number (n = 0.8), in order to highlight a magnetic (AFM/FM) transition. In this case the phase transition from AFM to FM phase is obtained at t/U = 0.13, the value at which the total FM energy becomes smaller than the AFM energy.

Zoom In Zoom Out Reset image size

In figure 5, we have drawn the band structures for both FM and AFM phases to highlight the differences with the PM case. FM bands are shifted rigidly with respect to PM bands by about ±Uδn, where δn is half the spin unbalance (see below). AFM bands are instead characterized by the opening of a gap that can be related to the difference $U(| \langle n_{1\bar{\sigma }}\rangle \!-\!\langle n_{2\bar{\sigma }}\rangle |)$ in (13). This gap leads to an insulating phase for n_e = 2.0. It is important to notice that such a gap is not related to a metal–insulator transition (MIT) of the Mott–Hubbard kind, as it is not related to the electronic correlations (that are absent by definition in a mean-field calculation), but to magnetism. This MIT is called a Slater MIT, in honour of Slater, who foresaw it in 1951 [8]. In fact, unlike Mott, who did not originally ascribe his MIT to magnetic interactions, Slater thought that the origin of the metal-to-insulator transition was determined by the onset of AFM long-range order, exactly as in the scheme described in the present paper. Therefore a Slater insulator is characterized by a band gap determined by a superlattice modulation of the magnetic periodicity. This is not the case of a Mott insulator.

Zoom In Zoom Out Reset image size

Our main result is the ground-state phase diagram, drawn in figure 6 as a function of the number of electrons per site and of t/U. Such a phase diagram has been already obtained in the literature in 1985 by Hirsch [9], though in a different context and without providing all the details of the derivation that can be found here. As a general feature, the phase diagram shows a clear symmetry around half filling, i.e., one electron per site, where antiferromagnetism is the lowest-energy configuration. This symmetry had to be expected, since it is a symmetry of the Hubbard Hamiltonian for nearest-neighbour hopping. Far from half filling, paramagnetism has the advantage of a high value of t/U, whereas a low value of t/U leads instead to ferromagnetism. This tendency for the PM/FM phases can be easily understood: the ground-state of n non-interacting electrons (U = 0) is PM because the minimum-energy constraint in combination with the Pauli principle forces to fill all the energy levels from the lowest (ε_min) to the highest (μ) with an equal number of n/2 up and n/2 down electrons. In the opposite extreme case, for U/t → ∞, the system can gain energy by a total magnetization (say, all electrons with spin up), in order to minimize the energy term U〈n_↓〉. In the intermediate t/U cases, only the numerical study of equations (11), (14)–(16) can provide us with the magnetic phase of the system.

Zoom In Zoom Out Reset image size

There exists however a criterion that allows us to foresee the stability of the PM phase versus the FM phase just on the base of the two parameters U and ρ(μ), the DOS at the Fermi energy (ε_F = μ at T = 0). This is the Stoner criterion [10]. Start from the PM phase, represented in grey in figure 7(a) (n_↑ = n_↓ = n/2) and move δn electrons from spin-down to spin-up states, so that n_↑ = n/2 + δn and n_↓ = n/2 − δn. This leads to a shift of the chemical potential by ±δε ≃ ±δn/ρ(μ) (if δn ≪ n). The system is not at equilibrium, as μ_↑ = μ + δn/ρ(μ) ≠ μ_↓ = μ − δn/ρ(μ). Yet, because of the change in δn and of the mean-field Coulomb energy (10), the minority-spin band shifts as a whole upwards by Uδn + U(δn)²/n and majority-spin band shifts as a whole downwards by Uδn − U(δn)²/n as shown in figure 7(b).

Zoom In Zoom Out Reset image size

In this situation, if ${\tilde{\mu }_{\uparrow }} = \mu + \delta n/\rho (\mu )-U\delta n + U (\delta n)^2/n$ is less than ${\tilde{\mu }_{\downarrow }}=\mu - \delta n/\rho (\mu )+U\delta n + U (\delta n)^2/n$, then it is favourable to still increase the number of spin-up electrons until ${\tilde{\mu }_{\uparrow }} = {\tilde{\mu }_{\downarrow }}$. Therefore, ferromagnetism appears when ${\tilde{\mu }_{\uparrow }} \le {\tilde{\mu }_{\downarrow }}$. This leads to the Stoner criterion for ferromagnetic stability:

$$\begin{equation} U\rho (\mu ) \ge 1. \end{equation} \tag{ 18 }$$

Equation (18) can also be written as $\frac{t}{U} \le t \rho (\mu )$. The corresponding equality, that marks the phase transition, has been reproduced in figure 6. All our calculated points for the PM/FM transition lie on the theoretical curve $\frac{t}{U} = t \rho (\mu )$ represented by a dotted line in figure 6. We infer from the Stoner criterion that an FM instability is expected in materials showing a high DOS at Fermi level. This is indeed the case for Fe, Co and Ni [11]. It is also possible to find a similar criterion for the AFM/PM and AFM/FM transitions, but its derivation is technically more involved because it is based on a Bogoliubov transformation. This leads to a gap equation formally equivalent to the one of the BCS theory of superconductivity. A rather detailed description of this generalization can be found in [12].