# Coexistence of antiferromagnetism and superconductivity within - model with strong correlations and nonzero spin polarization

###### Abstract

The coexistence of antiferromagnetism with superconductivity is studied theoretically within the - model with the Zeeman term included. The strong electron correlations are accounted for by means of the extended Gutzwiller projection method within a statistically-consistent approach proposed recently. The phase diagram on the band filling - magnetic field plane is shown, and subsequently the system properties are analyzed for the fixed band filling . In this regime, the results reflect principal qualitative features observed recently in selected heavy fermion systems. Namely, (i) with the increasing magnetic field the system evolves from coexisting antiferromagnetic-superconducting phase, through antiferromagnetic phase, towards polarized paramagnetic state, and (ii) the onset of superconducting order suppresses partly the staggered moment. The superconducting gap has both the spin-singlet and the staggered-triplet components, a direct consequence of a coexistence of the superconducting state with antiferromagnetism.

###### pacs:

71.27.+a, 74.20.-z, 74.70.Tx, 74.25.Ha^{†}

^{†}preprint: APS/123-QED

## I Introduction

The interplay of antiferromegnetism (AF) with superconductivity (SC) is one of the important topics in condensed-matter physics, Rev_Demler () as better understanding of this subject would improve our knowledge of a number of systems such as high-Tc, Rev_High_Tc () heavy-fermion, Pfleiderer () and organic Rev_Organic_SC () superconductors. In all those systems, superconductivity appears in the vicinity of magnetic phases (mostly antiferromagnetic, but also ferromagnetic Saxena (); Aoki ()). Moreover, magnetic interactions or fluctuations are very frequently considered to be the pairing mechanism in unconventional superconductors. Monthoux (); MonthouxSF (); MonthouxSF2 () Typically, antiferromagnetism and superconductivity are competing quantum phenomena because of the competition between the Meissner-supercurrent screening and the internal-field generation by magnetic ordering. This antagonism can be overcome by a spatial separation of the AF and the SC phases or by subdivision of the electrons into more localized (resulting in AF) and more itinerant (participating in SC) parts. However, especially interesting is the situation, when the same electrons are involved in both phenomena, as is the case for some heavy-fermion systems. There, SC and AF can coexist easily, when the periodicity of magnetic structure is much smaller than the coherence length for the Cooper pair. In other words, when , the staggered exchange field averages out to zero within the coherence volume. In this respect the Ce-based ’115’ heavy-fermion compounds - the family of CeMIn (with ) Sarrao (); Thompson (); Shishido () is the most promising, as both antiferromagnetism and superconductivity are believed to arise from electrons, where even the interplay of the two orders can be studied by tuning the system with pressure, magnetic field, or doping.

Also, recently a competitive coexistence of AF and SC has been reported in both CeRhIn Chen (); Knebel (); Knebel2 (); ParkPNAS () and CeCo(InCd). Nair (); Urbano () In the latter system, a mutual influence of AF and SC has been observed. Namely, it turns out that the onset of SC order with lowering temperature prevents any further increase of the antiferromagnetic magnetization. Nair () A similar type of coexistence has also been observed in CeRhSi. Kimura ()

Generally, in the heavy-fermion systems strong correlations among electrons are the reason for an emergence of new and nontrivial physics. Those nontrivial features should be properly accounted for when modeling those systems. In the present paper, an investigation of the coexistence of AF with SC in an applied magnetic field is presented. To account for strong electron correlations, the Gutzwiller-projected - model is used with the Zeeman term included. The extended Gutzwiller scheme proposed recently Fukushima () is utilized for calculation of statistical averages of the relevant operators. Our model, although at first sight seems too simplified to be related to heavy fermion systems, it nonetheless reflects qualitatively principal features observed recently in selected heavy fermion systems.

It is commonly believed that the minimal model for investigation of heavy-fermion systems should be the two-band Periodic Anderson Model (PAM) (see e.g. Ref. Sacramento, ) or the Kondo lattice model. Tsunetsugu () On the other hand, the one-band calculations have already proved fruitful in the analysis of AF and SC coexistence in CeRhIn, Alvarez () as well as in investigations of the high-field low-temperature unconventional superconducting phase of CeCoIn. Yanase (); Aperis () The narrow-band limit of PAM has been discussed theoretically also elsewhere (see Refs. JS_RSP1, and Maska, , Appendix A). Generally, it appears when only a single hybridized band is involved and in the heavy-fermion limit (i.e. when -level occupancy , with ). Simply put, the --type model reflects the physics of those hybridized and strongly correlated systems in the narrow -band limit.

## Ii Model, order parameters, and constraints

We start from the - model with the Zeeman term included, as represented by the Hamiltonian^{1}^{1}1For the discussion of various forms of - model see Ref. JJ2, below.

(1) |

where denotes the summation over bonds, and is the spin -component. Since its derivation, ChaoJS (); ZR () the - model represents an active field of research (see e.g. Ref. Moreno, for a recent analysis of the one-dimensional situation). The - model captures the essential ingredients of physics of the high-Tc superconductors. The advantage of using this model is that both AF and SC come from a microscopic parameter - antiferromagnetic exchange and therefore there are no phenomenological terms in the Hamiltonian (as opposed to some earlier studies of AF and SC coexistence). We neglect the orbital effects, as the Maki parameter Maki () in the systems of our interest here is high. Kumagai2 (); Knebel2 () The Gutzwiller projector eliminates double occupancies in real space. In the following we will use the more general correlator

(2) |

where are the so-called fugacity factors. Also, this correlator connects the correlated and uncorrelated wave functions, Anderson () via

(3) |

This allows to express average of any operator in the correlated state as

(4) |

where . With the above equation one can in principle calculate average value of Hamiltonian (1), namely

(5) |

but this is a nontrivial task, as after applying the Wick theorem too many terms appear (see Ref. Fukushima, , e.g. Eq. (8) and the discussion afterwards), and one has to resort to making approximations at this point. There are a few ways to perform this operation, and this is still an active field of research, so one can expect new calculation schemes to appear Bunemann (). Here, we use the scheme proposed recently by Fukushima Fukushima (); Fukushima2 () in the local-constraint version, which assumes that the average number of particles at any site and with any spin is unchanged by the projection,

(6) |

This formalism is known to reproduce the Variational Monte Carlo results better than the conventional Gutzwiller approximation (at least, for the projected uniform nonmagnetic -wave BCS superconductor - see Figs. 3 and 4 of Ref. Fukushima, ). The local-constraint version of the formalism is quite general in the sense that it is capable of accounting for antiferromagnetism, superconductivity, and the ferromagnetic polarization. The explicit expressions for all averages appearing in Eq. (5) are given in Ref. Fukushima, . To express them in terms of mean-fields of our interest, we need to assume what is the character of the uncorrelated wave function . Since our goal is the description of coexistence of AF and SC, we assume the corresponding mean-fields as nonzero at the level of as in the following. We start with the particle number in the form

(7) |

where is the band filling (assumed as constant), is the ferromagnetic (longitudinal) spin-polarization component, and is the antiferromagnetic (staggered) spin polarization. The factor (with ) is responsible for the sign reversal of the staggered magnetic moment when exchanging the two sublattices A and B. Fazekas () We also assume the superconducting order parameter can be decomposed into two components

(8) |

where ensures the -wave gap symmetry by setting for () respectively, and with , being the square-lattice basis vectors. The -wave solution (of the form) is taken throughout in the following analysis, as it is the most favorable energetically (cf. e.g. Ref. Ogata_tJ_Review, ) The superconducting order parameter can be rewritten in terms of the singlet and the staggered -triplet components, namely

(9) | |||||

with

(10) | |||||

(11) |

The superconducting order parameter is defined on bond (nearest-neighbor pair of sites). To define the gap per site, we make use of the standard Kyung () relation for -wave solution

(12) | |||||

(13) |

where denotes the nearest neighbors of site . The existence of the triplet component is inevitable even if there is no triplet channel in the pairing potential. Namely, the triplet component is dynamically induced by the singlet gap and antiferromagnetism. Psaltakis (); Kyung (); Tsonis (); Aperis () From a microscopic point of view, this is also not surprising at all (see Fig. 1 for an intuitive illustration).

An interesting feature of the superconducting gap defined by Eq. (9) is the nonzero momentum of Cooper pairs for the triplet component (it results from the term, in an analogy to center-of-mass momentum in the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) phase FF (); LO (); *LO2; Matsuda-Rev ()). A superconducting state with nonzero momentum has been investigated in a number of cases, Aperis (); FF (); LO (); *LO2; Ogata_Q () even in zero external magnetic field. Loder (); Maska2 () The one presented here is analogous to that from Ref. Aperis, , except we consider both the microscopic - model and the limit of strong correlations.

With the above assumptions, we can express the ground-state energy [cf. Eq. (5)] as a function of the band filling , the magnetization components and , the superconducting gaps and , and the hopping amplitudes . We assume as nonzero first- and second-nearest-neighbor hopping integrals and , which yields 6 different hopping amplitudes

(14) |

The resulting expression for is quite lengthy and has been presented in Appendix A. Next, as in the method proposed earlier JJ2 (); JJ1 (); Acta (); SGA () (our present formulation is analogous to that from Ref. JJ2, ), to solve the model in a statistically-consistent manner, we impose additionally the constraints on all introduced mean fields by means of the Lagrange multipliers method. In effect, to carry out the subsequent analysis, we use the following energy operator

(15) | |||||

This method of approach is equivalent in the () limit to that presented in Refs. Wang, and Yang, . The equivalence can be seen from the comparison of Eq. (15) and Eqs. (30)-(37) with the corresponding equations from Refs. Wang, and Yang, (e.g. Eq. (13) from Ref. Wang, provides an effective Hamiltonian with the operator part equivalent to our ). The Lagrange multipliers , , and have the same symmetries, as the corresponding mean fields , , and . We also assume they are spatially homogeneous. Namely,

(16) | |||||

(17) |

with

(18) | |||||

(19) |

After performing Fourier transformation of the operator part of we obtain

(20) | |||||

where the primed summation runs over the folded (magnetic) Brillouin zone, by we denote all the mean-fields, is the total number of sites, and the four-component operator has the following components

(21) |

The matrix is given as

(22) |

where

(23) | |||||

(24) | |||||

(25) | |||||

(26) | |||||

(27) |

We have also used the fact that . Note that in the present formulation corresponds to sum of the magnetic field and the correlation-induced field JKJS (); JKJS2 () (or equivalently the Lagrange multiplier in the slave-boson theory Korbel (); KR (); SBREVIEW ()). Namely, , what is evident from comparison of Eq. (24) with appropriate expressions taken from Refs. JKJS, ; JKJS2, ; Korbel, .

Next, we determine the eigenvalues of , as they correspond to quasiparticle excitations of the system. An analytic diagonalization of produces very long expressions, and more importantly, expressions with square roots of possibly negative numbers. Therefore, having in mind their subsequent implementation to calculate the physical properties, we diagonalize this matrix numerically. Next, having determined the eigenvalues , we determine the generalized grand potential functional for the system of fermions, which is

(28) | |||||

The physical (equilibrium) values of the mean fields and the Lagrange multipliers are obtained from the necessary conditions for to have a minimum subject to the constraints, i.e.,

(29) |

where by we denote collectively the Lagrange multipliers. Equations provide the explicit analytic expressions for the Lagrange multipliers, i.e.,

(30) | |||||

(31) | |||||

(32) | |||||

(33) | |||||

(34) | |||||

(35) | |||||

(36) | |||||

(37) |

The above expressions can be utilized to eliminate Lagrange multipliers from the solution procedure. Thus, we obtain 11 equations to be solved numerically for the mean fields , instead of 22 equations for both and . The equations for the mean fields (obtained from ) have the following form

(38) | |||||

(39) | |||||

(40) | |||||

(41) | |||||

(42) | |||||

(43) | |||||

(44) | |||||

(45) |

where

(46) |

The derivative is computed numerically with a 5-point stencil method (as it gives two-three orders of magnitude better precision than the standard 3-point stencil). For example

(47) | |||||

where we use the “equilibrium” values of as given by Eqs. (30)-(37). The step is typically equal to . Larger values of would cause greater error in the above formula. Smaller values would result in a loss of numerical precision. We have verified that at (where analytical formulas for the eigenvalues are available) the numerical computation of the derivatives according to Eq. (47) with the chosen step introduces error smaller than the precision of the procedure of solving the set of Eqs. (38)-(45). The value of the step has been chosen after an analysis of the error at and the numerical-precision loss (for very small the numerical-precision loss lead to impossibility of solving the set of Eqs. (38)-(45) with the given precision).

## Iii Results and physical discussion

The equations (38)-(45) are solved numerically with the use of GNU Scientific Library (GSL) GSL () on a grid of size . We use the `gsl_multiroot_fsolver_hybrids`

solver which implements the hybrids algorithm. We use the precision . Namely, the procedure converges when the relation is fulfilled (where the sum is taken over all equations, which have been brought to the form and divided by to ensure lattice-independent convergence conditions). We assume the following values of parameters: , , , and , what yields the temperature . In Table I the exemplary numerical values of the parameters have been provided for the sake of completeness. Numerical accuracy is on the level of the last digit specified. The energy scale has been set by taking the value of the exchange integral as unit, . For more details on the numerical procedure see Chapter 6 of Ref. JKPHD, .

Table I. Equilibrium values of mean-field variables, Lagrange multipliers, | |||
---|---|---|---|

free energy and grand potential functional at and . | |||

Variable | Value | Variable | Value |

3.2997360 | -5.1331996 | ||

0.0000000 | 0.3000001 | ||

0.8100315 | 2.5817963 | ||

0.0922998 | 0.2730140 | ||

0.0479298 | 0.4395630 | ||

0.1218625 | 0.4258402 | ||

0.1218625 | 0.4258402 | ||

-0.0167505 | -0.1031297 | ||

0.0275895 | -0.0120561 | ||

0.0275895 | -0.0120561 | ||

-0.0167505 | -0.1031297 | ||

-1.0110048 | -4.2117488 |

A number of stable phases emerge as solutions of the equations, depending on the physical condition (, ). As we work with constant number of particles , the stable phase is the one with the lowest free energy, defined by

(48) |

where all the optimal values of mean fields and Lagrange multipliers (i.e. those being solution to Eqs. (30)-(45)) are inserted in the functional and is the chemical potential.

The exemplary phase diagram on the band filling - magnetic field plane is exhibited in Fig. 2. It can be seen that antiferromagnetic phase is the predominant one in the low-field regime and above . For antiferromagnetism coexists with superconductivity, what amounts to a phase with nonzero three order parameters (similarly, as in e.g. Ref. Aperis, ). In the low- part of the phase diagram (for ) the saturated ferromagnetic (SFM) phase with becomes the stable state. This phase is stable even in the limit. This is an interesting result, which adds to the discussion of ferromagnetism in the - Koretsune (); *Koretsune2; Putikka () model. There is also a number of papers (see e.g. Refs. JJ2, and Marcin1, ) analyzing the - model (1) with the Gutzwiller-type of approach with the parameters in a similar range (i.e. with and similar values of and ). Some of those papers disregard completely the Zeeman-term influence, and this omission is justified when applying the model to high-Tc superconductors, where orbital effects dominate over the Pauli magnetism. We have shown, that even at the system may be completely spin-polarized and therefore, the inclusion of ferromagnetic polarization is important in treatment of the - model. Finally, our phase diagram can be compared (although this is not a direct comparison, as even at we have for the AF phase) to that obtained by other Gutzwiller approximation scheme, Ogata () in which the coexisting phase was stable up to the doping , and for higher doping levels the pure superconducting state was stable. In our approach the antiferromagnetic phase is stable for such dopings instead. We comment on the strong antiferromagnetism in the following analysis.

At band filling the phase diagram (or the phase sequence as a function of field ) resemble those observed recently in the heavy-fermion compounds CeCo(InCd) Nair () at doping and CeRhSi Kimura () at pressure .^{2}^{2}2Although in CeRhSi the phases have only been analyzed as a function of temperature. Namely, in low magnetic fields a phase with coexisting antiferromagnetic and superconducting orders (AF+SC) is stable, whereas for higher magnetic fields a continuous transition to a pure antiferromagnetic (AF) phase takes place, followed by a discontinuous transition to the polarized paramagnetic (PP) phase. The phases appearing at this band filling () are analyzed thus in detail in the following.

In Fig. 3 we have plotted the free energy curves for a choice of possible a priori phases. It can be seen (also from the following Figures) that the transition AF+SC AF is continuous, whereas the transition AF PP is of the first order. Also, pure superconducting (SC) solution is unstable, and this holds for other band fillings as well. It can be concluded from Fig. 3 that antiferromagnetism is the “dominating” phenomenon, since the energy gain from developing antiferromagnetic order (which can be seen from closer look at the difference ()) is much higher than the gain from developing superconducting order (). Moreover, the energy gain from developing AF order within SC phase () is much higher than that from developing SC order within the AF phase (). This observation can be compared to the results of the Variational Monte Carlo method, Himeda () in which the -wave solution is only slightly higher in energy than the coexisting phase (see Fig. 2 of Ref. Himeda, ).

In Fig. 4 we exhibit the magnetic moment per site of the system for different phases. Namely, we plot the staggered magnetization and the ferromagnetic magnetization (spin-polarization). The staggered magnetization is close to the saturation value of . Such overestimation of the staggered magnetization value over the Variational Monte Carlo results Himeda (); Giamarchi () is also present in the slave-boson approach. Inaba () This is not surprising, as the method we use JJ2 (); JJ1 (); Acta (); SGA () is similar in structure (the Lagrange multipliers are present in both methods) to the slave-bosons approach (for the discussion of the equivalence for the paramagnetic state see Ref. SGA, ). The obtained ferromagnetic spin-polarization for the pure AF phase is equal to at all magnetic fields. Also, it can be seen that development of the SC order within the AF phase alters by a small amount the staggered magnetization , which drops by approximately (see Fig. 4a).

In Fig. 5 various superconducting gaps are shown. Namely we exhibit both the “uncorrelated” gap for the wave function , as well as the gap for the correlated wave function , the latter defined by Eq. (4) and labeled as . In the pure SC phase the sublattice gaps are equal (), what amounts to the absence of the triplet component. Note that although the uncorrelated gaps (, ) are larger in the pure SC phase than in the AF+SC phase, the correlated gaps (, ) are much larger in the AF+SC phase than in the pure SC phase. This very important conclusion means that the presence of antiferromagnetism supports superconductivity in the present situation. The opposite is not true as the staggered moment is slightly larger in the AF phase than in the AF+SC phase. Finally, the renormalized gaps are more than an order of magnitude smaller than their bare (uncorrelated) correspondants.

The picture with large antiferromagnetic magnetization (Fig. 4) and small superconducting gap (Fig. 5) is consistent with the energy curves displayed in Fig. 3. To shift the energy balance towards the SC phase one could either decrease , or increase . By doing that within a wide parameter margin, the antiferromagnetic phase still remains a predominant phase. Other possibility to weaken antiferromagnetism is to include additionally the intersite attraction () in the starting Hamiltonian. This has been shown to stabilize the -wave superconducting state Yanase () (see Ref. Fukushima2, for the expression for the average value of this term within the extended Gutzwiller scheme we use). The strong antiferromagnetism may represent an apparent feature of the Gutzwiller scheme used, Fukushima () in which magnetization is not changed by the projection, as follows directly from Eq. (6).

Finally, in Fig. 6 we display the quasiparticle energies (the Slater subbands) for the phases discussed above for . The crossing of one of the bands with the zero energy line at the S point of the Brillouin zone in Figs. 6bc means that the quasiparticle excitations will be spontaneously created (are gapless), a circumstance leading to a nonzero spin-polarization (cf. Fig. 6), similarly as in the situation for the FFLO state. FF () Note also that for the AF+SC electronic structure is gapful for the -wave superconducting phase, because of the presence of the Slater-type (magnetic) splitting. This is not true anymore for (cf. Figs. 6bc and Fig. 4a) when a uniform ferromagnetic component appears. Also, the bands are sizably wider in the PP state.

## Iv Summary

In summary, we have carried out a detailed analysis of the coexistence of antiferromagnetism and superconductivity within a microscopic - model, with the Zeeman term included. The strong correlations were accounted for by means of the extended Gutzwiller projection method. Also, the constraints assuring the statistical-consistency of the results have been included. We have obtained the phase diagram on the band filling-magnetic field plane, in which for the band fillings in the range and with the increasing magnetic field, a series of phase transitions takes place. Namely, the system evolves from the coexisting (AF+SC) phase, through the antiferromagnetic (AF) phase, to the polarized paramagnetic (PP) phase. Also, the onset of superconducting order reduces the AF order parameter. By contrast, the superconducting gaps are enhanced by the presence of the AF order. In the AF+SC phase there are two superconducting gaps of an almost equal amplitude: the singlet and the staggered-triplet components. These features reflect in a qualitative manner the experimental findings in the CeCo(InCd) Nair (); Urbano () and CeRhSi Kimura () heavy fermion systems, although our model is too simplified to be quantitatively related to such complex heavy fermion systems. Additionally, both antiferromagnetism and superconductivity originate from the same electrons and are driven by the same kinetic-exchange interaction. Note that the real-space pairing is the pairing without “boson glue”, i.e. without paramagnons. It is the mechanism of pairing arising entirely from interelectron correlations which is particularly effective when renormalized hopping and exchange interaction are of comparable magnitude.

As mentioned earlier, it would be very interesting to perform similar analysis within the Periodic Anderson Model, as this might allow for a comparison with the experiments [work along this line is in progress in our group Olga ()]. Also, testing other Gutzwiller schemes seems crucial to verify if the strong antiferromagnetism and the tendency towards saturated ferromagnetism are only the characteristic features of the utilized renormalization scheme, or represent a universal tendency of the projected - model. For that purpose, the inclusion of realistic, orbitally degenerate level structure, not just pseudospin doublet of Ce, would be desirable.

One should also note that the present approach includes the effect of applied magnetic field only via Zeeman term (the Pauli limit). For discussion of high-temperature superconductivity for the orbital effects should be incorporated.

## Acknowledgements

The work was supported by Ministry of Higher Education and Science, Grants Nos. N N202 173735 and N N202 128736, as well as by the Foundation for Polish Science under the ”TEAM” program for the years 2011-2014. Discussions with Olga Howczak and Jakub Jȩdrak are appreciated.

## Appendix A Explicit expression for

We provide here the expression for . This expression can be divided into parts coming from different terms of Hamiltonian with and , as follows . The expressions for and are given by

(49) | |||||

and

(50) | |||||

## References

- (1) E. Demler, W. Hanke, and S.-C. Zhang, Rev. Mod. Phys. 76, 909 (2004)
- (2) P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006)
- (3) C. Pfleiderer, Rev. Mod. Phys. 81, 1551 (2009)
- (4) J. Wosnitza, J. Low Temp. Phys. 146, 641 (2007)
- (5) S. Saxena, P. Agarwal, K. Ahilan, F. Grosche, R. Haselwimmer, M. Steiner, E. Pugh, I. Walker, S. Julian, P. Monthoux, G. Lonzarich, A. Huxley, I. Sheikin, D. Braithwaite, and J. Flouquet, Nature 406, 587 (2000)
- (6) D. Aoki, A. Huxley, E. Ressouche, D. Braithwaite, J. Flouquet, J.-P. Brison, E. Lhotel, and C. Paulsen, Nature 413, 613 (2001)
- (7) P. Monthoux, D. Pines, and G. Lonzarich, Nature 450, 1177 (2007)
- (8) P. Monthoux and G. G. Lonzarich, Phys. Rev. B 59, 14598 (1999)
- (9) P. Monthoux and G. G. Lonzarich, Phys. Rev. B 63, 054529 (2001)
- (10) J. L. Sarrao and J. D. Thompson, J. Phys. Soc. Jpn. 76, 051013 (2007)
- (11) J. Thompson, R. Movshovich, Z. Fisk, F. Bouquet, N. Curro, R. Fisher, P. Hammel, H. Hegger, M. Hundley, M. Jaime, P. Pagliuso, C. Petrovic, N. Phillips, and J. Sarrao, J. Magn. Magn. Mater. 226-230, 5 (2001)
- (12) H. Shishido, R. Settai, D. Aoki, S. Ikeda, H. Nakawaki, N. Nakamura, T. Iizuka, Y. Inada, K. Sugiyama, T. Takeuchi, K. Kindo, T. C. Kobayashi, Y. Haga, H. Harima, Y. Aoki, T. Namiki, H. Sato, and Y. Ōnuki, J. Phys. Soc. Jpn. 71, 162 (2002)
- (13) G. F. Chen, K. Matsubayashi, S. Ban, K. Deguchi, and N. K. Sato, Phys. Rev. Lett. 97, 017005 (2006)
- (14) G. Knebel, D. Aoki, D. Braithwaite, B. Salce, and J. Flouquet, Phys. Rev. B 74, 020501 (2006)
- (15) G. Knebel, D. Aoki, J.-P. Brison, L. Howald, G. Lapertot, J. Panarin, S. Raymond, and J. Flouquet, phys. stat. sol. (b) 247, 557 (2010)
- (16) T. Park, M. J. Graf, L. Boulaevskii, J. L. Sarrao, and J. D. Thompson, Proc. Natl. Acad. Sci. U. S. A. 105, 6825 (2008)
- (17) S. Nair, O. Stockert, U. Witte, M. Nicklas, R. Schedler, K. Kiefer, J. D. Thompson, A. D. Bianchi, Z. Fisk, S. Wirth, and F. Steglich, Proc. Natl. Acad. Sci. U. S. A. 107, 9537 (2010)
- (18) R. R. Urbano, B.-L. Young, N. J. Curro, J. D. Thompson, L. D. Pham, and Z. Fisk, Phys. Rev. Lett. 99, 146402 (2007)
- (19) N. Kimura, K. Ito, K. Saitoh, Y. Umeda, H. Aoki, and T. Terashima, Phys. Rev. Lett. 95, 247004 (2005)
- (20) N. Fukushima, Phys. Rev. B 78, 115105 (2008)
- (21) P. D. Sacramento, J. Phys.: Condens. Matter 15, 6285 (2003)
- (22) H. Tsunetsugu, M. Sigrist, and K. Ueda, Rev. Mod. Phys. 69, 809 (1997)
- (23) J. V. Alvarez and F. Yndurain, Phys. Rev. Lett. 98, 126406 (2007)
- (24) Y. Yanase and M. Sigrist, J. Phys. Soc. Jpn. 78, 114715 (2009)
- (25) A. Aperis, G. Varelogiannis, and P. B. Littlewood, Phys. Rev. Lett. 104, 216403 (2010)
- (26) J. Spałek, Phys. Rev. B 38, 208 (1988)
- (27) M. M. Maśka, M. Mierzejewski, J. Kaczmarczyk, and J. Spałek, Phys. Rev. B 82, 054509 (2010)
- (28) For the discussion of various forms of - model see Ref. \rev@citealpnumJJ2 below.
- (29) K. A. Chao, J. Spalek, and A. M. Oles, J. Phys. C: Solid State Phys. 10, L271 (1977)
- (30) F. C. Zhang and T. M. Rice, Phys. Rev. B 37, 3759 (1988)
- (31) A. Moreno, A. Muramatsu, and S. R. Manmana, Phys. Rev. B 83, 205113 (2011)
- (32) K. Maki, Phys. Rev. 148, 362 (1966)
- (33) K. Kumagai, M. Saitoh, T. Oyaizu, Y. Furukawa, S. Takashima, M. Nohara, H. Takagi, and Y. Matsuda, Phys. Rev. Lett. 97, 227002 (2006)
- (34) P. W. Anderson, P. A. Lee, M. Randeria, T. M. Rice, N. Trivedi, and F. C. Zhang, J. Phys.: Cond. Matter 16, R755 (2004)
- (35) J. Bünemann, private communication
- (36) N. Fukushima, J. Phys. A: Math. Theor. 44, 075002 (2011)
- (37) P. Fazekas, Lecture Notes on Electron Correlation and Magnetism (World Scientific Publishing, 1999) chapter 7
- (38) M. Ogata and H. Fukuyama, Reports on Progress in Physics 71, 036501 (2008)
- (39) B. Kyung, Phys. Rev. B 62, 9083 (2000)
- (40) G. C. Psaltakis and E. W. Fenton, J. Phys. C 16, 3913 (1983)
- (41) S. Tsonis, P. Kotetes, G. Varelogiannis, and P. B. Littlewood, J. Phys.: Cond. Matter 20, 434234 (2008)
- (42) P. Fulde and R. A. Ferrell, Phys. Rev. 135, A550 (1964)
- (43) A. Larkin and Y. Ovchinnikov, J. Exp. Theor. Phys. 47, 1136 (1964)
- (44) A. Larkin and Y. Ovchinnikov, Sov. Phys. JETP 20, 762 (1965)
- (45) Y. Matsuda and H. Shimahara, J. Phys. Soc. Jpn. 76, 051005 (2007)
- (46) M. Ogata, J. Phys. Soc. Jpn. 66, 3375 (1997)
- (47) F. Loder, A. P. Kampf, and T. Kopp, Phys. Rev. B 81, 020511 (2010)
- (48) A. Ptok, M. M. Maśka, and M. Mierzejewski, J. Phys.: Cond. Matter 21, 295601 (2009)
- (49) J. Jȩdrak and J. Spałek, Phys. Rev. B 83, 104512 (2011)
- (50) J. Jȩdrak and J. Spałek, Phys. Rev. B 81, 073108 (2010)
- (51) J. Kaczmarczyk, J. Jȩdrak, and J. Spałek, Acta Phys. Polon. A 118, 261 (2010)
- (52) J. Jȩdrak, J. Kaczmarczyk, and J. Spałek, arXiv:1008.0021
- (53) Q.-H. Wang, Z. D. Wang, Y. Chen, and F. C. Zhang, Phys. Rev. B 73, 092507 (2006)
- (54) K.-Y. Yang, W. Q. Chen, T. M. Rice, M. Sigrist, and F.-C. Zhang, New J. Phys. 11, 055053 (2009)
- (55) J. Kaczmarczyk and J. Spałek, Phys. Rev. B 79, 214519 (2009)
- (56) J. Kaczmarczyk and J. Spałek, J. Phys.: Condens. Matter 22, 355702 (2010)
- (57) P. Korbel, J. Spałek, W. Wójcik, and M. Acquarone, Phys. Rev. B 52, R2213 (1995)
- (58) G. Kotliar and A. E. Ruckenstein, Phys. Rev. Lett. 57, 1362 (1986)
- (59) J. Spałek and W. Wójcik, Spectroscopy of the Mott Insulators and Correlated Metals, Vol. 119 (Springer-Verlag, Berlin, 1995) pp. 41–65
- (60) M. Galassi, J. Davies, J. Theiler, B. Gough, G. Jungman, P. Alken, M. Booth, and F. Rossi, GNU Scientific Library Reference Manual (3rd Ed.), ISBN 0954612078.
- (61) J. Kaczmarczyk, Ph. D. Thesis, Jagiellonian University, Kraków, 2011 http://th-www.if.uj.edu.pl/ztms/download/phdTheses/Jan_Kaczmarczyk_dokt%orat.pdf
- (62) T. Koretsune and M. Ogata, Phys. Rev. Lett. 89, 116401 (2002)
- (63) T. Koretsune and M. Ogata, J. Phys. Soc. Jpn. 72, 2437 (2003)
- (64) W. O. Putikka, M. U. Luchini, and M. Ogata, Phys. Rev. Lett. 69, 2288 (1992)
- (65) M. Raczkowski, D. Poilblanc, R. Frésard, and A. M. Oleś, Phys. Rev. B 75, 094505 (2007)
- (66) M. Ogata and A. Himeda, J. Phys. Soc. Jpn. 72, 374 (2003)
- (67) Although in CeRhSi the phases have only been analyzed as a function of temperature.
- (68) A. Himeda and M. Ogata, Phys. Rev. B 60, R9935 (1999)
- (69) T. Giamarchi and C. Lhuillier, Phys. Rev. B 43, 12943 (1991)
- (70) M. Inaba, H. Matsukawa, M. Saitoh, and H. Fukuyama, Physica C 257, 299 (1996)
- (71) O. Howczak, Ph. D. Thesis, Jagiellonian University, Kraków, in preparation