Coexistence of spin-triplet superconductivity with magnetism within a single mechanism for orbitally degenerate correlated electrons: Statistically-consistent Gutzwiller approximation

An orbitally degenerate two-band Hubbard model is analyzed with inclusion of the Hund's rule induced spin-triplet paired states and their coexistence with magnetic ordering. The so-called statistically consistent Gutzwiller approximation (SGA) has been applied to the case of a square lattice. The superconducting gaps, the magnetic moment, and the free energy are analyzed as a function of the Hund's rule coupling strength and the band filling. Also, the influence of the intersite hybridization on the stability of paired phases is discussed. In order to examine the effect of correlations the results are compared with those calculated earlier within the Hartree-Fock (HF) approximation combined with the Bardeen-Cooper-Schrieffer (BCS) approach. Significant differences between the two used methods (HF+BCS vs. SGA+real-space pairing) appear in the stability regions of the considered phases. Our results supplement the analysis of this canonical model used widely in the discussions of pure magnetic phases with the detailed elaboration of the stability of the spin-triplet superconducting states and the coexistent magnetic-superconducting states. At the end, we briefly discuss qualitatively the factors that need to be included for a detailed quantitative comparison with the corresponding experimental results.


Introduction
The question of coexistence of magnetism and superconductivity appears very often in correlated electron systems. In this context, both the spin-singlet and the spintriplet paired states should be considered. A general motivation for considering here the spin-triplet pairing is provided by the discoveries of superconductivity in Sr 2 RuO 4 [1,2], UGe 2 [3,4], URhGe [5], UIr [6], and UCoGe [7,8,9]. In the last four compounds, superconductivity indeed coexists with ferromagnetism. Moreover, for both, the spin-singlet high-temperature superconductors and the heavy-fermion systems, the antiferromagnetism and the superconductivity can have the same origin. Hence, it is natural to ask whether ferromagnetism and spin-triplet superconductivity also have the same origin in the itinerant uranium ferromagnets. A related and a very nontrivial question is concerned with the coexistence of antiferromagnetism with triplet superconducting state as in UNi 2 Al 3 [10,11,12] and UPt 3 [13,14].
It has been argued earlier [15,17,18,19] that for the case of indistinguishable fermions, the intra-atomic Hund's rule exchange can lead in a natural manner to the coexistence of spin-triplet superconductivity with magnetic ordering -ferromagnetism or antiferromagnetism in the simplest situations. This idea has been elaborated subsequently by us [20,21,22] by means of the combined Hartree-Fock(HF)-Bardeen-Cooper-Shrieffer(BCS) approach.
In particular, the phase diagrams have been determined which contain regions of stability of the pure superconducting phase of type A (i.e., the equal-spin-paired phase), as well as superconductivity coexisting with either ferromagnetism or antiferromagnetism.
The HF approximation, as a rule, overestimates the stability of phases with a broken symmetry. Therefore, in this work, we apply the Gutzwiller approximation for the same selection of phases in order to examine explicitly the effects of interelectronic correlations. The extension of the Gutzwiller method to the multi-band case [23,24,25] provides us with the so-called renormalization factors for our degenerate two-band models. With these factors we construct an effective Hamiltonian by means of the statistically consistent Gutzwiller approximation, SGA, in which additional constraints are added to the standard Gutzwiller approximation (GA) and with the incorporation of which the single-particle state has been determined (see [26,27,28,29] for exemplary applications of the SGA method). The detailed phase diagram and the corresponding order parameters are determined as functions of the microscopic parameters such as the band filling, n, the Hund's rule exchange integral, J, and the Hubbard interaction parameters, U and U ′ . The obtained results are compared with those coming out from the Hartree-Fock approximation. In this manner, the paper extends the discussion of itinerant magnetism within the canonical (extended Hubbard) model, appropriate for this purpose, to the analysis of pure and coexisting superconducting-magnetic states within a single unified approach. Additionally, at the end, we dwell briefly on the applicability of our original concepts to more realistic systems. It should be noted that theoretical investigations regarding the spin-triplet pairing have been performed recently also for other systems [30,31,32,33,34,35].
The paper is composed as follows. In Sec. II we provide the principal aspects of real-space spin-triplet pairing induced by the Hund's rule coupling, and introduce the band-renormalization factors for our two-band model. Furthermore, in subsections A and B of Sec. II we explain how the effective Hamiltonian is constructed, according to the statistically consistent Gutzwiller approximation, for all the phases considered in this work. In Sec. III we discuss the phase diagram, and the principal order parameters in the considered phases, whereas Sec. IV contains the concluding remarks and outlook concerning future investigations to make the approach applicable to real systems.

Model and method
We consider the extended orbitally-degenerate Hubbard Hamiltonian, which has the formĤ where l = 1, 2 label the orbitals and the first term describes electron hopping between atomic sites i and j. For l = l ′ this term represents electron hopping with change of the orbital (i.e., hybridization in momentum space). The next two terms describe the Coulomb interactions between electrons on the same atomic site. However the second term contains also the contribution, originating from the exchange interaction (J). The last term expresses the Hund's rule i.e., the ferromagnetic exchange between electrons localized on the same site, but on different orbitals. This term contributes to magnetic coupling and is responsible for the spin-triplet pairing leading to magnetic ordering, superconductivity and coexistent magnetic-superconducting phases. In the Hamiltonian (1), we have disregarded the pair hopping term because it hardly influences the ordered phases which we analyze in this work. In our variational method we assume that the correlated state |Ψ G of the system can be expressed in the following manner where |Ψ 0 is the normalized non-correlated state to be determined later andP G is the Gutzwiller correlator selected in the following form Here, λ (i) I,I ′ are the variational parameters, which are assumed to be real. In the two-band situation the local basis consists of 16 states (see Table 4), which are defined as follows where γ = 1, 2, 3, 4 labels the four spin-orbital states (in the lσ notation: 1 ↑, 1 ↓, 2 ↑, 2 ↓, respectively) and |I| is the number of electrons in the local state |I . In general, an index I can be interpreted as a set in the usual mathematical sense. The creation operators in (4) are placed in ascending order, i.e., γ 1 < ... < γ |I| . In an analogous manner, one can define the product of annihilation operatorŝ which are placed in descending order γ 1 > ... > γ |I| .
The operator |I ii I ′ | can be expressed in terms ofĈ † I andĈ I in the following mannerm In the subsequent discussion, we write expectation values with respect to |Ψ 0 as while the expectation values with respect to |Ψ G will be denoted by The most important step within the Gutzwiller approach is to derive the formula for the expectation value of the HamiltonianK =Ĥ − µN with respect to |Ψ G . This can be done in the limit of infinite dimensions by a diagrammatic approach [25] which uses the variational analog of Feynmann diagrams. By applying this method to the interaction part of the Hamiltonian (1), which is completely of intra-site character, one obtains whereĒ and L is the number of atomic sites. In (10) we have assumed that our system is homogeneous. Note that, with the use of Wick's theorem, the purely local expectation values m I 1 ,I 4 0 can be expressed in terms of the local single-particle density matrix elements ĉ α iγĉ α ′ iγ ′ 0 . Here,ĉ α iγ are either creation or annihilation operators. The expectation value of the single-particle part in the Hamiltonian (1) can be cast to the form where we have assumed that the renormalization factors q andq are real numbers and t γγ ′ = t γ ′ γ . Moreover, in the equation above we have neglected the part containing the inter-site pairing terms ĉ † i,γĉ † j,γ ′ 0 and ĉ i,γĉj,γ ′ 0 as we are going to concentrate on the Hund's rule induced intra-site spin-triplet paired states. The inter-site pairing amplitudes are much smaller than the intra-site terms, in the considered model. The renormalization factors, introduced in (12), have the form The parameters c * I 1 ,I 2 |γ in (13) are defined as where we introduced the fermionic sign function The renormalization factorsq γγ have to be included in (12), when there are nonzero gap parameters ( ĉ αĉα 0 = 0) in |Ψ 0 , which is the case considered here. The form ofq γγ is as follows The remaining part of K G that has to be derived is the expectation value N G . Also in this case, the diagrammatic evaluation in infinite dimensions gives the proper formula, where n iγ G = and The pairing densities in the correlated state that are going to be useful in the subsequent discussion can be expressed in the following way where Using (10), (12), and (18) one can express K G in terms of the variational parameters λ I,I ′ , local and non-local single particle density matrix elements , ĉ α iγĉ α ′ iγ ′ 0 , ĉ † i,γĉ j,γ ′ 0 , and the matrix elements of the atomic part of the atomic Hamiltonian represented in the local basis I|Ĥ at |I ′ .
The formula for K G , obtained in the way described above, can be written as an expectation value of an effective HamiltonianK GA , evaluated with respect to |Ψ 0 where q s γ = n iγ G / n iγ 0 . There is no guarantee that the condition is fulfilled. It turns out that it is fulfilled for the paramagnetic and the magnetically ordered phases of our two-band system, however it is not for the superconducting phases.
Physically it is most sensible to fix n G instead of n 0 , during the minimization. This is the reason why we include the term −µN already at the beginning of our derivation in K G . In this manner the chemical potential µ refers to the initial correlated system, not to the effective non-correlated one (for which the chemical potential can be different).
Having in mind that there are 16 states in the local basis there could be up to 16 × 16 = 256 variational parameters λ I,I ′ . However, for symmetry reasons many of these parameters are zero. The finite parameters can be identified by the following rule It should also be noted that, as shown in [25], the variational parameters are not independent since they have to obey the constrains which are going to be used to fix some of the parameters λ I,I ′ .
The results presented in this work have been obtained for the case of a square lattice with the band dispersions and also where β h ∈ [0, 1]. The orbital degeneracy and spatial homogeneity allow us to write Similar expressions as in (29) can be introduced for the expectation values in the noncorrelated state |Ψ 0 . Before discussing the principal magnetic and/or spin-triplet superconducting phases, we introduce first the exact expression of the full exchange operator (the last term of our Hamiltonian) via the local spin-triplet pairing operators (Â † im ,Â im ) namely We see that those two representations are mathematically equivalent, so the phase with S z G = Ŝ z il G = 0 and that with the corresponding off-diagonal order parameter should be treated on equal footing.

Statistically-consistent Gutzwiller method for superconducting and coexistent superconducting-ferromagnetic phases
In this subsection we will describe the SGA approach as applied to the selected phases characterized by the following order parameters • Superconducting phase of type A1 coexisting with ferromagnetism (A1+FM): G|u refers to the uniform magnetic moment and are the spin-triplet local gap parameters which are assumed as real here. The (correlated) order parameters which have been used above to define the relevant phases can also be defined for the non-correlated state |Ψ 0 . With these, we can determine which of the matrix elements m I,I ′ 0 are equal to zero for the considered phases. The assumption (25) then allows us to choose the non-diagonal variational parameters, λ I,I ′ , that have to be taken into account during the calculations. We list their indexes (I, I ′ ) in Table 2. Table 2. Nonzero, off-diagonal local variational parameters (λ I,I ′ = λ I ′ ,I ) that are used in the calculations for the considered phases. As one can see from Table 2, the off-diagonal variational parameters correspond to the creation or annihilation of the Cooper pair in the proper spin-triplet states |1 ↑, 2 ↑ and |1 ↓, 2 ↓ (phase A). Because in the A1 phase only electrons with spin-up are paired one can assume that λ 1,16 , λ 2,15 , λ 3,14 , λ 8,16 , λ 8,9 are zero (and their transoposed corespondantsλ I,I ′ = λ I ′ ,I ). For the FM and NS unpaired states only λ 10,11 and λ 11,10 are nonzero. They correspond to the two non-diagonal matrix elements of the atomic Hamiltonian, I|Ĥ at |I ′ . With the information contained in Table 2, one obtains the following relations regarding the band-narrowing renormalization factors where we have again used the γ = lσ notation. Due to the degeneracy of our bands we find Using the equations above we can rewrite the Hamiltonian (23) in the more explicit form, in reciprocal spacê where the renormalization factors Q σ are defined as Having the formula forK GA , given by (36), one can introduce next the so-called statistically-consistent Gutzwiller approximation (SGA). In this method, the mean fields (such as the expectation values for magnetization or superconducting gaps) are treated as variational mean-field order parameters with respect to which the energy of the system is minimized. However, in order to make sure that they coincide with the corresponding values calculated self-consistently, additional constraints have to be introduced with the help of the Lagrange-multiplier method [26,27,28,29]. This leads to supplementary terms in the effective Hamiltonian of the following form where the Lagrange multipliers λ m , λ s , and λ n are introduced to assure that the averages Â km , Ŝ kl and n klσ calculated either from the corresponding self-consistent equations or variationally, coincide with each other [29]. Introducing the four-component representation of single-particle operatorŝ we can write down the effective Hamiltonian in the following form whereM kσ is a 4x4 orthogonal matrix Here we introduced λ ↑ and λ ↓ which correspond to the Lagrange parameters λ m=1 and λ m=−1 , respectively. The bare quasiparticle energiesǫ klσ are defined as The diagonalization of the matrix (41) yields the quasiparticle eigen-energies in the paired states of the following form The first two energies correspond to the doubly degenerate spin-split quasiparticle excitations in the A phase, whereas the remaining two are their quasihole correspondents. Even though the Gutzwiller approach was derived for zero temperature, we may still construct the grand-potential function F λ (per atomic site) that corresponds to the effective Hamiltonian (40), i.e., The values of the mean fields, the variational parameters, and the Lagrange multipliers are found by minimizing the F λ functional, i.e., the necessary conditions for minimum are where A, Λ V , Λ L denote collectively the mean fields in the non-correlated state, the variational parameters and the Lagrange multipliers respectively. Additionally, the chemical potential, µ enters through the relation After solving the complete set of equations, one still has to calculate the mean fields in the correlated state with the use of their analogs in the non-correlated state and the variational parameters using (19) and (21).
With the SGA method one minimises the variational ground state energy K G with respect to the variational parameters λ I,I ′ and the single-particle states |Ψ 0 . Note that an alternative way for this minimization has been introduced, e.g., in [37]. Beyond the ground-state properties ofK one is often also interested in the (effective) single-particle Hamiltonian (40) because its eigenvalues are interpreted as quasi-particle excitation energies [38].

Statistically-consistent Gutzwiller method for the coexistent antiferromagnetic-spin-triplet superconducting phase
To consider antiferromagnetism in the simplest case, we divide our system into two interpenetrating sublattices A and B. In accordance with this division, we define the annihilation operators on the sublatticeŝ The same holds for the creation operators. Next, the Gutzwiller correlator can be expressed in the formP If we assume that charge ordering is absent, we have where the states I 1 and I 2 have opposite spins to those in the I 3 and I 4 states, respectively. The same division has to be made for the renormalization factors q,q and q s . They fulfill the transformation relations where γ andγ are spin-orbitals with opposite spins. The coexistent superconductingantiferromagnetic phase (SC+AF) can be defined in the following way Considerations analogical to those presented in subsection 2 lead to the conclusion that for both sublattices the non-diagonal variational parameters, λ A I,I ′ and λ B I,I ′ , that have to be used in the calculations, appropriate for the SC+AF phase, are the same as those listed in Table 2. This fact, and the degeneracy of our bands, allow us to apply (35) for both sets of renormalization factors (for A and B sublattices), as we have whereσ represents the spin opposite to σ. Now, we can write down the Hamiltonian K GA for the case of SC+AF phasê where It should be noted that the sums in (56) are taken over all L/2 independent k states. As before, we apply the SGA method which leads to the effective Hamiltonian with the statistical-consistency constraints of the form Introducing now the eight-component composite operator we can write down the effective HamiltonianK λ in the following form where the explicit form of the 8x8 matrix iŝ and Diagonalization of (61) leads to the quasi-particle energies E klσ (l = 1, 2, 3, ..., 8). The corresponding grand potential function F λ per atomic site now has the form As before, we minimize the F λ function to determine the values of the mean fields, the variational parameters and the Lagrange parameters. The necessary conditions for the minimum are again expressed by (45) and (46). In the subsequent discussion we consider also the pure antiferromagnetic phase (AF), for which S z G|s = 0 but ∆ + = ∆ − ≡ 0. The number of equations that need to be solved is different for different phases considered in this work. In Table III we show how many equations are included in (45) and (46) for all phases discussed.

Results and discussion
Equations (45) and (46) have been solved numerically for all phases by means of the hybrd1 subroutine from the MINPACK library, which performs a modification of the Powell hybrid method. The maximal estimated error of the procedure was set to 10 −7 .
The derivatives in Eq. (45) and (46) were computed by using a 5-step stencil method with the step equal to x = 10 −4 .
We concentrate now on the detailed numerical analysis of the phase diagram and the microscopic characterization of the stable phases. Having in mind that for 3d orbitals U ′ = U − 2J, one obtains the HF condition for the pairing to occur, U < 3J (see [22]). We discuss thus first and foremost the limit U < 3J, as it allows for a direct comparison of SGA with the HF solution. In this manner we can single out explicitly the role of correlations in stabilizing the relevant phases. One should note that in the considered regime (U < 3J) we have a model with intraatomic interorbital attractions leading to spin-triplet pairs. As the main attractive force is of intraatomic nature, we focus here on the local (s-wave) type of pairing only. In other words, as we discuss the situation with no or small hybridization, the intersite part of the pairing can be disregarded.
The calculations have been carried out for U ′ = U − 2J, U = 2.2J, k B T /W = 10 −4 . This leaves us still with three independent microscopic parameters in our model: n G , J, and β h . All the energies have been normalized to the bare band-width W = 8|t| (as we consider the square lattice with nearest neighbor hopping). For comparison, we also show the results calculated by means of the combined HF-BCS≡HF approximation. This method is described in detail in our previous paper for the same model as considered here. We can also reproduce the HF results by using the Gutzwiller method described in this work and setting λ I,I ′ = δ I,I ′ .
In Fig. 1 we display the free energy, superconducting gaps, and magnetic moments for the two values J = 0.299 and J = 0.4545. As one can see from the free-energy plots (Figs. 1a and 1d), below some certain value of band filling, the pure superconducting phase of type A is stable for the SGA method. The increase of the number of electrons in the system, enhances the gap in this region (Figs. 1b and 1e). Above the critical band filling n c , the staggered moment structure is created and a division into two gap parameters (∆ + and ∆ − ) appears, as can be seen in Figs. 1b, 1e, 1c, and 1f. In this regime the SC+AF phase becomes stable.
When approaching half filling, both gaps gradually approach zero and for n = 2 we are left with a pure AF phase, which is of Slater insulating type evolving towards the Mott-Hubbard insulating state with the increasing U. As the staggered magnetic moment is rising (with the increase of n G ), the renormalization factor is approaching unity (cf. Insets to Fig. 1a and 1d). This is a consequence of the fact that for large values of S z G , the configurations with two electrons of opposite spin, on the same orbital, are ruled out. Comparing Figs. 1a, 1b, 1c with Figs. 1d, 1e, 1f one sees that by increasing J we make the value of n c smaller. However, the decrease in n c is not as significant in SGA as it is in the HF case. In general the results presented Figs. 1b, 1c, 1e, 1f look similar from the qualitative point of view for both methods. For SGA, the onset of antiferromagnetically ordered phase appears closer to half filling than for the HF method. Another difference between HF and SGA is that for the former the staggered moment in the SC+AF phase is increased by the appearance of SC for the whole range of band fillings, whereas in SGA calculations the staggered moment is slightly stronger in the AF phase than in the SC+AF phase for a small region close to the half-filled situation (inset of Fig. 1f). Significant differences between HF and SGA can be seen in Figs. 1c and 1f. While changing the band filling from 0 to 2, in the case of SGA calculations we move consecutively through the regions of stability of NS (for J = 0.299), A, SC+AF phases, and for n = 2 we have pure antiferromagnetism. The situation is different in the HF approximation, where in between the regions of stability of A and SC+AF phase, we have also the stable A1+FM phase. It should be also noted that the free energy calculated in SGA is lower than the one for the HF situation, as one should expect, since the correlations are accounted more accurately in the former method. It is also   very interesting that having the system with U < 3J, no pure ferromagnetism appears in this canonical model of itinerant magnetism. In Fig. 2, we present the results for the case with nonzero hybridization parameter, β h = 0.1. In this case there are no superconducting solutions below some certain value of the band filling (cf. Fig. 2a) and an extended region of NS stability occurs. The influence of the hybridization on the antiferromagnetically ordered phases is weak, as can be seen more clearly in Fig. 3. The changes in the superconducting gap and the magnetic moment triggered by the hybridization, are quite small even for larger values of β h .
Next, we discuss the J dependence of the superconducting gap, the free energy and the magnetic moment for selected values of band filling. As in the case of n-dependences the gap parameters and the magnetic moments in both SGA and HF approximation are qualitatively similar. In Fig. 4 we can see that for n = 1.9 even the free-energy plots and regions of stability of certain phases are comparable for both calculation schemes used. For the quarter-filled case (cf. Fig. 5) the A1+FM phase is stable above some value of J, according to the HF results. However, this is not the case in the SGA approximation, where the A phase has lower free energy even than the saturated ferromagnetic phase coexisting with superconductivity. Comparing Figs. 5b and 5d (as well as 1d and 2b) one sees that the region of stability of the A phase narrows down in favor of the NS phase, due to the influence of hybridization.
It is important to check whether the itinerant magnetic phases are stable in the regime U ′ > J (U > 3J), i.e., when the superconductivity is absent in the HF approximation. For this purpose, in Fig. 6 we provide the band-filling dependence of the free energy corresponding to stable phases for U = 4J. Indeed, the paramagnetic and the magnetically ordered phases are stable for both methods of calculations. Therefore, for U > 3J we recover the magnetic phase diagram for this model, which was considered originally only in the context of magnetism. The free energy of the saturated ferromagnetic phase calculated by the SGA is very close to the one obtained by the HF approach. This is again caused by the circumstance that in the saturated state all of the spins are parallel and the double occupancies on the same orbital are absent. In this situation, the intra-orbital Coulomb interaction is automatically switched off. It would be interesting to determine the stability of the coexistent phases in this regime (U ′ > J). Work along this line is in progress.

Conclusions
The principal purpose of this paper was to formulate a many-particle method which allows to investigate the spin-triplet real-space pairing in correlated system with an orbital degeneracy. To this end, we have carried out a detailed analysis using the statistically-consistent Gutzwiller approximation (SGA) for the two-band degenerate Hubbard model with the spin-triplet superconductivity and itinerant magnetism included, both treated on equal footing. The results were compared with those coming from the Hartree-Fock approximation amended with the Bardeen-Cooper-Schrieffer (BCS) approach. The obtained Hund's coupling and band filling dependences of the magnetic moment and the superconducting gap parameters are often similar from the Table 4. Exemplary values of the order parameters, the chemical potential, the free energy, and the band renormalization factors corresponding to the considered phases, for two different sets of values of the microscopic parameters n and J. The underlined values correspond the stable phases. The numerical accuracy is on the level of the last digit specified. However, the stability regions of the considered phases are significantly different for the two applied methods. In SGA, the stable coexisting superconducting-ferromagnetic phase is absent while it appears in the HF approximation in a certain range of J and n values. Furthermore, the coexistence of the paired state with antiferromagnetism appears much closer to the half-filled situation in SGA than in HF approximation. For n = 2 the superconductivity disappears and only the pure antiferromagnetism survives; this state can be termed a correlated Slater-insulator state, which evolves gradually into the Mott-Hubbard insulating state with increasing U > 1 and S z G|s → 1/2. The influence of hybridization for both approximations is similar. With an increase of the β h parameter, the region of stability of the superconducting type-A phase narrows down in favor of the NS state. On the other hand, the antiferromagnetic phase is not affected in any significant manner by an increase of β h .
The band renormalization factors approach unity as the interaction constants J, U and U ′ tend to zero, what represents an additional test of our numerical results correctness. Generally, in the low-coupling limit our present results reduce to those obtained in HF approximation analysed by us in [22], as it should be.
It is important to emphasize that for both approaches the phase diagrams have been obtained for U < 3J, i.e. for relatively low value of the Hubbard interaction U, or equivalently, for a relatively high value of the Hund's rule exchange integral. A complete analysis of the present model would require studying the stability of the spin-triplet superconductivity and its coexistence with magnetic ordering in the complementary regime U > 3J, where the magnetism is favored against superconductivity. This regime has been the subject in a number of earlier papers [24,45,46], as then both the intraorbital, as well as the interorbital interaction is repulsive, and lead in a natural manner to magnetic ordering.

Outlook: Extension and application to real systems
In connection with the remarks provided above, we would like to characterize briefly the possibility of extending the present model (1) to the uranium systems in which superconductivity and ferromagnetism coexist in an unambiguous manner [36]. First of all, the magnetic moment in those systems, particularly in UGe 2 and URhGe, is quite large, with an associated molecular field in the megagauss range, which most probably rules out any spin-singlet character of pairing (note that the Curie temperature (T C ) to the superconducting transition temperature (T S ) ratio reaches in the uranium compounds the value T C /T S ∼ 10 2 ). In spite of those circumstances, our solution does not provide any extended regime (for the studied parameter range) for the ferromagnetism-spin-triplet superconductivity coexistence. Instead, in a wide range of band fillings, the coexisting SC+AF phase is stable (cf. Fig. 1c and 1f), as well as the pure spin-triplet superconducting phase of type A (the equal-spin-paired phase). The pure A phase seems to be realized in Sr 2 RuO 4 , though then a detailed three-orbital structure of the order parameter seems to be relevant [40]. A direct application of our SGA scheme to a realistic three-band system is more involved, as the number of parameters to minimize would lead to a computing time-consuming procedure, but still possible to tackle.
The extension of the present model to the uranium system such as UGe 2 would require considering orbitally degenerate and correlated 5f 2 − 5f 3 quasi-atomic states due to U and hybridized with the uncorrelated conduction band states. This means that we must have minimally a three-orbital system with two partially occupied 5f quasiatomic states (so the Hund's rule becomes operative) and at least one extra conduction band. Such situation may lead to a partial Mott-localization phenomenon, i.e., to a spontaneous decomposition of 5f n (n > 1) configuration of electrons into the localized and the itinerant parts [41]. In such a situation, it is possible that the localized electrons are the source of ferromagnetism, whereas the itinerant particles are paired [9]. This is not the type of coexistent phase we have in mind here, since in the model considered by us all the system electrons are indistinguishable in the quantum-mechanical sense.
These considerations lead to the conclusion that one would require minimally a periodic Anderson model with degenerate 5f n states, to mimic the uranium-based ferromagnetic superconductors. This variant of the multiple-band model is also very useful in the discussion of heavy-fermion compounds. Moreover, in the systems represented by this model, the coexistence of antiferromagnetism and superconductivity has been shown to appear in both experiment [48] and theory [49]. More specifically for the systems UPt 3 and UNi 2 Al 3 the coexistence of the spin-triplet pairing and the antiferromagnetism has already been suggested to appear [50,51], although not elaborated in any detail. One specific feature should be mentioned. Namely, in the situation when we have antiferromagnetic superconductor, then there is also a strong theoretical indication that there is a spin-triplet component even for the pure spin-singlet mechanism of pairing [27,49]. The spin-triplet important admixture results simply from a decomposition of the system into two sublattices with staggered magnetic moment. These and related features must be taken into quantitative analysis before any realistic consideration of concrete systems is carried out.
In connection with the whole discussion, it is intriguing to ask if a symmetric model system of the type exemplified by Hamiltonian (1) could be experimentally realized in the optical lattice. Some model systems (e.g., the Hubbard model system) have been experimentally achieved in this manner [42].
In relation to the spin-triplet real-space pairing induced by the Hund's rule, one should also mention the spin fluctuations (SF) as a possible mechanism of spin-triplet pairing in both magnetic [42] and liquid 3 He systems [43]. Within the present approach the spin fluctuations should be treated as quantum fluctuations around the present selfconsistently renormalized mean field state [44]. The real-space and the spin-fluctuation contributions may become of comparable magnitude in the close vicinity of the quantum critical point, where the ferro-or antiferro-states disappear under e.g. pressure. This is, however, a completely separate topic of studies.

Acknowledgements
Disscusions with Jakub Jędrak and Jan Kaczmarczyk are gratefully acknowledged. M.Z. has been partly supported by the EU Human Capital Operation Program, Polish Project No. POKL.04.0101-00-434/08-00. J.S. acknowledges the financial support from the Foundation for Polish Science (FNP) within project TEAM. The grant MAESTRO from the National Science Center (NCN) was helpful for the PL-DE cooperation within the present project on a unified approach to magnetism and superconductivity in correlated fermion systems.