Magnetic phase diagram of the Hubbard model with next-nearest-neighbour hopping

We calculate the magnetic phase diagram of the Hubbard model for a Bethe lattice with nearest neighbour (NN) hopping $t_1$ and next nearest neighbour (NNN) hopping $t_2$ in the limit of infinite coordination. We use the amplitude $t_2/t_1$ of the NNN hopping to tune the density of states (DOS) of the non-interacting system from a situation with particle-hole symmetry to an asymmetric one with van-Hove singularities at the lower ($t_2/t_1>0$) respectively upper ($t_2/t_1<0$) band edge for large enough $|t_2/t_1|$. For this strongly asymmetric situation we find rather extended parameter regions with ferromagnetic states and regions with antiferromagnetic states.


Introduction
The phenomenon of spontaneous magnetism is one of the oldest topics in physics. That lodestone can attract iron is known for over 2500 years. In contrast, a rigorous understanding of the microscopic processes which lead to magnetism still is a matter of present day research [1]. In order to microscopically describe the phenomenon "magnetism" quantum mechanics, in particular the spin of the electrons, and the inclusion of interactions respectively many-body correlations are mandatory.
A further typical property of materials which show magnetic behaviour is that they possess partially filled d-or f-shells. In this case, orbital degrees of freedom usually quite dramatically influence the existence and nature of magnetically ordered states. A rather notorious example are the manganites, which show a rather complex phase diagram due to an interplay of orbital and spin degrees of freedom [2,3].
A much simpler situation occurs when, for example, in a crystal with low symmetry due to lattice distortions, only one of the d-or f-states effectively plays a role at the Fermi energy. In this case one can think of an effective one-band model as appropriate description. A well-known example for such a situation are the cuprate superconductors [4]. Here, too, magnetic order can occur. However, while for materials with orbital degrees of freedom the existence of both antiferromagnetism and ferromagnetism can easily be accounted for [4], the one-band situation prefers the formation of antiferromagnetic order [5]. While ferromagnetic states are known to exist under certain extreme conditions [6], their possible occurrence and stability regions in physically relevant model parameter regimes is still an intensively investigated research topic.
In this paper we therefore want to focus on the one-orbital situation. A suitable model for describing strong correlation physics in such a single band is the Hubbard model [7,8,9] The operator c † iσ creates an electron with spin σ at site i, t ij describes the "hopping" amplitude from site i to j and µ is the chemical potential, which can be used for tuning the occupation of the system. The two particle interaction is purely local and only entering via a product of two density operators n i↑ = c † i↑ c i↑ with amplitude U. In recent years progress in understanding the physics of this model in dimensions larger than one was mostly gained from calculations using the dynamical mean field theory (DMFT) [10,11] or cluster variants of it [12]. The DMFT relates the lattice model to an impurity model in an effective medium representing the lattice, which must be solved self-consistently. It can be shown that this mapping is exact in the limit of infinite spatial dimensions or infinite coordination of the lattice [10,13]. Note that the remaining (effective) impurity problem represents a quantum-impurity, which by itself is complicated to solve. From the methods available we here use the numerical renormalisation group (NRG) [14,15], because it is by far the most efficient and accurate technique for single-band problems. For the calculation of spectral functions we employ the complete Fock space variant [16,17] of the NRG.
For real three dimensional materials the DMFT is, of course, only an approximation. Nevertheless, the Hubbard model within DMFT describes a lot of strong correlation physics, which can be seen in real materials, at least qualitatively correct. In this sense it is therefore justified to study for example magnetic properties of the Hubbard model within this approximation. As the DMFT can be seen as a thermodynamically consistent mean-field theory [10,18], one can expect that the phase diagram obtained at least gives an account for potential phases, albeit not necessarily the correct phase boundaries.
The aim of the present paper is to give an account of the possible antiferromagnetic and ferromagnetic phases of the doped single-band Hubbard model. For a particlehole symmetric density of states (DOS) the model has an antiferromagnetically ordered ground state at half filling for every finite value of U, which phase separates upon doping [19,20,22]. Ferromagnetism can also be found in the single band Hubbard model, but only for very high interaction parameter and close to half filling [6,21,22,23], or for a pronounced asymmetric DOS also for moderate values of U [24,25].
Deviations from particle-hole symmetry in the single-band model leading to such asymmetries in the DOS are achieved by inclusion of longer-range single-particle hopping processes. It is important to stress that in DMFT the actual lattice structure only enters via the DOS. As we are interested in a qualitative investigation of the possible magnetic phases, it is permissible to work with a computationally convenient DOS, which is the one obtained from an infinitely-coordinated Bethe lattice [10] with nearest neighbour (NN) and next-nearest neighbour (NNN) hopping t 1 respectively t 2 . For t 2 = 0 one obtains the well-known semicircular DOS [10], which for values t 2 > 0 becomes asymmetric and can even develop a singularity at one of the band edges [26,27]. From this point of view, the Bethe lattice in the limit of infinite coordination has all typical features of the DOS of a real lattice -compact support, van-Hove singularities -and one can hope that results obtained with it give a reasonable qualitative account of true three-dimensional systems.
The paper is organised as follows. In the next section we introduce the DOS of the t 1 -t 2 Bethe lattice with infinite coordination, which will be used throughout the paper. Section three focuses on the antiferromagnetic phase, which is realised near half filling. In section four we present the results for the ferromagnetic calculations. Quite surprisingly, for strong enough t 2 we observe regions, where both antiferromagnetic and ferromagnetic states are stable. A summary and discussion will conclude the paper.

Density of States
Early studies of the Bethe lattice with longer-ranged hopping usually focused on the simplified variant proposed by Georges et al. [10,35,36]. While in this approximation one introduces frustration to magnetic correlations, the resulting DOS retains particlehole symmetry, which of course is somewhat artificial. The proper form of the DOS was deduced by Kollar et al. [26,27]. Figure 1 shows the result for different ratios of t 2 /t 1 . The non-interacting Green's function G t 1 ,t 2 (ζ) and by this the DOS ρ(ω) = −1/πℑG t 1 ,t 2 (ω + iη) for the Bethe lattice with nearest neighbour hopping t 1 and next nearest neighbour hopping t 2 in the limit of infinite coordination is given by the formula Analysing this formula shows that there appears a singularity in the DOS for t 2 > 1 4 t 1 . The singularity is due to the factor 1/b and thus is a square root singularity. For t 2 < 1 4 t 1 the band edges lie at ω 1,2 = 3t 2 ± 2t 1 . For t 2 > 1 4 t 1 the lower band edge is upper band edge ω 2 = 3t 2 + 2t 1 . Thus the bandwidth is It should be emphasised that by tuning the NNN hopping t 2 , the DOS change from a particle-hole symmetric semi-ellipse to a strongly asymmetric DOS with singularity for t 2 /t 1 > 1 4 . This is a rather important feature expected to occur also in real materials. On the other hand, previous investigations of frustration effects within DMFT used the so-called two sub-lattice fully frustrated model [10,33,34,35,36], which misses this particular asymmetry and the van-Hove singularity.

Magnetic phases close to half filling
Before discussing the magnetic phases within DMFT of the system with finite t 2 , let us briefly review the results for the case t 2 = 0. Figure 2 shows the Néel-and the paramagnetic state around half filling. The Néel state does only exist exactly at half filling. For interaction strengths below the critical value of the paramagnetic metal insulator transition U M IT (black line in the left panel) we find phase separation between the Néel state and the paramagnetic state, which can be seen in the right panel of figure 2. There tuning the chemical potential leads to a jump in magnetisation and occupation. For larger values of the interaction U > U M IT there is a parameter region, where our calculations do not converge (c.f. also [22]). If one looks at the occupation and magnetisation as function of DMFT-iteration, they show an oscillatory behaviour with periods longer than two. Motivated by similar previous observations [37] we interpret such a behaviour as indication that an incommensurate spin spiral is the proper magnetic state. Note that within a simple AB lattice calculation such a spin-spiral cannot be stabilised, and consequently calculations do not converge in this parameter region. As we cannot determine the nature of the magnetic order, we left this region blank in figure 2. Apparently, where for the paramagnet at half filling the metal insulator transition occurs, the magnetic state of the doped system also changes from phase separated to an incommensurate structure.
A ferromagnetic state, on the other hand, cannot be stabilised for the Bethe lattice at t 2 = 0. Note that this is strikingly different from the hypercubic lattice, where for large U and small doping a Nagaoka ferromagnet occurs [22]. The explanation is that Nagaoka's state needs closed loops on the lattice, which are available for the hypercube (leading to the exponential tails), but are absent for the Bethe lattice. Thus, although in DMFT only the DOS enters the calculations, subtle differences in the structure and support may matter quite strongly for certain aspects of the physics.
As the DOS is particle-hole symmetric for t 2 = 0, the phase diagram is completely symmetric with respect to half filling.
3.2. 0 < t 2 ≤ 1/4t 1 As t 2 becomes finite the DOS becomes asymmetric and consequently the magnetic phase diagram becomes asymmetric with respect to half filling, too. However, for sufficiently small values of t 2 it will still look very similar to figure 2, with two notable exceptions: For the hole doped side of the phase diagram, the incommensurate magnetic phase sets in at smaller values of the interaction, while on the electron doped side it starts for larger values of the interaction. Thus, for electron doping, phase separation between the antiferromagnetic state at half filling and the paramagnetic state at n > 1 prevails for stronger interaction strengths. Already for t 2 /t 1 = 0.2 we found no incommensurate phase on the electron doped side for U/W < 3. As already stated previously [33,34,35,36,37], in order to stabilise the antiferromagnetic phase for a finite next-nearest neighbour hopping one needs a finite interaction strength U c > 0.
For 1/4t 1 < t 2 ≤ t 1 one obtains according to Figure 1 a strongly asymmetric DOS showing a square-root singularity at the lower band edge. Here we expect, and observe, a radically different phase diagram. As can be seen for t 2 /t 1 = 0.8 in figure 3 the Néel state can now be hole doped and does extend to large values of the doping, i.e. strong frustration seems to stabilise the Néel state. The incommensurate phase, on the other hand, completely vanished from the phase diagram. If one inspects figure 3 more closely, one sees that the antiferromagnetic state actually sets in away from half filling for increasing interaction strength. At half filling we find for this values of interaction a paramagnetic metal. On the electron doped side, we only find a paramagnetic state, which is still phase separated to the antiferromagnetic state at half filling.
As discussed in our previous work for half filling [37], for very large t 2 /t 1 > 0.96 there appears a new phase which, motivated by a 120 • order expected for a classical spin system at this level of frustration, we interpreted as such a 120 • order.  shows the phase diagram for t 2 = t 1 , i.e. a with respect to antiferromagnetic order fully frustrated spin system. The parameter region for large interaction left blank denotes precisely this 120 • state, which also can be hole doped. What is most remarkable and rather mysterious, even for the fully frustrated system we found a stable Néel state for fillings between 0.55 < n < 0.8. To ensure that this result is not a numerical artifact, we performed several calculations at different temperatures and with different NRG parameters like discretization or states kept. However, for low enough temperatures we always found this antiferromagnetic island. We will come back to this point in the last section.

Ferromagnetism
As already mentioned in the introduction, while antiferromagnetism is the "natural" order occurring in single-band systems as studied here, ferromagnetic order is usually only obtained under more restrictive conditions. In this section we therefore want to focus on possible ferromagnetic solutions in our system.
One of the first heuristic treatments of metallic ferromagnetism was by E. Stoner [38]. He gave the criterion UD F > 1 for stabilising ferromagnetism, where U is the value of the on site Coulomb interaction and D F is the value of the density of states at the Fermi level. Already in this criterion one sees that ferromagnetism is created by the interplay of the kinetic energy, characterised by D F , and the Coulomb interaction, characterised by U. A rigorous result was obtained by Nagaoka [6], who proved the existence of ferromagnetism at U = ∞ and "one hole" for certain lattices.
In the beginning of the 1990's, Mielke and Tasaki proved the existence of ferromagnetism under certain conditions on the dispersion, known as "flat band ferromagnetism" [39,40]. Here the ferromagnetic groundstate appears due to a dispersionless (flat) lowest lying band. This flat band introduces a huge degeneracy of the groundstate at U = 0, which is lifted by the Coulomb interaction. A nice overview about this topic and other rigorous results for ferromagnetism can be found in [41].
Remembering the singularity in the DOS for t 2 /t 1 > 0.25 (see figure 1), the situation present in our system is very similar to the "flat band" scenario. Former studies for an asymmetric DOS [24,25,42,43] already showed the existence of ferromagnetism in such a situation. Consequently, we have to expect ferromagnetism in our system, too. Indeed, Figure 5 shows the ferromagnetic polarisation n ↑ −n ↓ n ↑ +n ↓ colour encoded over the occupation n ↑ + n ↓ and the interaction strength at low temperature (T /W = 2 · 10 −4 ). The NNN hopping for this system is t 2 /t 1 = 0.6. One sees that the singularity in the DOS alone cannot create ferromagnetism. Here one again needs a finite interaction strength of approximately U/W ≈ 0.3, which however is a realistic number for transition metal compounds of both the 3d and 4d series. In the right panel of figure 5 we depict the lower and upper critical occupation between which the ferromagnetic state is stable as function of the interaction strength. Below the lower critical occupation, our DMFT simulations do not converge independent of the interaction strength. We believe that this is a numerical problem due to the singularity in the DOS: If the Fermi level lies very close to the singularity, the slope of the DOS at the Fermi level is very large. Small differences in the position of structures in the interacting Green's function will consequently have a great influence. We however cannot rule out the possibility of the existence of another phase in this regime. The occupation number jumps in this region between almost zero and a larger value, and cannot be stabilised. The behaviour can only be seen at low temperatures and for t 2 /t 1 > 0.25, where the singularity in the DOS is sufficiently strong and not smeared by temperature broadening.
At the upper critical occupation and low interaction strengths the system jumps from a fully polarised ferromagnet to a paramagnetic phase. For strong interaction the upper occupation is large enough such that the system directly changes from a ferromagnetic state into the incommensurate phase or the Néel phase.
As we already noted, the "flat band" scenario indicates that the ferromagnetic state is intimately connected to the appearance of the van-Hove singularity at the lower band edge. Let us therefore look somewhat closer on the relation of the formation of a ferromagnetic state and the appearance of the singularity in the DOS. Figure 6 shows the polarisation versus the occupation for different NNN hopping t 2 /t 1 and interaction strengths. The upper panels represent a situation where there is no singularity present in the DOS. The interaction needed to stabilise a ferromagnetic state in these systems without singularity is strongly increased. For the case of t 2 /t 1 = 0.2 we found no ferromagnetic phase for interactions as strong as U/W ≈ 10. As soon as t 2 /t 1 > 0.25, the critical interaction strength lies below U/W = 1. Increasing NNN hopping t 2 as well as increasing the interaction strength favours the ferromagnetic state as the region in occupation gets more and more extended. In the DMFT/QMC study of Wahle et al. [24] a peak at the lower band edge was enough to stabilise a ferromagnetic phase at moderate interaction strengths. In our calculations the tendency towards ferromagnetism dramatically decreases for a DOS without singularity.

Competition between ferromagnetism and antiferromagnetism
A careful look at the phase diagrams reveals that there are parameter regions where one seemingly can obtain both an antiferromagnetic as well as a ferromagnetic solution to the DMFT equations. This is rather unusual because conventionally DMFT will show oscillating behaviour if one performs a ferromagnetic calculation in a regime with antiferromagnetic ground state and vice versa.
To decide which of the two solutions is the thermodynamically stable one, one has to compare their respective free energies. As the calculations were done practically at T = 0, we calculate the energy of the system, given by where H T is the kinetic energy and N the number of sites. The interaction term is purely local and thus can be taken from the converged impurity calculation. The kinetic energy on the other hand can be calculated from the expression where Σ(z, θ) is the lattice self-energy, θ a suitable variable to label the single-particle energies on the lattice under consideration and µ the chemical potential. Within DMFT, the lattice self-energy is approximated by a local self-energy, i.e. we may set Σ(z, θ) = Σ(z). Furthermore, for the Bethe lattice with infinite coordination ǫ(θ) = t 1 θ + t 2 (θ 2 − 1) and ρ(θ) = 1 2π √ 4 − θ 2 holds. Substituting ǫ(θ) by ǫ in the integral, the resulting DOS takes on the form given in section 2.
Since the Néel state is defined on an AB lattice, one has to distinguish between the inter-and intra-sublattice hopping terms, and the formula for the kinetic energy takes on the form Note that with the definition of the matrix Green function this formula can be put into the compact matrix form The energies of the converged solutions for t 2 = t 1 and U/W = 2.5 can be seen in figure 7. The antiferromagnetic solution could be stabilised in this parameter region for Energies for the converged paramagnetic, ferromagnetic and antiferromagnetic solution for t 1 = t 2 and U/W = 2.5. The lines are meant as guide to the eye. occupations 0.55 < n < 0.8. From figure 7 it becomes clear now that the ferromagnetic state has the lowest energy for n < 0.6. For 0.6 < n < 0.75 the antiferromagnetic state takes over as the groundstate, but is nearly degenerate with the paramagnetic state. For fillings larger than 0.8 no staggered magnetisation can be stabilised any more.
Thus, the energy calculations reveal two things. Firstly, an antiferromagnetic Néel state indeed seems to form away from half filling in the fully frustrated system. Secondly, the energy differences are extremely small, in particular the antiferromagnet and paramagnet are de facto degenerate over the full parameter regime where the former exists.
To understand this at first rather irritating observation let us recall the well-known fact that in strongly frustrated systems it is a common feature to have a large number of degenerate groundstate configurations, which also can include magnetically ordered ones [41]. Thus, the degeneracy of the antiferromagnet and the paramagnet hints towards the possibility that there may exist a larger number of other magnetically ordered states in this parameter region. Unfortunately we are not able to search for and in particular stabilise those magnetic phases with the technique at hand. Further investigations using different methods to solve the DMFT equations are definitely necessary.

Conclusions
In conclusion, we have calculated the magnetic phase diagram for the Bethe lattice with NN-and NNN-hopping in the limit of infinite coordination. For this purpose we have used the proper expression for the DOS of this lattice as deduced by Kollar et al. . By varying the NNN hopping one can tune the DOS from a symmetric semi-ellipse to a very asymmetric shape with a square-root van-Hove singularity at the lower band edge. While the electron doped side of the phase diagram tends to phase separate between the Néel state and a paramagnetic metal just like at the particle-hole symmetric point, the hole doped side reveals a surprisingly rich phase diagram.
We first note that the regimes with phase separation respectively incommensurate spin-spiral states are replaced by a doped Néel state. As expected, we need a finite interaction U c to allow the existence of the Néel state, which for larger t 2 has its minimum at finite doping, i.e. the Néel state is first formed away from half filling.
In addition, with increasing NNN hopping t 2 a ferromagnetic phase at low fillings can be found. For large t 2 and strong interaction U this ferromagnetic phase can extend to occupations n > 0.7. The dependence of the appearance of this phase on the parameter t 2 shows that it is related to Mielke's and Tasaki's notion of "flat-band" ferromagnetism rather than Nagaoka's ferromagnetism found at low doping and U → ∞ in the hypercubic lattice.
Quite amazingly, we found that for t 2 ≈ t 1 and large enough interaction U a doped antiferromagnet can also be stabilised in the same filling region. Calculating the groundstate energies of both magnetic states and the paramagnetic solution, we find that the ferromagnet is the ground state below some critical filling n c . For n > n c , the Néel state and the paramagnet are degenerate within numerical accuracy and lower than the ferromagnet.
Finding both magnetic states stable within DMFT and the near degeneracy of them could be an effect of the strong frustration, where a large number of degenerate or nearly degenerate groundstates is a common feature. This would also mean that in this parameter region in our model more magnetically ordered states should be observable. As we are however only able to look for homogeneous or Néel states, this is only speculative, nevertheless motivating further studies of magnetic order in the singleband Hubbard model with different methods to solve the DMFT equations. However, for these studies the Bethe lattice may not be a suitable choice any more, as the definition of a wave vector Q to identify the various possible spin structures is not possible here.