text:
# DDW Order and its Role in the Phase Diagram of Extended Hubbard Models

## I Introduction

The high-$`T_c`$ cuprates exhibit peculiar behavior when underdoped: the density-of-states is depleted at low energies, as if some of the degrees of freedom of the system were developing a gap. This behavior, observed in optical conductivity,optics  NMR,nmr  angle-resolved photoemission, arpes  $`c`$-axis tunneling,tunneling  and specific heat measurements,spec-heat  was dubbed the ‘pseudogap’. The emergence of the ‘pseudogap’ mimics somewhat the impoverishment of the low-energy excitation spectrum which accompanies the development of $`d_{x^2y^2}`$ superconductivity and resembles, more generally, the type of gap formation which is concomitant with a large class of order parameters. However, it does not – at first glance – appear to be connected with the formation of an ordered state. Consequently, it was initially believed that the pseudogap was a crossover phenomenon, and attempts to describe it depended on various approximate methods of treating states with local, fluctuating order.phase-fluct ; gauge

However, it has recently been proposed that the ‘pseudogap’ state *is* actually a broken-symmetry ordered state, and that the signatures of the order are subtle enough that the state was able to appear incognito.Chakravarty01 ; Varma99 ; Kivelson98  In Ref. Chakravarty01, , the $`d_{x^2y^2}`$ density-wave (DDW) state DDW  was advanced as a candidate order. The realization that this is a realistic possibility has led to a re-examination of the experimental circumstances. Recent *elastic* neutron scattering experiments, which directly probe the symmetries broken by DDW order – time-reversal and translation by one lattice spacing – appear to have observed it.Mook  A number of other experiments are consistent with the proposal,Chakravarty01  especially measurements of the superfluid density as a function of doping.Uemura89

The experimental situation seems promising, which is strong incentive to reconsider the theoretical state of affairs. If the ‘pseudogap’ state is, indeed, an ordered state, then we should be able to study it within mean-field theory, as we would study the antiferromagnetic state, superconducting state, or any other ordered state. Mean-field theory is unlikely to explain the detailed shape of the phase boundary, but one can hope that it will capture the broad features of the phase diagram, such as its topology and the basic temperature scales. Deep within any phase, with $`T0`$ and far from any quantum phase transitions, the mean-field Hamiltonian should be the correct Hamiltonian, although the parameters in it may need to be renormalized from their mean-field values. Thus it seems natural to simultaneously study the antiferromagnetic (AF), $`d_{x^2y^2}`$ -wave superconducting (DSC), and DDW order parameters in mean-field theory. The interplay and possible coexistence of these orders should be qualitatively and semi-quantitatively explained by mean-field theory. Phase transitions, quantum or thermal, may not be accurately described in their asymptotic limits, but the AF, DDW, and DSC phases will, as will possible phases with coexisting AF, DDW, and DSC orders.

However, there is an immediate problem faced by such a program. What microscopic Hamiltonian should be used? In the early days of high-$`T_c`$, it was hoped that the important physics of strong local repulsion and superexchange, which is present in the simplest models, such as the Hubbard and $`tJ`$ models, would be sufficient to explain all of the interesting physics of the cuprates. This appears not to be the case. Monte Carlo studies have not found superconductivity in the Hubbard model, MonteCarlo  while Monte Carlo calculations, exact diagonalization, and density-matrix renormalization-group (DMRG) calculations give conflicting results for the $`tJ`$ model.t-J  DMRG studies have found that the behavior of $`n`$-leg ladders depends sensitively on the strength of, for instance, second-neighbor hopping,White99  as have Monte Carlo studies.Poilblanc00  Indeed, some numerical results are sensitively dependent on boundary conditions,Hellberg99  which is further indication of the instability of many of these models to relatively small changes in the parameters. Furthermore, the physics of charge-ordering is probably not correctly described by the $`tJ`$ model without near-neighbor (and possibly long-range) Coulomb repulsion.Kivelson96 ; Zaanen99 ; Arrigoni00  Indeed, it is also clear from experiments that relatively small changes – such as those associated with substituting Nd for La,Tranquada95  which is off the radar screen of the $`tJ`$ and Hubbard models – can radically change at least some aspects of the behavior of these materials. In short, the detailed form of the underlying Hamiltonian matters.

Fortunately, we are not completely in the dark about the nature of the microscopic Hamiltonian. Local Coulomb repulsion, both on-site and near-neighbor, is clearly an important part of the physics. This is known from microscopic calculations of the Hubbard parameters $`t`$, $`U`$, and also from the fact that the undoped parent compounds are antiferromagnetic insulators. The other important clue – which derives entirely from experiments – is that the cuprates superconduct. The correct microscopic model (or models) must support $`d`$-wave superconductivity when doped away from half-filling. If the Hubbard and $`tJ`$ do not – and it appears that they do not for $`t/U`$ small – then they cannot describe the cuprates fully.

Our strategy will be to take a generalization of the Hubbard model which includes next-neighbor repulsion and, most importantly, pair-hopping (or correlated hopping). The pair-hopping term favors superconductivity. Even when it is relatively small, it stabilizes superconductivity in the Hubbard model, as we will see. There are a variety of ways in which such a term – or another term with similar effect – could arise, either from quantum chemistry Appel93 ; Hirsch89  or in the passage to an effective description such as the $`tJ`$ model; in both cases, it is essentially a result of strong local Coulomb repulsion, as superexchange is. In any event, it appears that such physics is necessary to stabilize superconductivity, so we will incorporate it in our model. We will find that such a term also leads to DDW order.

To summarize, we consider a model which is chosen so that it incorporates the basic physics of strong local repulsion *and* so that will have a phase diagram which includes AF at half-filling and DSC at some finite doping. We find that it naturally supports DDW order. In mean-field theory, we find a phase diagram in the temperature-doping plane which resembles the experimental phase diagram of the cuprates, with the DDW phase boundary playing the role of the experimental pseudogap onset line. This DDW line continues into the DSC state, so that the underdoped superconducting state is characterized by both DSC and DDW orders. At low doping, there is also a region of coexistence between AF and DDW orders. We comment on the interpretation of experiments vis-à-vis our findings.

## II Model Hamiltonian

We consider the following bilayer lattice model of interacting electrons: Nayak00a

$$=_{\text{kin}}+_{\text{int}},$$
(1)
where

$$\begin{array}{c}_{\text{kin}}=t_{ij}\underset{i,j}{}\left(c_{i\sigma }^{(\lambda )}c_{j\sigma }^{(\lambda )}+\text{h.c.}\right)\hfill \\ \hfill +\frac{t_{}}{16}\underset{i}{}(c_{i+\widehat{x}+\widehat{y},\sigma }^{(1)}c_{i\sigma }^{(2)}+c_{i+\widehat{x}\widehat{y},\sigma }^{(1)}c_{i\sigma }^{(2)}\\ \hfill c_{i\sigma }^{(1)}c_{i\sigma }^{(2)}c_{i+2\widehat{x},\sigma }^{(1)}c_{i\sigma }^{(2)}+xy+12+\text{h.c.}),\end{array}$$
(2)
and

$$\begin{array}{c}_{\text{int}}=U\underset{i}{}n_i^{(\lambda )}n_j^{(\lambda )}+V\underset{i,j}{}n_i^{(\lambda )}n_j^{(\lambda )}\hfill \\ \hfill t_c\underset{\begin{array}{c}i,j,i^{},j\\ ii^{}\end{array}}{}c_{i\sigma }^{(\lambda )}c_{j\sigma }^{(\lambda )}c_{j\sigma }^{(\lambda )}c_{i^{}\sigma }^{(\lambda )}.\end{array}$$
(3)
In the formulas above $`t_{ij}`$ is hopping with $`t_{ij}=t`$ for nearest neighbors, $`t_{ij}=t^{}`$ for next nearest neighbors and $`t_{ij}=0`$ otherwise. The other parameters are the tunneling $`t_{}`$, the on-site repulsion $`U`$, the nearest-neighbor repulsion $`V`$, and next-nearest-neighbor correlated hopping $`t_c`$. The indices $`i,j`$ correspond to a lattice site, $`\sigma `$ to the spin, and $`\lambda `$ to the layer.

The next-nearest-neighbor correlated hopping term is physically kinetic, but since it is also quartic, we are going to treat it as interaction. It hops an electron from $`i^{}`$ to $`j`$ when $`j`$ is vacated by an electron hopping to $`i`$. These two hops are correlated by virtue of Coulomb interaction between the electrons. The presence of this term in the cuprates has been shown in band-structure calculations.Appel93  Correlated hopping has been discussed in Refs. Wheatley88, ; Hirsch89, ; Chakravarty93, ; Assaad97,  as a possible mechanism of superconductivity, but it has also been foundNayak00a  that it favors DDW order as well.

The tunneling term is momentum conserving.Chakravarty93 ; Andersen95  We consider a CuO2 bilayer because the pseudogap has been best characterized in bilayer materials such as YBCO and Bi2212.

To derive a mean-field theory, it is convenient to take the Fourier transform of Eq. (1) and regroup the terms. This task would be particularly simple if there were only one phase at a given set of parameters. For example, a DDW reduced Hamiltonian would look like

$$_{\text{DDW}}=g_{\text{DDW}}_{k,k^{}}f\left(k\right)f\left(k^{}\right)c_{k+Q,\sigma }^{(\lambda )}c_{k\sigma }^{(\lambda )}c_{k^{}\sigma ^{}}^{(\lambda )}c_{k^{}+Q,\sigma ^{}}^{(\lambda )},$$
(4)
where $`f\left(k\right)=\mathrm{cos}k_x\mathrm{cos}k_y`$ (the lattice spacing has been set to unit) and the DDW mean-field coupling constant is

$$g_{\text{DDW}}=8V+24t_c.$$
(5a)
Similar values of the mean filed coupling constants can be derived for other phases as well. Thus, for antiferromagnesm, $`d`$-wave superconductivity and $`(\pi ,\pi )`$ charge-density wave we derive:
$`g_{\text{AF}}`$ $`=2U,`$ (5b)
$`g_{\text{DSC}}`$ $`=12t_c8V,`$ (5c)
$`g_{\text{CDW}}`$ $`=16V+24t_c2U.`$ (5d)
In fact, the interaction part of the Hamiltonian Eq. (3) can be further generalized to include the *interlayer* Coulomb interactions:

$$_{\text{int}}^{}=U^{}\underset{i}{}n_i^{(\lambda )}n_j^{(\lambda ^{})}+V^{}\underset{i,j}{}n_i^{(\lambda )}n_j^{(\lambda ^{})},$$
(6)
where $`\lambda \lambda ^{}`$. Then for the given interlayer configuration of the order parameters (antisymmetric for AF and DDW and symmetric for DSC), the mean field coupling constants become:

$`g_{\text{DDW}}`$ $`=8V+8V^{}+24t_c,`$ (7a)
$`g_{\text{AF}}`$ $`=2U+2U^{},`$ (7b)
$`g_{\text{DSC}}`$ $`=12t_c8V+8V^{}.`$ (7c)
For the opposite configuration (symmetric for AF and DDW and antisymmetric for DSC) the contributions of $`U^{}`$, $`V^{}`$ would be negative, which is the main reason why such configurations have generally higher energy and are not observed. On the other hand, the fact that five interaction terms produce only three phases means that we can have the same phase diagrams (corresponding to a given set of $`g_p`$’s) for a range of values of the interaction constants. In the following section we will assume that $`U^{}=V^{}=0`$ so that each phase diagram will correspond to a unique set of $`U,V,t_c`$.

The total Hamiltonian contains the reduced parts corresponding to these phases as well as the interactions between the order parameters. However, since we expect $`g_{\text{CDW}}`$ to be negative so that the corresponding order parameter is always zero, we will ignore the term corresponding to this phase. The final form of the reduced Hamiltonian is

$$\begin{array}{c}_{\text{red}}=_kϵ_{k\lambda \lambda ^{}}c_{k\sigma }^{\left(\lambda \right)}c_{k\sigma }^{\left(\lambda ^{}\right)}\hfill \\ \hfill g_{\text{AF}}_{k,k^{}}c_{k+Q,\sigma }^{(\lambda )}c_{k\sigma }^{(\lambda )}c_{k^{}\sigma ^{}}^{(\lambda )}c_{k^{}+Q,\sigma ^{}}^{(\lambda )}\\ \hfill g_{\text{DDW}}_{k,k^{}}f\left(k\right)f\left(k^{}\right)c_{k+Q,\sigma }^{(\lambda )}c_{k\sigma }^{(\lambda )}c_{k^{}\sigma ^{}}^{(\lambda )}c_{k^{}+Q,\sigma ^{}}^{(\lambda )}\\ \hfill g_{\text{DSC}}_{k,k^{}}f\left(k\right)f\left(k^{}\right)c_k^{(\lambda )}c_k^{(\lambda )}c_k^{}^{(\lambda )}c_k^{}^{(\lambda )},\end{array}$$
(8)
where $`ϵ_{k11}=ϵ_{k22}=ϵ_k=2t\left(\mathrm{cos}k_x+\mathrm{cos}k_y\right)4t^{}\mathrm{cos}k_x\mathrm{cos}k_y`$, $`ϵ_{k12}=ϵ_{k21}=ϵ_k=\left(t_{}/4\right)f\left(k\right)^2`$.

The standard Hubbard-Stratonovich mean-field-theoretical treatment of Eq. (8) is to assume the presence of a bosonic mean field, defined as an order parameter, neglect the fluctuations, find the eigenvalues of the Hamiltonian and finally, integrate out the fermion degrees of freedom to derive the free energy.

We define the order parameters of DDW, AF and DSC phases as follows:

$`\varphi _\lambda `$ $`=g_{\text{DDW}}{\displaystyle _k}f\left(k\right)c_{k+Q,\sigma }^{(\lambda )}c_{k\sigma }^{(\lambda )},`$ (9a)
$`M_\lambda `$ $`=g_{\text{AF}}{\displaystyle _k}c_{k+Q,\sigma }^{(\lambda )}c_{k\sigma }^{(\lambda )},`$ (9b)
$`\mathrm{\Delta }_\lambda `$ $`=g_{\text{DSC}}{\displaystyle _k}f\left(k\right)c_k^{(\lambda )}c_k^{(\lambda )}.`$ (9c)
We assume that $`\varphi _\lambda `$ and $`M_\lambda `$ are anti-symmetric in the bilayer index. Then, the free energy of the system is

$$\begin{array}{c}f=\frac{\left|M\right|^2}{g_{\text{AF}}}+\frac{\left|\varphi \right|^2}{g_{\text{DDW}}}+\frac{\left|\mathrm{\Delta }\right|^2}{g_{\text{DSC}}}\hfill \\ \hfill +\underset{s_1,s_2,s_3=\pm 1}{}_{\begin{array}{c}k_x>0\\ k_y>k_x\end{array}}[s_1ϵ_k+s_2ϵ_k\mu \\ \hfill 2T\mathrm{ln}(2\mathrm{cosh}\{\frac{1}{2T}[\left\{f\left(k\right)\mathrm{\Delta }\right\}^2+(s_1\{[ϵ_k+s_2ϵ_k]^2\\ \hfill +[f\left(k\right)\varphi +s_3M]^2\}^{1/2}\mu )^2]^{1/2}\})].\end{array}$$
(10)
As we expand this expression for small values of the order parameters, we can construct a Landau-Ginzburg theory:

$$\begin{array}{c}f\left(T\right)=f_0\left(T\right)+\underset{p}{}a_p\left|\mathrm{\Phi }_p\right|^2+\underset{p}{}b_p\left|\mathrm{\Phi }_p\right|^4\hfill \\ \hfill +\underset{pp^{}}{}c_{pp^{}}\left|\mathrm{\Phi }_p\right|^2\left|\mathrm{\Phi }_p^{}\right|^2,\end{array}$$
(11)
where $`p`$ denotes the kind of the order parameter (AF, DDW or DSC) and $`\mathrm{\Phi }_p`$ is the order parameter ($`M`$, $`\varphi `$ or $`\mathrm{\Delta }`$, respectively). The $`a_p`$ coefficients cross zero at the transitions so that $`a_p=0`$ are the equations that determine critical temperature $`T_c`$:

$$a_p=\frac{1}{g_p}\underset{s_1,s_2,s_3=\pm 1}{}_{\begin{array}{c}k_x>0\\ k_y>k_x\end{array}}K_p\left(k\right),$$
(12)
where

$`K_{\text{AF}}\left(k\right)`$ $`={\displaystyle \frac{1}{2\left|\overline{ϵ}_k\right|}}\mathrm{tanh}\left({\displaystyle \frac{\left|s_1\overline{ϵ}_k\mu \right|}{2T}}\right),`$ (13a)
$`K_{\text{DDW}}\left(k\right)`$ $`=f\left(k\right)^2K_{\text{AF}}\left(k\right),`$ (13b)
$`K_{\text{DSC}}\left(k\right)`$ $`={\displaystyle \frac{f\left(k\right)^2}{2\left|s_1\overline{ϵ}_k\mu \right|}}\mathrm{tanh}\left({\displaystyle \frac{\left|s_1\overline{ϵ}_k\mu \right|}{2T}}\right).`$ (13c)
Here $`\overline{ϵ}_k=ϵ_k+s_2ϵ_k`$. The $`b_p`$ coefficients are positive:

$$b_p=\underset{s_1,s_2,s_3=\pm 1}{}_{\begin{array}{c}k_x>0\\ k_y>k_x\end{array}}K_p^{}\left(k\right),$$
(14)
where

$`K_{\text{AF}}^{}\left(k\right)`$ $`={\displaystyle \frac{\eta _1\left(k\right)}{8\left|\overline{ϵ}_k\right|^3}},`$ (15a)
$`K_{\text{DDW}}^{}\left(k\right)`$ $`=f\left(k\right)^4K_{\text{AF}}^{}\left(k\right),`$ (15b)
$`K_{\text{DSC}}^{}\left(k\right)`$ $`={\displaystyle \frac{f\left(k\right)^4\eta _2\left(k\right)}{8\left|s_1\overline{ϵ}_k\mu \right|^3}},`$ (15c)
where

$`\eta _1\left(k\right)`$ $`=\mathrm{tanh}\left({\displaystyle \frac{\left|s_1\overline{ϵ}_k\mu \right|}{2T}}\right){\displaystyle \frac{\left|\overline{ϵ}_k\right|/2T}{\mathrm{cosh}\left(\frac{\left|s_1\overline{ϵ}_k\mu \right|}{2T}\right)^2}},`$ (16a)
$`\eta _2\left(k\right)`$ $`=\mathrm{tanh}\left({\displaystyle \frac{\left|s_1\overline{ϵ}_k\mu \right|}{2T}}\right){\displaystyle \frac{\left|s_1\overline{ϵ}_k\mu \right|/2T}{\mathrm{cosh}\left(\frac{\left|s_1\overline{ϵ}_k\mu \right|}{2T}\right)^2}}.`$ (16b)
Finally, the $`c_{pp^{}}`$ coefficients are

$$c_{pp^{}}=\underset{s_1,s_2,s_3=\pm 1}{}_{\begin{array}{c}k_x>0\\ k_y>k_x\end{array}}K_{pp^{}}^{\prime \prime }\left(k\right),$$
(17)
where

$`K_{\text{AF,DSC}}^{\prime \prime }\left(k\right)`$ $`={\displaystyle \frac{f\left(k\right)^2}{4\left|s_1\overline{ϵ}_k\mu \right|^2\left|\overline{ϵ}_k\right|}}\eta _2\left(k\right),`$ (18a)
$`K_{\text{DDW,DSC}}^{\prime \prime }`$ $`=f\left(k\right)^2K_{\text{AF,DSC}}^{\prime \prime }\left(k\right),`$ (18b)
$`K_{\text{AF,DDW}}^{\prime \prime }\left(k\right)`$ $`={\displaystyle \frac{3f\left(k\right)^2}{4\left|\overline{ϵ}_k\right|^3}}\eta _1\left(k\right).`$ (18c)
The fact that all $`K_{pp^{}}^{\prime \prime }>0`$ implies that the phases compete with each other.

## III Phase diagram

The mean field phase diagram can be derived by minimizing Eq. (10) at fixed doping.

Since there is a large number of parameters in our model, there is substantial variety in the possible diagrams. One such diagram, generated with $`t=0.5`$ eV, $`t^{}=0.025`$ eV, $`t_{}=0.05`$ eV, $`g_{\text{AF}}=0.06`$ eV, $`g_{\text{DDW}}=0.02`$ eV, $`g_{\text{DSC}}=0.01`$ eV, is shown on the figure 1. The corresponding values of the interaction constants are $`U0.03`$ eV, $`V=0`$ eV, $`t_c0.8\times 10^3`$ eV. Note that for these values of the constants, $`g_{\text{CDW}}=0.04`$ eV $`<0`$, which is consistent with our assumption that a $`(\pi ,\pi )`$ charge-density wave is not energetically favorable.

Another diagram, shown on the figure 2, was generated with $`t=0.5`$ eV, $`t^{}=0`$, $`t_{}=0.1`$ eV, $`g_{\text{AF}}=0.084`$ eV, $`g_{\text{DDW}}=0.038`$ eV, and $`g_{\text{DSC}}=0.017`$ eV. The corresponding values of the interaction constants are $`U0.042`$ eV, $`V1.7\times 10^4`$ eV, $`t_c1.5\times 10^3`$ eV, and also $`g_{\text{CDW}}0.045`$ eV.

As we can see, in both diagrams the antiferromagnetic transition temperature at half-filling is close to 1000 K. This should be understood as the scale at which two-dimensional antiferromagnetic correlations develop locally. Due to the Mermin-Wagner-Coleman theorem, which states that a continuous symmetry cannot be broken spontaneously at finite-temperature in $`2D`$, the transition temperature is zero for a single bilayer. The coupling between different bilayers (which is not included in our single-bilayer calculation) stabilizes the antiferromagnetic phase with a transition temperature around 410 K. In lightly-doped cuprates, the presence of impurities causes the misalignment of locally ordered antiferromagnetic clusters, thereby forming a spin glass. Thus, if we interpret our $`T_N`$ as the scale of local $`2D`$ antiferromagnetic order, which could become $`3D`$ antiferromagnetic order or spin glass order, then the phase diagrams of Figs. 1 and 2 are very reasonable, indeed.

Experiments might lead us to expect that DDW order would occur in the range of doping between 0.07 and 0.19. This range is smaller than one shown on Fig. 2 and a bit larger than that shown in Fig. 1. The temperature scale for this phase on Fig. 1 is very reasonable; it is almost three times higher on Fig. 2. This change occurred primarily as a result of the increased value of $`t_c`$. If we further increase $`t_c`$ to $`1.9\times 10^3`$ eV, the DDW phase will begin to suppress the AF phase and will expand up to half-filling at finite temperatures. In general, varying the interaction constants by less that 20 – 30% does not change the phase diagram qualitatively. However, larger variations lead to completely different classes of phase diagrams, such as those with the AF phase suppressed or without a DDW phase at all. For example, Fig. 3 shows the case when due to the smaller value of correlation hopping, both DDW and DSC phases disappear and only AF phase remains in the diagram.

The DSC phase occupies a doping range away from half-filling primarily as a result of band structure effects associated with the bilayer splitting. In the absence of other orders, it would extend all the way to half-filling, but it is suppressed at low doping by DDW and AF order. In a more realistic calculation, superconductivity would be suppressed close to half-filling by no-double-occupancy constraint, i.e. by strong local Coulomb repulsion. However, the DSC phase never even makes it that close to half-filling because the DDW phase intervenes.

An important feature common to both diagrams is the existence of regions with two simultaneous kinds of order. Namely, there is a region with DDW+AF order and a region with DDW+DSC order. The system is an insulator in the AF state at half-filling, a metal in the DDW and DDW+AF states, and a superconductor in the DSC and DDW+DSC states.

All of the transitions are of second order at the mean-field level because of the signs of the $`c_{pp^{}}`$ couplings between the order parameters in the Ginzburg-Landau theory Eq. (11).

The calculated dependence of the chemical potential $`\mu `$ on the doping $`x`$ inside the DDW phase and in its proximity is nonmonotonic. This is due to the rapid development of the DDW gap, which causes the chemical potential to be lower than in the normal state. The thermodynamic inequality $`(\mu /x)_{T,V}0`$ implies that when this is violated, mean-field theory should be corrected using Maxwell’s construction, which signals that fluctuations drive the transitions first-order as a function of $`\mu `$. Consequently, we would expect the underdoped side of the DSC phase to be characterized by a smaller than expected chemical shift, as has been observed.Ino97  A first-order phase transition as a function of chemical potential is manifested as phase separation in a two-phase coexistence region spanning a range of dopings when the doping is held fixed instead. It has been argued that such phase separation will be precluded by Coulomb interactions, thereby leading to stripe formation.Kivelson96 ; Zaanen99

## IV Conclusion

We have studied the phase diagram of a bilayer lattice model using mean field theory. Since we have focused on ordered phases, this should be a valid approximation. We found that for certain ranges of the values of the interaction constants the phase diagram agrees well with the experimentally observed phase diagram of YBCO if the ‘pseudogap’ is associated with DDW order. The diagram remains in qualitative agreement with the experimental data when the parameters of our model vary by less than 20 – 30% and becomes qualitatively different for larger variations. Clearly, such a phase diagram is reasonably robust, but is hardly inevitable. This is reassuring because high-temperature superconductivity is stable, but only appears in a special class of materials (to the best of our current knowledge).

There are some systematic errors associated with mean-field theory, on which we now comment. It underestimates the effect of fluctuations. Thus, the Néel temperature is very large in mean-field theory, while it should actually be zero in any strictly two-dimensional system. However, the Néel temperature which we find should be regarded as the temperature below which a renormalized classical description is valid.Chakravarty89  The Néel temperature observed in experiments is associated with the crossover from 2D to 3D. Mean-field theory also overestimates the coupling which drives antiferromagnetism, which it takes to be essentially $`U`$. For small $`U`$, this is correct, but for large $`U`$, it should be replaced by $`Jt^2/U`$. Indeed, the large-$`U`$ limit is generally somewhat problematic near half-filling since the Gutwiller constraint is not enforced in mean-field theory. The $`d_{x^2y^2}`$ symmetry of the DDW and DSC states lead one to the erroneous conclusion that they are completely unaffected by large $`U`$. This cannot, of course, really be true; clearly, mean-field theory underestimates the tendency of large-$`U`$ to push these ordered states away from half-filling. The seemingly small value of $`U`$ taken in our calculation should be interpreted in light of these observations. Other mean-field treatments which incoporate strong local Coulomb repulsion more prominently have also found DDW order in a generalization of the $`tJ`$ model Cappelluti99  and in the Hubbard model with nearest neighbor attraction.Stanescu01

We find that the scale associated with superconductivity is largely determined by the strength of correlated hopping. At the moment, this is rather *ad hoc*, but we had little choice but to introduce some term of this sort in order to have a phase diagram which includes superconductivity. It is possible that the superexchange coupling $`J`$ plays a more important role than we have accorded it in setting $`T_c`$, but superexchange is beyond a mean-field treatment.

As we have seen, the very term which stabilizes superconductivity also supports the development of DDW order. One way of interpreting our results begins with the observation that the DDW order parameter, when combined with the real and imaginary parts of the DSC order parameter form a triplet under an $`SU(2)`$ group of transformations.Nayak00a ; Nayak00b  If this ‘pseudospin’ $`SU(2`$) is a symmetry of the Hamiltonian, then DDW and DSC orders will be equally favored. Thus, one can envision that the important order-producing term in the Hamiltonian is $`SU(2)`$-symmetric while small symmetry-breaking terms drive the system into either the DDW or DSC states. Our result shows that pair-hopping is of this form. Are all physically reasonable mechanisms for $`d_{x^2y^2}`$ superconductivity similarly invariant under pseudospin $`SU(2)`$? This is an open problem; we have answered in the affirmative for one particular class of Hamiltonians.

###### Acknowledgements.

We would like to thank S. Chakravarty for discussions. C.N. is supported by the National Science Foundation under Grant No. DMR-9983544 and by the Alfred P. Sloan Foundation. E.P. is supported in part by the Department of Energy under grant DE-FG03-92-ER 40701.


format: text

text:
# Cosmological Parameters from Redshift-Space Correlations

## 1 Introduction

Recently, independent cosmological observations have turned out to be consistent with a canonical model of the structure formation based on the assumption that the main ingredients of our universe are the cold dark matter (CDM), and the dark energy component such as a cosmological constant. Introducing the unknown dark matter and the peculiar dark energy could have been thought of as somewhat artificial at the time when observations were not constraining enough, now this model has become more and more convincing with remarkable observational support.

To make a cosmological model truly convincing, it is important to make sure that any possible observations are consistent with that model. This can be accomplished by determining cosmological parameters as precisely as possible from each observation, and then comparing the results with one another. So far the individual observations were only able to constrain the parameters into some region in a multi-dimensional parameter space. Therefore, it is quite common to combine different observations to obtain the tight constraint on the parameters. A striking example is a combination of fluctuations in cosmic microwave background (de Bernardis et al., 2000) and a Hubble diagram of supernova Ia (Riess et al., 1998; Perlmutter et al., 1999) which is considered to be the prominent evidence of the current canonical model.

However, one can test the consistency of the theory only when each observation can independently constrain the parameters in a sufficiently small region. In this paper, we explore the capabilities of redshift surveys in this respect. In redshift surveys, there are several degrees of freedom in planning survey strategies: the area and the depth of a survey, and the choice of specific targets for the spectroscopic measurements. Each of these affects the analysis in a different and complex fashion. Thus, it is important to estimate the impact of a given strategy on the cosmological analysis of the data.

In this paper, we apply our considerations to the Sloan Digital Sky Survey (SDSS, York et al., 2000). It is not only the largest ongoing redshift survey, but it also contains three different types of spectroscopic targets within the same area of the sky that can be used to infer cosmological information.

In deep, wide redshift surveys like SDSS, the apparent clustering properties in observable redshift space are complicated because of the redshift-evolution and redshift-distortion effects. In linear regime, with a standard initial Gaussian density field, the correlation function $`\xi (r)`$ or the power spectrum $`P(k)`$ completely characterizes the statistical property of the density field. Spatial homogeneity and isotropy in real space guarantees these to be one-dimensional functions of the scale, $`r`$ or $`k`$. So it is usual to estimate the real-space power spectrum with certain correction methods or ansatz from observed redshift-space data, in order to compare with the theoretical prediction of various cosmological models (e.g., Davis & Peebles, 1983; Feldman, Kaiser & Peacock, 1994).

However, the redshift-evolution and redshift-distortion also have invaluable information on cosmology. The redshift-evolution of the density fluctuations depends on the growth factor which is a function of the density parameter $`\mathrm{\Omega }_\mathrm{M}`$ and the normalized cosmological constant parameter $`\mathrm{\Omega }_\mathrm{\Lambda }`$. It is also a function of the galaxy bias $`b`$ which contains the information on the unknown galaxy formation process. The redshift-distortion depends on the density-velocity relation which is a function of $`\mathrm{\Omega }_\mathrm{M}^{0.6}/b`$ (Kaiser, 1987), and also on $`\mathrm{\Omega }_\mathrm{M}`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ through geometrical distortions (Alcock & Paczyński, 1979; Ballinger, Peacock & Heavens, 1996; Matsubara & Suto, 1996).

When our main purpose is not only determining the power spectrum of initial fluctuations, but also constraining above parameters, it is desirable to directly analyze the data to obtain cosmological information without any intermediate statistical quantities like the power spectrum. Since the apparent clustering in redshift space is quite different from that in real space, we can take advantage of the maximum-likelihood method to directly constrain cosmological models from the apparent clustering data in redshift space.

In our previous work (Matsubara & Szalay, 2001), we focused on how near-future redshift surveys can constrain the cosmological constant by the maximum-likelihood method. Since the redshift surveys become larger and larger, even the multi-dimensional parameter space can be constrained by redshift surveys alone. Combining with other cosmological observations like the cosmic microwave background fluctuations, Hubble diagram of the type Ia supernova, and so on, one can constrain the parameters more accurately. However, independent observational constraints can serve as a consistency check of the fundamental cosmological model. In this respect, we investigate how one can constrain the multi-dimensional parameter space of the standard CDM model with cosmological constant from redshift surveys, extending our previous work.

## 2 Maximum-likelihood Method for the Apparent Redshift-space Clustering

In maximum-likelihood methods, we directly deal with the density fluctuations in redshift space. In order to do this, first we place smoothing cells in redshift space and consider the number of objects $`n_i`$ in each cell $`i`$. Each count represents the smoothed density field contaminated by shot noise. For given mean number densities, $`\overline{n}_i=n_i`$, the density fluctuations $`a_i=n_i/\overline{n}_i1`$ are the fundamental quantities to be analyzed. In the following, the most important statistical quantity is the correlation matrix, $`C_{ij}=a_ia_j`$. This matrix is related to a continuous two-point correlation function in redshift space $`\xi (𝒔,𝒔^{})=\delta (𝒔)\delta (𝒔^{})`$, where $`\delta (𝒔)`$ is the underlying continuous density contrast which survey objects trace. The relation between the correlation matrix and the correlation function is given by

$`C_{ij}={\displaystyle d^3sd^3s^{}W^{(\mathrm{s})}(𝒔_i𝒔)W^{(\mathrm{s})}(𝒔_j𝒔^{})\xi (𝒔,𝒔^{})}+{\displaystyle d^3s\frac{W^{(\mathrm{s})}(𝒔_i𝒔)W^{(\mathrm{s})}(𝒔_j𝒔)}{\overline{n}(𝒔)}},`$ (2.1)
where $`W^{(\mathrm{s})}(𝒔)`$ is the smoothing kernel in redshift space, normalized as $`d^3sW^{(\mathrm{s})}(𝒔)=1`$ and $`\overline{n}(𝒔)`$ is the mean number density. The last term represents contribution from the shot noise effect. In linear regime, the two-point correlation function in redshift space is analytically given, including high redshift effect.

Although the general form is fairly complicated (Matsubara, 2000), the distant-observer approximation (Kaiser, 1987; Matsubara & Suto, 1996), which is applicable in many cases, simplify the analysis below. The corresponding formula (Matsubara & Suto, 1996; Matsubara, 2000) is given by

$`\xi (𝒔_1,𝒔_2)=b(z_1)b(z_2)D(z_1)D(z_2){\displaystyle \frac{d^3k}{(2\pi )^3}e^{i𝒌(𝒙_1𝒙_2)}\left[1+\beta (z_1)\frac{k_{}^{}{}_{}{}^{2}}{k^2}\right]\left[1+\beta (z_2)\frac{k_{}^{}{}_{}{}^{2}}{k^2}\right]P(k)},`$ (2.2)
where $`z_1`$ and $`z_2`$ are the redshifts of the two points, $`k_{}`$ is the line-of-sight component of the wave number, $`P(k)`$ is the linear mass power spectrum at $`z=0`$, $`D(z)`$ is the linear growth rate normalized as $`D(0)=1`$, $`b(z)`$ is the bias factor at redshift $`z`$, and $`\beta (z)`$ is the redshift distortion parameter, which is accurately fitted by the redshift-dependent mass density parameter $`\mathrm{\Omega }_\mathrm{M}(z)`$ and normalized cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }`$ as $`\beta (z)\{\mathrm{\Omega }_\mathrm{M}^{4/7}(z)+\mathrm{\Omega }_\mathrm{\Lambda }(z)[1+\mathrm{\Omega }_\mathrm{M}(z)/2]/70\}/b(z)`$ (Lahav et al., 1991; Lightman & Schechter, 1990). The vectors $`𝒙_1`$ and $`𝒙_2`$ are the comoving positions of the two points which are labeled by $`𝒔_1`$ and $`𝒔_2`$ in redshift space (see Matsubara & Suto, 1996).

Except for the nearby universe, the comoving distance is not proportional to the redshift so that the smoothing kernel $`W^{(\mathrm{s})}`$ in redshift space is distorted along the line of sight. In a distant-observer regime, a line-of-sight component $`s_{}`$ and a transverse component $`s_{}`$ of redshift-space separation $`𝒔`$ are given by

$`s_{}=H(z)x_{},s_{}={\displaystyle \frac{z}{d_A(z)}}x_{},`$ (2.3)
where $`x_{}`$ and $`x_{}`$ are corresponding line-of-sight and transverse components of real-space separation, respectively, and $`z`$ is the mean redshift. We have chosen the unit system $`H_0=1`$ here for simplicity. This unit system simplifies the equations in this section. We mention that another unit system for length is introduced in §4 for a practical reason which is explained there. In the above transform, the redshift-dependent Hubble’s parameter and the comoving angular diameter distance are explicitly given by

$`H(z)=\sqrt{\mathrm{\Omega }_\mathrm{M}(1+z)^3+(1\mathrm{\Omega }_\mathrm{M}\mathrm{\Omega }_\mathrm{\Lambda })(1+z)^2+\mathrm{\Omega }_\mathrm{\Lambda }},`$ (2.4)
$`d_A(z)=\{\begin{array}{cc}(K)^{1/2}\mathrm{sinh}\left[(K)^{1/2}\chi (z)\right]\hfill & (K<0),\hfill \\ \chi (z)\hfill & (K=0),\hfill \\ K^{1/2}\mathrm{sin}\left[K^{1/2}\chi (z)\right]\hfill & (K>0),\hfill \end{array}`$ (2.8)
where

$`K=1\mathrm{\Omega }_\mathrm{M}\mathrm{\Omega }_\mathrm{\Lambda },\chi (z)={\displaystyle _0^z}{\displaystyle \frac{dz^{}}{H(z^{})}},`$ (2.9)
are the curvature of the universe and comoving distance to the redshift $`z`$. Thus, the transform from the comoving space to the redshift space is the function of $`\mathrm{\Omega }_\mathrm{M}`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$

As a result of the anisotropic transform of equation (2.3), the spherical kernel in redshift space, for example, ends up with oscillating two-dimensional numerical integration (Matsubara & Szalay, 2001) in evaluating the correlation matrix $`C_{ij}`$. Unfortunately, such numerical integration increases CPU time and decreases accuracy. This is crucial for the analysis of the real data, because the position of kernels can hardly be regular so that one has to calculate the correlation matrix pair by pair. The shape of the kernel is not important for largely separated pairs, but it affects nearby and identical pairs.

Although one can not exactly know what is the distortion unless the parameters $`\mathrm{\Omega }_\mathrm{M}`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ are fixed, it turns out to be a good approximation to set an approximate spherical kernel in real space with approximate values of these parameters. In fact, we found that the latter approximation is better than doing the exact numerical integration which inevitably has integration errors. When the parameters $`\mathrm{\Omega }_\mathrm{M}`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ should also be determined in real analyses, one can iterate the analysis to use better estimates for those parameters, but approximate values are good enough. Thus, suppose that we can approximately set a spherical kernel in real space, so that we have $`W^{(\mathrm{s})}(𝒔)=J^1W_R(|𝒙|)`$, where $`J=(𝒔)/(𝒙)`$ is the Jacobian from $`x`$-space to $`s`$-space, and $`R`$ is the (effective) smoothing radius. In the distant-observer approximation, $`J(z)=z^2H(z)/d_{A}^{}{}_{}{}^{2}(z)`$.

This approximation dramatically simplifies the evaluation of the correlation matrix. In fact, when the kernels are spherical in real space, at least approximately, equations (2.1) and (2.2) give

$`C_{ij}=b_ib_jD_iD_j{\displaystyle \frac{d^3k}{(2\pi )^3}e^{i𝒌(𝒙_i𝒙_j)}\left[1+\beta _i\frac{k_{}^{}{}_{}{}^{2}}{k^2}\right]\left[1+\beta _j\frac{k_{}^{}{}_{}{}^{2}}{k^2}\right]W^2(kR)P(k)}`$
$`+{\displaystyle \frac{1}{\sqrt{J_iJ_j\overline{n}_i\overline{n}_j}}}{\displaystyle \frac{d^3k}{(2\pi )^3}e^{i𝒌(𝒙_i𝒙_j)}W^2(kR)},`$ (2.10)
where $`b_i=b(z_i)`$, $`D_i=D(z_i)`$, $`J(z_i)=J_i`$, etc., and $`W(kR)`$ is the three-dimensional Fourier transform of the spherical kernel $`W_R(|𝒙|)`$. In our notation, $`\overline{n}_i`$ is the mean number density in redshift space so that $`J_i\overline{n}_i`$ is the mean number density in real space. We assume that the variation of the selection function within a kernel is negligible. The above equation means that the distant-observer redshift-distortion is commutable with spherical smoothing in real space, comparing with equation (2.2) besides shot-noise term. Angular integration of the first term can be analytically performed (see Hamilton, 1992; Matsubara & Suto, 1996; Matsubara, 2000), resulting in

$`C_{ij}=b_ib_jD_iD_j\{[1+{\displaystyle \frac{1}{3}}(\beta _i+\beta _j)+{\displaystyle \frac{1}{5}}\beta _i\beta _j]P_0(\mu _{ij})\xi _0(x_{ij},R)`$
$`[{\displaystyle \frac{2}{3}}(\beta _i+\beta _j)+{\displaystyle \frac{4}{7}}\beta _i\beta _j]P_2(\mu _{ij})\xi _2(x_{ij},R)+{\displaystyle \frac{8}{35}}\beta _i\beta _jP_4(\mu _{ij})\xi _4(x_{ij},R)\}`$
$`+{\displaystyle \frac{1}{\sqrt{J_iJ_j\overline{n}_i\overline{n}_j}}}{\displaystyle \frac{k^2dk}{2\pi ^2}j_0(kx_{ij})W^2(kR)},`$ (2.11)
where $`\beta _i=\beta (z_i),\beta _j=\beta (z_j)`$, $`x_{ij}=|𝒙_i𝒙_j|`$ is the comoving separation between the centers of the kernels, $`\mu _{ij}=x_{ij}/x_{ij}`$ is the direction-cosine of the separation $`x_{ij}`$ relative to the line of sight, $`j_l`$ is the spherical Bessel function, and $`P_l(\mu )`$ is the Legendre polynomial. We also define the quantity,

$`\xi _l(x,R)={\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{k^2dk}{2\pi ^2}}j_l(kx)W^2(kR)P(k).`$ (2.12)
We need to numerically integrate the equation (2.12). Once we make a table of the integrated results for various $`x`$ and $`R`$ of interest, we can use this table to interpolate $`\xi _l`$, which makes the evaluation of the correlation matrix $`C_{ij}`$ a really fast operation. This method is very effective because the number of elements of a correlation matrix is the square of the number of cell sites $`(N_{\mathrm{cell}})^2`$, where $`N_{\mathrm{cell}}1,000`$$`10,000`$, or more. We tried exact 2-dimensional integrations for identical and nearby cell-pairs, and found that they agree with the above approximation within the integration errors, even though the assumed $`\mathrm{\Omega }_\mathrm{M}`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ are different from the true values by a factor of about 0.5.

The shot noise term is diagonal when the kernels do not overlap each other, as one can see from the equation (2.1). Using a finite-volume kernel without any overlapping region, the last term of equation (2.1) or (2.11) is given by

$`{\displaystyle \frac{\delta _{ij}}{\overline{n}_i}}{\displaystyle d^3s\left[W^{(\mathrm{s})}(𝒔)\right]^2}={\displaystyle \frac{\delta _{ij}}{J_i\overline{n}_i}}{\displaystyle \frac{k^2dk}{2\pi ^2}W^2(kR)}.`$ (2.13)
The top-hat kernel, $`W_R(x)=3\mathrm{\Theta }(Rx)/(4\pi R^3)`$, where $`\mathrm{\Theta }`$ is the step function, is the most popular example of the finite-volume kernel, in which case the above shot noise term simply reduces to $`\delta _{ij}/\overline{N}_i`$, where $`\overline{N}_i=4\pi R^3J_i\overline{n}_i/3`$ is the expected mean number count in a cell.

Although we only use the top-hat kernel in this paper, it is also useful to consider another possibility as digression, since this fast numerical approximations will also be used in the analysis of the real data. The top-hat kernel is not a smooth function, thus it is sometimes desirable to have an smoother kernel. An example of the smoother, yet finite-volume kernels is the Epanechnikov kernel, which is defined by $`W_R(x)=15(1x^2/R^2)\mathrm{\Theta }(Rx)/(8\pi R^3)`$. While the top-hat kernel is discontinuous on the edge $`x=R`$, the Epanechnikov kernel is continuous. The derivative is still discontinuous on the edge. One can also use the $`m`$-weight Epanechnikov kernel defined by

$`K_m(x;R)={\displaystyle \frac{(2m+3)!!}{m!2^{m+2}\pi R^3}}\left(1{\displaystyle \frac{x^2}{R^2}}\right)^m\mathrm{\Theta }(Rx).`$ (2.14)
The top-hat kernel and the Epanechnikov kernel are the special cases of this general Epanechnikov kernel with $`m=0`$ and $`m=1`$, respectively. The $`(m1)`$-derivative of the $`m`$-weight Epanechnikov kernel is continuous on the edge. The larger the rank in the general Epanechnikov kernel is, the more regular it is on the edge. However, the effective volume of the kernel is smaller for a larger-weight kernel when a smoothing length $`R`$ is fixed. Taking limits $`m\mathrm{},R\mathrm{}`$ with $`m^{1/2}R`$ fixed, one can show that the higher-weight Epanechnikov kernel reduces to a Gaussian kernel with a variance $`R/\sqrt{2m}`$. The Fourier window function for the $`m`$-weight Epanechnikov kernel has a simple form,

$`W_m(kR){\displaystyle _0^{\mathrm{}}}d^3xe^{i𝒌𝒙}K_m(x;R)={\displaystyle \frac{(2m+3)!!j_{m+1}(kR)}{(kR)^{m+1}}}.`$ (2.15)
The larger the rank $`m`$ is, the faster the window function drops off with the wavenumber $`k`$, which reflects the smoothness on the edge of the kernel. It is desirable to have a degree of freedom to choose an appropriate rank $`m`$, because one can optimize the analysis depending on the nature of an individual data set.

The observational data is the set of discrete density fluctuations $`a_i`$. The cosmological models are characterized by a given set of parameters $`(\theta _\alpha )`$. Besides $`\mathrm{\Omega }_\mathrm{M}`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, and $`b`$, which is explicitly appeared in the expression of the correlation matrix, the power spectrum contains several other parameters when we characterize it by a standard cold dark matter model. The latter parameters are the normalization $`\sigma _8`$, baryon mass fraction $`f_\mathrm{B}\mathrm{\Omega }_\mathrm{B}/\mathrm{\Omega }_\mathrm{M}`$, spectral index $`n`$, and Hubble constant $`H_0`$, etc.

According to the Bayes’ rule, the direct likelihood function of the cosmological parameters for a data set is given by

$`(𝜽|𝒂)={\displaystyle \frac{(𝒂|𝜽)(𝜽)}{(𝒂)}},`$ (2.16)
where $`(𝜽|𝒂)`$ indicates the conditional likelihood function of the parameter set $`𝜽=(\theta _1,\theta _2,\mathrm{},\theta _M)`$ given the date set $`𝒂=(a_1,a_2,\mathrm{},a_N)`$, $`(𝜽)`$ indicates the prior likelihood function of the parameter set $`𝜽`$, and so on. Normalization determines the denominator of the right hand side, and the prior likelihood $`(𝜽)`$ comes from the other observations, or just a constant if we do not combine our prier knowledge on the cosmological parameters. The remaining factor $`(𝒂|𝜽)`$ is the distribution function of the data set given a cosmological model. In linear regime where the density fluctuations are considered to be random Gaussian, this function is simply the multi-variate Gaussian distribution function:

$`(𝒂|𝜽)={\displaystyle \frac{1}{\sqrt{(2\pi )^Ndet𝑪(𝜽)}}}\mathrm{exp}\left({\displaystyle \frac{1}{2}}𝒂^\mathrm{T}𝑪^1𝒂\right),`$ (2.17)
where $`𝑪^1`$ is the inverse matrix of the correlation matrix $`𝑪=(C_{ij})`$. The matrix $`𝑪`$ is the function of cosmological parameters $`𝜽`$ as noted above. Thus, calculation of the likelihood function is reduced to the calculation of equation (2.11) for a given survey geometry and selection function for various cosmological models.

## 3 Evaluation of the Fisher Information Matrix

The Fisher information matrix,

$`F_{\alpha \beta }(𝜽)={\displaystyle \frac{\mathrm{ln}(𝒂|𝜽)}{\theta _\alpha }}{\displaystyle \frac{\mathrm{ln}(𝒂|𝜽)}{\theta _\beta }}={\displaystyle \frac{^2\mathrm{ln}(𝒂|𝜽)}{\theta _\alpha \theta _\beta }},`$ (3.1)
where $`\mathrm{}=(𝒂|𝜽)\mathrm{}d^Na`$ represents the averaging over the possible data realization given the model parameters, is the key quantity in estimation theory, because the Cramér-Rao theorem (Kendall & Stuart, 1969; Therrien, 1992) states that an unbiased estimator constrains the model parameters with a minimum variance

$`\mathrm{\Delta }\theta _\alpha \mathrm{\Delta }\theta _\beta \text{}\text{}\left(F^1\right)_{\alpha \beta },`$ (3.2)
where $`F^1`$ is the inverse matrix of the Fisher matrix and $`\mathrm{\Delta }\theta _\alpha =\widehat{\theta }_\alpha \theta _\alpha `$ is the deviation of the estimated value $`\widehat{\theta }_\alpha `$ from the “true” value $`\theta _\alpha `$. In the above equation, the notation $`\text{}\text{}`$ means that the matrix whose $`\alpha \beta `$-elements are given by $`\mathrm{\Delta }\theta _\alpha \mathrm{\Delta }\theta _\beta \left(F^1\right)_{\alpha \beta }`$ is positive semidefinite. Moreover, the equation (3.2) is satisfied with equality if one uses the maximum likelihood estimate, which is obtained by a maximization of the likelihood function $`(𝒂|𝜽)`$ as in the previous section:

$`\mathrm{\Delta }\theta _\alpha \mathrm{\Delta }\theta _\beta =\left(F^1\right)_{\alpha \beta }\text{(maximum likelihood estimate)}.`$ (3.3)
Thus, the inverse of the Fisher matrix gives us the expected constraint on the model parameters $`𝜽`$. A fiducial model for the “true” values of $`𝜽`$ are required to calculate the Fisher matrix.

In the linear regime where the distribution function is given by equation (2.17), the Fisher matrix is analytically reduced to the following form (see, e.g., Vogeley & Szalay, 1996; Tegmark, Taylor & Heavens, 1997):

$`F_{\alpha \beta }={\displaystyle \frac{1}{2}}\mathrm{Tr}\left(𝑪^1{\displaystyle \frac{𝑪}{\theta _\alpha }}𝑪^1{\displaystyle \frac{𝑪}{\theta _\beta }}\right)`$ (3.4)
Using this equation, the calculation of the Fisher matrix is straightforward once the calculation of the correlation matrix for a given model is established in the previous section.

The ellipsoid defined by an equation $`\mathrm{\Delta }𝜽^\mathrm{T}𝑭\mathrm{\Delta }𝜽=\nu ^2`$ defines the concentration ellipsoid which is interpreted as error bounds for a given observation. If the parameter likelihood function $`(𝜽|𝒂)`$ is Gaussian, the concentration ellipsoid with a fixed value of $`\nu `$ represents the $`\nu \sigma `$ confidence level when one uses the maximum likelihood estimate. Although in general the parameter likelihood function is not Gaussian, the concentration ellipsoids still give the rough idea of expected error bounds with a given confidence level of $`\nu `$.

When the parameter space is just two-dimensional, one can easily plot the concentration ellipses. The presentation of higher-dimensional space is not easy. We have at most seven parameters in our model: $`(\mathrm{\Omega }_\mathrm{M},\mathrm{\Omega }_\mathrm{\Lambda },b,f_\mathrm{B},h,\sigma _8,n)`$. To depict the higher dimensional concentration ellipsoids, we can plot all the possible two-dimensional parameter spaces which are marginalized over all the other parameters. This gives the idea how one can impose the constraints in the multi-dimensional parameter space. Although the form of the parameter likelihood function is not known, we can assume it is approximated by a multi-variate Gaussian function so that $`(𝜽|𝒂)\mathrm{exp}(\frac{1}{2}\mathrm{\Delta }𝜽^\mathrm{T}𝑭\mathrm{\Delta }𝜽)`$. In which case, the marginalization over the other parameters but two can be performed. For example, from the standard Gaussian integration, marginalized 2-dimensional likelihood function for two parameters $`\theta _\alpha `$ and $`\theta _\beta `$ is given by

$`(\theta _\alpha ,\theta _\beta |𝒂)\mathrm{exp}\left\{{\displaystyle \frac{1}{2}}\left(\mathrm{\Delta }\theta _\alpha \mathrm{\Delta }\theta _\beta \right)\left[\begin{array}{cc}\left(F^1\right)_{\alpha \alpha }& \left(F^1\right)_{\alpha \beta }\\ \left(F^1\right)_{\alpha \beta }& \left(F^1\right)_{\beta \beta }\end{array}\right]^1\left(\begin{array}{c}\mathrm{\Delta }\theta _\alpha \\ \mathrm{\Delta }\theta _\beta \end{array}\right)\right\},`$ (3.9)
up to normalization factor, where $`(F^1)_{\alpha \beta }`$ indicates the $`\alpha \beta `$-element of the inverse matrix of the original higher-dimensional Fisher matrix, and so on. Therefore, plotting the concentration ellipses in the marginalized parameter space is straightforward: they are given by contours of the equation (3.9). Generally, concentration ellipsoids in the marginalized parameter sub-space are deduced from the sub-matrix of the full inverse Fisher matrix.

## 4 Results for the SDSS geometry

The largest redshift map which will be available in next several years is the SDSS survey. Therefore, it is of great interest how the maximum-likelihood method described in §2 can constrain various cosmological parameters. We have seen in §3 that it is straightforward to estimate the expected constraints on parameters given the survey geometry. We examine the maximum-likelihood method with the SDSS geometry as a typical application of our method.

Redshift-space observations in the SDSS survey consist of three samples of objects (York et al., 2000). The first sample is the set of approximately 1,000,000 “Main Galaxies”, which are selected by limiting an $`r^{}`$-band magnitude. The redshift of main galaxies extends up to $`z0.2`$. The second sample is the set of approximately 120,000 “Luminous Red Galaxies” (LRG’s), which are selected on the basis of color and magnitude (Eisenstein et al., 2001). Because of the clever selection criteria, the LRG’s are essentially a distance-limited sample whose redshift extends up to $`z0.45`$. The third sample is the set of approximately 100,000 quasars, which are selected from color-color diagrams. The redshift of the quasar sample extends up to $`z3`$.

In evaluating Fisher information matrices, the precise geometry of the survey is not needed, and one can approximate the survey geometry to reduce the numerically intense procedure of inverting huge matrices. We setup a generic cubic box with 1,000 cells for each samples. The redshift of the box is set to be a “mean” redshift of each sample: $`z=0.1`$ for main galaxies, $`z=0.3`$ for LRG’s, and $`z=1.8`$ for quasars. The precise value of this mean redshift is not a decisive factor in the Fisher matrix evaluations. The finiteness of the spatial volume and the number density of the objects are the main source of the statistical uncertainty, i.e., cosmic variance and shot noise. From the sky coverage of the three samples of the SDSS spectroscopic survey, which is about pi-steradian, and from the redshift distribution of each samples, one can obtain expected spatial volumes and number densities for each sample.

Although the unit system $`H_0=1`$ used in the previous sections simplifies the equations, it is not convenient in mentioning the clustering scales which become small figures in this system. Thus, we use a different unit system when we have to mention the length scale in the following paragraphs. Following Matsubara & Szalay (2001), we adopt a coordinate system in which radial distance $`s`$ to the objects is measured by redshift times Hubble length: $`s=cz/H_0=(2997.9h^1\mathrm{Mpc})\times z`$, directly observable quantity, instead of the comoving distance. The latter depends on cosmological models through parameters $`\mathrm{\Omega }_\mathrm{M}`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$, which we would like to measure. We introduce the factor $`c/H_0`$ in the definition of $`s`$ simply in order to avoid the small figures when we mention 10-100 Mpc scales. This measure of the radial distance $`s`$ is a linear extrapolation of the distance-redshift relation of the nearby universe. We save the unit $`h^1\mathrm{Mpc}`$ for the usual comoving coordinates, and we use a new notation $`h^1\mathrm{Mpc}_z`$ for the measure $`s`$. The relation of the distance $`s`$ and the comoving distance $`x`$ is redshift-dependent and is given by a differential form, $`ds=H(z)/H_0dx`$. In the unit system $`c=1`$ and $`H_0=1`$ which we used in previous sections, the distance $`s`$ is simply the redshift $`z`$.

We put the 1,000 cells in a generic cubic box of $`L^3`$ on regular $`10\times 10\times 10`$ sites in redshift space. The top-hat kernel \[$`m=0`$ in equation (2.14)\] is used and the radius $`R`$ of the kernel is just a half of the cell separation: $`R=L/20`$, so that there are no overlapping regions to ensure the independence of the cell volumes. The Fisher matrix is straightforwardly calculated for this generic cubic region. This generic box is smaller than the SDSS survey volume, therefore the Fisher matrix can be rescaled as if the whole survey volume consisted of many independent generic cubic boxes. One can easily see that the Fisher matrix of the composite sample is just given by a sum of Fisher matrices of the individual samples. Although we will miss the information from correlations between the boxes, the signals from the correlation at such large separations are very small.

The box size $`L`$ is chosen so that the cell radius be large enough to track the linear regime and not to have large shot noise. On the other hand, too large $`L`$ makes the cosmic variance large. Therefore, we should take a smallest $`L`$ which satisfies the above two conditions. We take $`L=200h^1\mathrm{Mpc}_z`$ for main galaxies and LRG’s, and $`L=1000h^1\mathrm{Mpc}_z`$ for quasars. As a result, we have the number densities $`35/(20h^1\mathrm{Mpc}_z)^3`$ for main galaxies, $`0.5/(20h^1\mathrm{Mpc}_z)^3`$, for LRG’s, and $`0.12/(100h^1\mathrm{Mpc}_z)^3`$ for QSO’s.

We assume fiducial values of model parameters such as $`(\mathrm{\Omega }_\mathrm{M},\mathrm{\Omega }_\mathrm{\Lambda },\mathrm{\Omega }_\mathrm{B}/\mathrm{\Omega }_\mathrm{M},h,\sigma _8,n)=(0.3,0.7,0.13,0.7,1,1)`$. The bias parameters are expected to be different among each samples. We take fiducial values $`b=1,2,2`$ for main galaxies, LRG’s and quasars, respectively. For an evaluation of the power spectrum, we use an useful fitting formula<sup>1</sup><sup>1</sup>1http://background.uchicago.edu/~whu/transfer/transferpage.html given by Eisenstein & Hu (1998) instead of running a Boltzmann code such as the publicly available CMBfast<sup>2</sup><sup>2</sup>2http://physics.nyu.edu/matiasz/CMBFAST/cmbfast.html (Seljak & Zaldarriaga, 1996).

In Figure 1 we plot the expected errors from the fiducial model, when the mass density parameter $`\mathrm{\Omega }_\mathrm{M}`$ and the cosmological constant parameter $`\mathrm{\Omega }_\mathrm{\Lambda }`$ are simultaneously evaluated.

Three contours for each sample represent the concentration ellipses of $`\nu =1,2,3`$ which correspond to expected $`1\sigma `$, $`2\sigma `$, $`3\sigma `$ constraints in the parameter space, as described in the previous section. One can notice that the LRGs impose the best constraint on the $`\mathrm{\Omega }_\mathrm{M}`$$`\mathrm{\Omega }_\mathrm{\Lambda }`$ plane. The main galaxies can only probe a shallow redshift region, so that the cosmological constant does not strongly contribute to the redshift-space clustering, and therefore is less constrained. Although the QSOs can probe deep redshifts, they are too sparse to detect the redshift-space clustering itself. The LRGs nicely have intermediate properties and are balanced between the depth and the density of the survey objects.

In Figure 2 we plot the expected errors when three parameters, $`\mathrm{\Omega }_\mathrm{M}`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`\mathrm{\Omega }_\mathrm{B}/\mathrm{\Omega }_\mathrm{M}`$ are simultaneously evaluated.

In each panel, marginalized errors are plot for each set of two parameters, as described in the previous section. In this case, the QSOs do not constrain each parameter well, because of their sparseness. Main galaxies and LRGs can impose similar constraints for the baryon fraction and the mass density parameter.

In Figure 3 we plot the expected errors when four parameters, $`\mathrm{\Omega }_\mathrm{M}`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`\mathrm{\Omega }_\mathrm{B}/\mathrm{\Omega }_\mathrm{M}`$, $`b`$ are simultaneously evaluated.

The results of the QSOs are not plotted just because the concentration ellipses are too large and comparable to the bounding boxes. In this Figure, the bias parameter $`b`$ is normalized by the fiducial values $`b_{\mathrm{fid}}=1,2,2`$ for each sample, just for convenience. The bias parameters in this Figure are assumed to be constants without redshift-evolution. If this is not the case, one can divide the sample into redshift bins, and the expected errors of bias parameter will be worse, according to the decrease of the sample volume. However, the error estimates of the other parameters in this Figure are approximately independent of the redshift evolution of the bias due to the marginalization over the bias parameter. This is because the information from cross terms in the redshift bins of $`0.1`$ are already ignored in our analysis, and the likelihood function with a fixed bias parameter is not too sensitive to the bias parameter Matsubara & Szalay (2001).

In Figure 4 we plot the expected errors, while simultaneously evaluating five parameters, $`\mathrm{\Omega }_\mathrm{M}`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`\mathrm{\Omega }_\mathrm{B}/\mathrm{\Omega }_\mathrm{M}`$, $`\sigma _8`$, and $`b`$.

The bias parameter $`b`$ and the normalization $`\sigma _8`$ are strongly correlated because their main contribution to the correlation matrix is through the amplitude of the correlations which is proportional to $`b^2\sigma _{8}^{}{}_{}{}^{2}`$. In Figure 5 we plot the expected errors, while simultaneously evaluating another set of five parameters, $`\mathrm{\Omega }_\mathrm{M}`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, $`\mathrm{\Omega }_\mathrm{B}/\mathrm{\Omega }_\mathrm{M}`$, $`n`$, and $`b`$.

Any similar plot can be generated for an arbitrary set of simultaneously determined parameters. There are $`\left(\genfrac{}{}{0pt}{}{7}{1}\right)+\left(\genfrac{}{}{0pt}{}{7}{2}\right)+\mathrm{}+\left(\genfrac{}{}{0pt}{}{7}{7}\right)=127`$ ways to choose which parameters are fixed, so that we need 122 figures more to be complete. Instead we have shown only several representative cases above. In Table 1 are given all the elements of our estimated seven-dimensional Fisher matrices from which all the information can be extracted according to the way described in the previous section.

The listed matrix elements are normalized by fiducial values as $`F_{\alpha \beta }/(\theta _\alpha \theta _\beta )_{\mathrm{fid}}`$. Thus, the inverse of square-root of each diagonal element gives the fractional error relative to the fiducial value of each parameter when a single parameter is determined knowing the other parameters. The square-root of the diagonal elements of the inverse of any sub-matrix gives the fractional error of the parameters when corresponding parameters are simultaneously determined.

Basically, the LRG’s can impose constraints on each parameter comparable to, or even better than main galaxies. The reason that the main galaxies do not work so well (even though their number is 10 times larger than that of LRG’s) is that we only use the linear regime of the clustering. Therefore, the depth of the LRG sample is more critical for the parameter estimation than the higher density of the main galaxies. The quasars are hopeless when many parameters are simultaneously determined, because of their extreme sparseness. In this way, our analysis gives the idea how the dense sampling affects cosmological parameter estimations.

## 5 Conclusions

In this paper, we described a maximum-likelihood method to determine the cosmological parameters from apparent redshift-space clustering. We took advantage of the analytical formula for linear two-point correlations in redshift space with cosmological distortions. We introduced an accurate approximation in calculating the elements of the correlation matrix among cells in redshift-space. This step makes the method practically feasible in spite of the large dimensionality of the correlation matrices.

In the maximum-likelihood method, one can easily evaluate the expected parameter estimation errors in any sample from the Fisher information matrix. We have calculated the seven-dimensional Fisher matrix for three types of objects in the SDSS survey: main galaxies, LRGs and quasars. To illustrate the behavior of the multi-dimensional Fisher matrix we have used concentration ellipses in marginalized two-dimensional parameter space.

In this paper, we divided the survey volume into generic boxes in order to simplify the Fisher matrix estimation. We ignored the correlations between these sub-regions, so the constraints will improve somewhat if those correlations are properly included. However, the inversion of the resulting huge matrices can become extremely time-consuming. The use of the Karhunen-Loève (KL) transform is a practical strategy in this case. Such methods can also be used in a targeted data-compression role to find linear combinations of counts which retain as much information about the parameters as possible (Vogeley & Szalay, 1996; Tegmark, Taylor & Heavens, 1997; Matsubara, Szalay & Landy, 2000; Taylor et al., 2001).

The choice of the cell radius is somewhat arbitrary in this work. We choose the spherical cell with radius of $`10h^1\mathrm{Mpc}`$ for galaxies and LRGs, which is the border of the linear regime. With larger cell radius, the validity of the linear theory increases and the shot noise is reduced. The cosmic variance, however, increases with cell radius. The parameter estimation is dominated by the highest signal-to-noise modes, which are at large wavelength, in particular for the case of the LRG sample. The high frequency modes close to pixel scales mostly contain shot noise after the KL transformation. As a result, we believe our conclusions are not sensitive to the choice of the cell radius. A fully accurate determination of the optimal choice of the cell radius depends on the behavior of the nonlinear effects, so that a comparison with numerical simulation is needed, beyond the scope of the current work.

We have considered three subsets of the SDSS redshift data, spanning a wide range of depth, sampling density and intrinsic clustering strength. We found, that for measuring cosmological parameters in the linear regime there is a clear optimum, represented by the intermediate-redshift LRG’s. The low spatial density of quasars is not overcome by their much larger depth, and the relatively small depth of the main SDSS galaxies is not compensated by their high sampling density – the redshift is not high enough to test curvature, and their cosmic variance is too large. The LRG sample, much smaller in numbers than the main sample, and much shallower than the quasars is an excellent compromise between sampling density and cosmological depth. The constraints derived from the LRGs are much tighter than for the other two samples.

The advantage of these intermediate-redshift objects, and the logic behind this optimum goes beyond the SDSS. In designing future redshift surveys, it is important to find the right balance between the density of objects and the survey depth. Their interplay can be quite complex, as we have shown here. The relation between accuracy and sky coverage is simple and can be estimated analytically.

We would like to acknowledge useful discussions with Dan VandenBerk, Daniel Eisenstein and Adrian Pope. TM acknowledges support from the Ministry of Education, Culture, Sports, Science, and Technology, Grant-in-Aid for Encouragement of Young Scientists, 13740150, 2001. AS acknowledges support from grants NSF AST-9802 980 and NASA LTSA NAG-53503.


text:
# Segmented Band Mechanism for Itinerant Ferromagnetism

## Abstract

We introduce a novel mechanism for itinerant ferromagnetism, which is based on a simple two-band model, and using numerical and analytical methods, we show that the Periodic Anderson Model (PAM) contains this mechanism. We propose that the mechanism, which does not assume an intra-atomic Hund’s coupling, is present in both the iron group and some $`f`$ electron compounds.

Even though itinerant ferromagnetism was the first collective electronic phenomena studied quantum mechanically, the microscopic mechanisms driving this phase are still unknown . In 1928, Heisenberg  formulated his spin model to address this issue, but as Bloch  pointed out, a model of localized spins cannot explain the metallic ferromagnetism observed in Fe, Co and Ni. After seven decades of intense effort we still do not know what is the minimal model of itinerant ferromagnetism and, more importantly, the basic mechanism of ordering.

In 1963, Hubbard  and others introduced the Hubbard model to explain the ferromagnetic (FM) properties of the iron group, incorporating the kinetic energy in a single nondegenerate band with an intra-atomic Coulomb repulsion $`U`$. With the exception of Nagaoka’s  and Lieb’s  theorems, subsequent theoretical approches were not controlled enough to determine whether the Hubbard model has a FM phase. The central issue is the precise evaluation of the energy for the paramagnetic (PM) phase. Because it does not properly incorporate the correlations, mean field theory overestimates this energy and predicts a large FM region . In contrast, numerical calculations have narrowed the extent of this phase to a small region around the Nagaoka point .

Going beyond the simple one-band Hubbard model is advocated, for example, by Vollhardt et. al. . They note that the inclusion of additional density-density interactions, correlated hoppings, and direct exchange terms favors FM ordering. A very simple analysis shows that increasing the density of states $`𝒟(E)`$ below the Fermi energy $`E_F`$ and placing $`E_F`$ close to the lower band edge increases the FM tendency. One can achieve this by including additional hopping terms. The effectiveness of a next nearestneighbor hopping $`t^{}`$ was studied numerically by Hlubina et al  for the Hubbard model on a square lattice. They found a FM state when the van Hove singularity in $`𝒟(E)`$ occurred at $`E_F`$. However, this phase was not robust against very small changes in $`t^{}`$.

Years ago, Slater  and van Vleck  speculated that band degenerancy is an essential precondition for itinerant ferromagnetism. They suggested that the intra-atomic Hund coupling in open shells could be transmitted from one atom to another by the conduction electrons. However, there are FM metals like Ni where the influence of the Hund’s coupling is not clear. In Ni, Hund’s coupling is associated with the $`3d^8`$ configuration which has low probability as the main configurations are roughly $`40\%`$ of $`3d^{10}`$ and $`60\%`$ of $`3d^9`$ . The relevant question for some transition metals is thus whether a model involving just the two configurations is sufficient, or are other orbitals and Hund’s exchange necessary to explain ferromagnetism. Furthermore, there are $`f`$ electron itinerant ferromagnets, like CeRh3B2 , whose only local magnetic coupling is Kondo like, i.e. antiferromagnetic.

The novel mechanism we now introduce emerges from a two-band model, such as the periodic Anderson model (PAM). The basic ingredients are an uncorrelated dispersive band hybridized with a correlated and narrow band. Missing is an explicit intra-atomic Hund’s exchange. We show the PAM supports our mechanism by interpreting the results of quantum Monte Carlo (QMC) simulations with an effective model derived from it.

We will discuss our mechanism in the context of:

$`H`$ $`=`$ $`t_b{\displaystyle \underset{𝐫,𝐫^{},\sigma }{}}(b_{𝐫\sigma }^{}b_{𝐫^{}\sigma }+b_{𝐫^{}\sigma }^{}b_{𝐫\sigma })+V{\displaystyle \underset{𝐫,\sigma }{}}(b_{𝐫\sigma }^{}a_{𝐫\sigma }+a_{𝐫\sigma }^{}b_{𝐫\sigma })`$ (1)
$``$ $`t_a{\displaystyle \underset{𝐫,𝐫^{},\sigma }{}}(a_{𝐫\sigma }^{}a_{𝐫^{}\sigma }+a_{𝐫^{}\sigma }^{}a_{𝐫\sigma })+ϵ_a{\displaystyle \underset{𝐫,\sigma }{}}n_{𝐫\sigma }^a+{\displaystyle \frac{U}{2}}{\displaystyle \underset{𝐫,\sigma }{}}n_{𝐫\sigma }^an_{𝐫\overline{\sigma }}^a,`$ (2)
where $`b_{𝐫\sigma }^{}`$ and $`a_{𝐫\sigma }^{}`$ create an electron with spin $`\sigma `$ in $`b`$ and $`a`$ orbitals at lattice site $`𝐫`$ and $`n_{𝐫\sigma }^a=a_{𝐫\sigma }^{}a_{𝐫\sigma }`$. The $`t_b`$ and $`t_a`$ hoppings are only to nearest-neighbor sites. When $`t_a=0`$, the Hamiltonian is the standard PAM. For the $`f`$ electron compounds, the $`a`$ and $`b`$ orbitals play the role of the $`f`$ and $`d`$ orbitals, and $`t_a0`$. For transiton metals, they correspond to the $`3d`$ and $`4s`$ orbitals. For $`U=0`$, the resulting Hamiltonian $`H_0`$ is easily diagonalized:

$`H_0={\displaystyle \underset{𝐤,\sigma }{}}\left(E_𝐤^+\beta _{𝐤\sigma }^{}\beta _{𝐤\sigma }^{}+E_𝐤^{}\alpha _{𝐤\sigma }^{}\alpha _{𝐤\sigma }^{}\right)`$ (3)
$`E_𝐤^\pm ={\displaystyle \frac{1}{2}}\left[e_𝐤^b+e_𝐤^a\pm \sqrt{(e_𝐤^be_𝐤^a)^2+4V^2}\right]`$ (4)
with $`e_𝐤^b=2t_b_{i=1}^D\mathrm{cos}k_{x_i}`$ and $`e_𝐤^a=ϵ_a2t_a_{i=1}^D\mathrm{cos}k_{x_i}`$ for a hypercubic lattice in dimension $`D`$. The operators which create quasi-particles in the lower and upper bands are:

$`\alpha _{𝐤\sigma }^{}=u_𝐤a_{𝐤\sigma }^{}+v_𝐤b_{𝐤\sigma }^{},\beta _{𝐤\sigma }^{}=v_𝐤a_{𝐤\sigma }^{}+u_𝐤b_{𝐤\sigma }^{}`$ (5)
$`u_𝐤={\displaystyle \frac{E_𝐤^+e_𝐤^a}{\sqrt{(E_𝐤^+e_𝐤^a)^2+V^2}}},v_𝐤={\displaystyle \frac{V}{\sqrt{(E_𝐤^+e_𝐤^a)^2+V^2}}}.`$ (6)
In Fig. 1, we illustrate the (one-dimensional) non-interacting bands for the case of interest: $`ϵ_a`$ close to $`E_F`$ and above the bottom of the $`b`$ band. If $`|V||t_b|`$, we can identify two subspaces in each band where the states have either predominantly $`b`$ ($`\psi `$ subspace) or $`a`$ ($`\varphi `$ subspace) character. The size of the crossover region around the points where the original unhybridized $`b`$ and $`a`$ bands crossed is proportional to $`|V/t_b|`$; that is, it is very small. The creation operators for the Wannier orbitals $`\psi _{𝐫\sigma }^{}`$ and $`\varphi _{𝐫\sigma }^{}`$ associated with each subspace are:

$`\psi _{𝐫\sigma }^{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N}}}\left[{\displaystyle \underset{𝐤𝐊^>}{}}e^{i𝐤𝐫}\beta _{𝐤\sigma }^{}+{\displaystyle \underset{𝐤𝐊^<}{}}e^{i𝐤𝐫}\alpha _{𝐤\sigma }^{}\right]`$ (7)
$`\varphi _{𝐫\sigma }^{}`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{N}}}\left[{\displaystyle \underset{𝐤𝐊^>}{}}e^{i𝐤𝐫}\alpha _{𝐤\sigma }^{}+{\displaystyle \underset{𝐤𝐊^<}{}}e^{i𝐤𝐫}\beta _{𝐤\sigma }^{}\right].`$ (8)
where $`N`$ is the number of sites. The subsets $`𝐊^>`$ and $`𝐊^<`$ are defined by: $`𝐊^>=\{𝐤:|u_𝐤||v_𝐤|\}`$ and $`𝐊^<=\{𝐤:|v_𝐤|>|u_𝐤\}`$. In this new basis:

$`H_0=H_0^\varphi +H_0^\psi ={\displaystyle \underset{𝐫,𝐫^{},\sigma }{}}\tau _{𝐫𝐫^{}}^\varphi \varphi _{𝐫\sigma }^{}\varphi _{𝐫^{}\sigma }^{}+{\displaystyle \underset{𝐫,𝐫^{},\sigma }{}}\tau _{𝐫𝐫^{}}^\psi \psi _{𝐫\sigma }^{}\psi _{𝐫^{}\sigma }^{},`$ (9)
with $`\tau _𝐫^\varphi =\frac{1}{N}[_{𝐤𝐊^>}e^{i𝐤𝐫}E_𝐤^{}+_{𝐤𝐊^<}e^{i𝐤𝐫}E_𝐤^+]`$ and $`\tau _𝐫^\psi =\frac{1}{N}[_{𝐤𝐊^>}e^{i𝐤𝐫}E_𝐤^++_{𝐤𝐊^<}e^{i𝐤𝐫}E_𝐤^{}]`$.

Because the $`U`$ term in $`H`$ involves only the $`a`$ orbitals, the matrix elements of $`H`$ connecting the $`\varphi `$ and $`\psi `$ subspaces are small compared to the characteristic energy scales of the problem (the matrix elements of $`H`$ within the subspaces). To see this we express $`a_{𝐫\sigma }^{}`$ as a function of $`\varphi _{𝐫\sigma }^{}`$ and $`\psi _{𝐫\sigma }^{}`$ by first inverting Eqs. 6 and 8 to find:

$$a_{𝐫\sigma }^{}=\underset{𝐫^{}}{}W_{𝐫𝐫^{}}\varphi _{𝐫^{}\sigma }^{}+w_{𝐫𝐫^{}}\psi _{𝐫^{}\sigma }^{}$$
(10)
with $`W_𝐫=\frac{1}{N}[_{𝐤𝐊^>}e^{i𝐤𝐫}u_𝐤_{𝐤𝐊^<}e^{i𝐤𝐫}v_𝐤]`$ and $`w_𝐫=\frac{1}{N}[_{𝐤𝐊^>}e^{i𝐤𝐫}v_𝐤+_{𝐤𝐊^<}e^{i𝐤𝐫}u_𝐤]`$. However, because the $`\varphi `$ orbitals have predominantly $`a`$ character, while the $`\psi `$ orbitals have predominantly $`b`$ character, $`a_{𝐫\sigma }^{}_𝐫^{}W_{𝐫𝐫^{}}\varphi _{𝐫\sigma }^{}`$. If $`|V||t_b|`$, then $`|w_𝐫||W_𝐫|`$. Consequently, the $`a`$ subspace becomes invariant under the application of $`H`$. In addition, because $`|W_\mathrm{𝟎}^{}||W_{𝐫\mathrm{𝟎}}^{}|`$, we can establish a hierarchy of terms where the lowest order one corresponds to a simple on-site repulsion:

$$H_{eff}^U=\stackrel{~}{U}\underset{𝐫}{}n_𝐫^\varphi n_𝐫^\varphi $$
(11)
with $`\stackrel{~}{U}=U|W_\mathrm{𝟎}|^4`$ and $`n_{𝐫\sigma }^\varphi =\varphi _{𝐫,\sigma }^{}\varphi _{𝐫,\sigma }^{}`$. The next order terms, containing three and two $`W_\mathrm{𝟎}`$ factors, are much smaller and are essentially the same as the intersite interactions which in the past were added to the Hubbard model to enhance the ferromagnetism . Adding $`H_{eff}^U`$ to $`H_0`$ we get the effective Hamiltonian:

$`H_{eff}={\displaystyle \underset{𝐫,𝐫^{},\sigma }{}}(\tau _{𝐫𝐫^{}}^\psi \psi _{𝐫\sigma }^{}\psi _{𝐫^{}\sigma }^{}+\tau _{𝐫𝐫^{}}^\varphi \varphi _{𝐫\sigma }^{}\varphi _{𝐫^{}\sigma }^{})+\stackrel{~}{U}{\displaystyle \underset{𝐫}{}}n_𝐫^\varphi n_𝐫^\varphi `$ (12)
The $`\psi `$ and $`\varphi `$ orbitals form uncorrelated and correlated non-hybridized bands: $`H_{eff}=H^\psi +H^\varphi `$. For the $`\varphi `$ orbitals we obtain an effective one band Hubbard model with the peculiar double shell like dispersion relation shown by the thick lines in Fig. 1.

Particularly for $`t_a=0`$, $`H^\varphi `$ has a very large density of states in the lower shell of the $`\varphi `$ band  which is located near $`ϵ_a`$. From Fig. 1 it is also clear that the electrons first doubly occupy the uncorrelated $`\psi `$ band states which are below $`ϵ_a`$. However, when $`E_F`$ gets close to $`ϵ_a`$, i.e. the system is in the mixed valence regime, the electrons close to the Fermi level go into some of the correlated $`\varphi `$ states. Then, the interaction term $`H_{eff}^U`$, combined with the double shell band structure of $`H_0^\varphi `$, gives rise to a FM ground state (GS): The electrons close to $`E_F`$ spread to higher unoccupied $`𝐤`$ states and polarize, which causes the spatial part of their wave function to become antisymmetric, eliminating double occupancy in real space and reducing the Coulomb repulsion to zero. The cost of polarizing is just an increase in the kinetic energy proportional to $`\delta _a\mathrm{}v_F\delta _k`$, where $`v_F`$ is the Fermi velocity and $`\delta _k`$ is the interval in $`𝐤`$ space in which the electrons are polarized.

To determine the stability of this unsaturated FM state, we compare its energy with that of the PM (nonmagnetic) state. If we were to build a nonmagnetic state with only the states of the lower $`\varphi `$ shell, we would find a restricted delocalization for each electron because of the exclusion of the finite set of band states ($`𝐤`$-states) in the upper shell. To avoid the Coulomb repulsion $`U`$ for double occupying a given site, the electrons need to occupy all $`𝐤`$-states. This means they have to occupy the $`\varphi `$ states in the upper and lower shells. This restricted delocalization is a direct consequence of Heisenberg’s uncertainty principle, and the resulting localization length depends on the wave vectors, where the original $`b`$ and $`a`$ bands crossed, that define the size ($`\mathrm{\Delta }𝐤`$) of each shell. The energy cost for occupying the $`\varphi `$ states in the upper shell is proportional to the hybridization gap $`\mathrm{\Delta }`$. Therefore if $`U`$ is the dominant energy scale in the problem and $`\mathrm{\Delta }\delta _a`$, the FM state lies lower in energy than the nonmagnetic state. Under these conditions, the effective FM interaction is proportional to the hybridization gap $`\mathrm{\Delta }`$.

This mechanism for ferromagnetism on a lattice is analogous to intra-atomic Hund’s mechanism polarizing of electrons in atoms. In atoms, we also have different degenerate (the equivalent of $`\delta _a`$ is zero) shells separated by an energy gap. If the valence shell is open, the electrons polarize to avoid the short range part of the Coulomb repulsion (again reflecting the Pauli exclusion principle). The energy of an eventual nonmagnetic state is proportional either to the magnitude of the Coulomb repulsion or to the energy gap between different shells. The interplay between both energies set the scale of Hund’s intra-atomic exchange coupling.

The FM mechanism just described applies to any finite dimension. For a chain of 16 sites we calculated, by the Lanczos method, the exact GS of $`H_{eff}`$ for $`ϵ_a=t_b`$, infinite $`\stackrel{~}{U}`$ ($`\stackrel{~}{U}|\tau _𝐫^\varphi |`$), and different values of $`V`$ as a function of electron concentration $`n=N_e/4N`$. We found that the GS is a nonsaturated ferromagnet between $`n=1/4`$ (one electron per site in the PAM) and $`n=3/8`$. In the local momentum regime QMC  and DMRG  calculations report a FM phase in a very similar range of $`n`$, in contrast to the much broader range found by dynamical mean field theory calculations . The largest magnetization $`M`$ is obtained when $`n`$ is such that the lower shell of the $`\varphi `$ band is completely filled and the upper one is nearly empty; i.e., when $`E_Fϵ_a`$.

In Fig. 2a we plot the PAM’s energy per site $`E/N`$, computed by our QMC method , as a function of the total spin per site $`S/N`$ for chains of varying length $`N`$ and fixed electron density $`n`$ . Over these chain lengths the data collapse, with $`E(S)/N`$ showing the minimum $`E_{GS}/N`$ at a non-zero value of $`S/N`$ that represents a very good estimate of the magnetization $`M/N`$ for the PAM in the thermodynamic limit (TL). In Figs. 2b and 2c, we show $`\mathrm{\Delta }E/N(E_{GS}E(S=0))/N`$ and $`M/N`$ for one and two dimensional systems as a function of $`1/N`$. In both cases $`M/N`$ smoothly varies to non-zero values in the TL. The non-smooth variation of $`\mathrm{\Delta }E/N`$ in two dimensions is a consequence of shell effects still present in a finite sized system, but is clearly suggestive of its likely extrapolation to a non-negative value for very large $`N`$. In Fig. 2d we plot our QMC results for $`(E(FM)E(0))/E(0)`$ as a function of $`ϵ_f`$ for a $`8\times 8`$ cluster and two different electron densities. As we increase $`ϵ_a`$ (starting from below the bottom of the $`b`$-band) $`\mathrm{\Delta }E`$ decreases and then increases. The most stable FM state occurs when $`E_Fϵ_a`$ as expected from our discussion of the effective model. Increasing $`ϵ_a`$ even further, $`\mathrm{\Delta }E`$ approaches zero indicating that the spin of the GS decreases (see Fig. 2e) and the system becomes a paramagnet.

By changing $`ϵ_a`$ we can evolve the system from the localized $`(ϵ_aE_F)`$ to the itinerant case $`(ϵ_aE_F)`$. From Figs. 2d and 2e we see that the itinerant and the localized FM phases are apparently continuously connected. With decreasing $`ϵ_a`$, $`|\mathrm{\Delta }E|`$ decreases while the zero temperature magnetization $`M`$ increases. The strong reduction of $`|\mathrm{\Delta }E|`$ is a result of the very small effective magnetic interaction in the localized limit ($`J_{RKKY}`$ is order $`V^4`$ ). Decreasing $`ϵ_a`$ increases the population of the lower $`\varphi `$-shell. Since most of the electrons in the lower $`\varphi `$ shell are polarized when $`ϵ_aE_F`$, this leads to an increase of $`M`$. These results are consistent with the observed behavior of $`M`$ and the Curie temperature $`T_c`$ in LaxCe1-xRh3B2 as function of $`x`$ .

Even though we cannot do finite temperature calculations with our present version of the QMC method, we can discuss, at least qualitatively, the predictions of our mechanism for finite temperatures. If we move $`ϵ_a`$ above $`E_F`$, the number of $`a`$-electrons decreases together with the magnetization $`M`$. A new energy scale $`ϵ_aE_F`$ emerges. We propose that this scale is responsible for the finite temperature peak in the magnetization of Ce(Rh1-xRux)3B2  that suggests an ordered state with high entropy. At $`T=0`$, $`M`$ is small because of the reduced number of $`a`$-electrons. When the temperature is of the order of $`ϵ_aE_F`$, electrons are promoted from the doubly occupied $`b`$ states to the unoccupied $`a`$ states which have a large entropy (large density of states). These $`a`$-electrons polarize because of the energy considerations discussed above. The source of the large entropy is thus associated with charge and not with spin degrees of freedom, which explains why a state with larger $`M`$ has a higher entropy. From this analysis we predict that the entropy below $`T_c`$ contains a considerable contribution from the the charge degrees of freedom.

We can also connect our mechanism with the hydrostatic pressure dependence of $`T_c`$. To do this we calculated $`|\mathrm{\Delta }E|/N`$ by the QMC method as function of increasing $`t_b`$ (Fig. 3a). Here we are assuming that the main effect of the hydrostatic pressure is to increase $`t_b`$ and to leave the other parameters unchanged. The order of magnitude of $`|\mathrm{\Delta }E|/N`$, which should be proportional to $`T_c`$, and its qualitative behavior in Fig. 3a are in good agreement with the experimental results in CeRh3B2 . We see from Fig. 3a that for the itinerant FM case, $`|\mathrm{\Delta }E|`$/N is of the order of 100K. This scale is much larger than the magnitude of the RKKY interaction  ($`1^{}`$K) which is commonly used to explain the origin of the magnetic phase when the $`a`$ electrons are localized. We also find that the FM state appears close to quarter filling and disappears for $`n`$ close to $`3/8`$.

We used a nonzero $`t_a`$ to study the stability of the FM phase when $`\delta _a`$ is varied. This study is also important because, in contrast to $`f`$ electron compounds, in the FM transition metals both bands are dispersive. In Fig. 3b we see that the FM phase is even more stable for $`t_a0.1t_b`$ than for $`t_a=0`$ and becomes unstable for $`t_a0.05t_b`$. The reason for this asymmetric behavior is easy to understand in terms of the variation of $`\delta _a`$: If $`t_a`$ is negative, then the effect of $`t_a`$ on the dispersion of the $`\varphi `$ band is opposite to that of the hybridization $`V`$. When $`t_a0.1t_b`$ we get, for the given $`ϵ_a`$ and $`V`$, the minimum value for $`\delta _a`$ and therefore the most stable FM case. When we depart from this value of $`t_a`$, $`\delta _a`$ increases, $`|\mathrm{\Delta }E|`$ decreases, and the FM state becomes less stable.

In summary, we introduced a novel mechanism for itinerant ferromagnetism which is present in a simple two band model. The picture just presented, combined with our previous results , allows a reconciliation of the localized and delocalized ferromagnetism pictures painted by Heisenberg  and Bloch . The hybridization between bands plays a crucial role. We have also considered the case relevant for the iron group where the dispersion of the lower band is not negligible. The fact that the ferromagnetism is even more stable for finite values of $`t_a`$ indicates that our mechanism is relevant to explain the ferromagnetism of the transition metals, like Ni, where a correlated and narrow $`3d`$ band is hybridized with the $`4s`$ band. It suggests that the ferromagnetism in the transition metals can originate, at least in part, in the interplay between the correlations and the particular band structure, and not solely in the intra-atomic Hund’s exchange . In addition, our results explain several qualitative features observed in the ferromagnet CeRh3B2 . Elsewhere we will discuss the relation of our mechanism to RE(Co1-xSix)2 (RE=Ho, Er) and the uranium monochalcogenides (US, USe, and UTe) which are ferromagnets .

Acknowledgements. This work was sponsored by the US DOE. We acknowledge useful discussions with A. J. Arko, B. Brandow, J. J. Joyce, J. M. Lawrence, S. Trugman, G. Ortiz, and J. L. Smith. We thank J. M. Lawrence for pointing out the experimental work on the Ce compounds. J. B. acknowledges the support Slovene Ministry of Education Science and Sports and FERLIN.