text:
# Binary hard-sphere fluids near a hard wall

## I Introduction

A variety of experimental techniques has emerged which allow one to resolve the inhomogeneous density distributions of fluids at interfaces, a subject which enjoys broad scientific interest. In this context the ability to manufacture highly monodisperse colloidal suspensions has turned out to be particular useful as they provide the possibility to tune the effective interactions in these systems such that, e.g., the colloidal particles closely resemble hard-sphere fluids  . Since many of these experimental probes are indirect scattering techniques there is a substantial demand to guide them theoretically. Computer simulations and integral theories  are important tools of statistical physics to address these issues. Density functional theory (DFT)  has emerged as an additional approach which is capable to capture interfacial phase transitions and to sweep the thermodynamic and interaction parameter space of the system under consideration. The potential to combine these two possibilities poses already a major challenge for the other techniques. If DFT acquires in addition the same accuracy as the other two techniques, it could gain a clear competitive edge.

Although there is no recipe for constructing systematically a reliable DFT in spatial dimensions $`d2`$, the constant flux of developments over many years has led to a rather high level of sophistication. Among these theories for hard-sphere fluids, which act as paradigmatic systems and stepping-stones for more complicated models, the Rosenfeld functional has emerged as a particular powerful theory which resorts to the fundamental geometrical measures of the individual sphere . For the standard test case of the highly inhomogeneous density distribution of a one-component hard-sphere fluid near a planar hard wall, the predictions of the Rosenfeld functional are very close to those of numerical simulations serving as benchmarks. For this case the mean square deviations \[see, c.f., Eq. (22)\] of Rosenfeld DFT results from the simulation data from Ref.  are at most $`1\times 10^3`$ at high packing fractions, otherwise less than $`3\times 10^4`$.

Another virtue of the Rosenfeld functional is that is easily lends itself to the generalization to multi-component hard-sphere fluids. This opens the door to investigate rich new physical phenomena as particles of different size compete for interfacial positions . Even for the simplest multi-component system, the binary hard-sphere fluid, there are only relatively few theoretical studies which determine their structural properties near a planar hard wall, using Monte Carlo simulations , integral equation theories , and various kinds of density functional theory , as well as in spherical pores . Here we analyze this problem by using the corresponding Rosenfeld functional both in its original version  as well as for sophistications thereof . By comparing these results with published simulation data  we assess to which extent the quantitative reliability of the Rosenfeld functional for the one-component hard-sphere fluid remains valid for the corresponding binary system. Moreover we determine concentration profiles, the excess coverage, and the surface tension of the binary hard-sphere fluids at a hard wall. We describe the DFT in Sec. II and report our results in Sec. III followed by a summary and our conclusions in Sec. IV. The Appendix contains important technical details.

## II Density functional theory

The Rosenfeld functional for the excess (over the ideal gas) Helmholtz free energy of a mixture of hard spheres with number density profiles $`\{\rho _i(𝐫)\},i=1,\mathrm{},N`$, can be written as

$$\beta _{ex}[\{n_\alpha \}]=\mathrm{d}^3r\mathrm{\Phi }(\{n_\alpha (𝐫)\})$$
(1)
which is a functional of the four scalar weighted densities $`n_\alpha (𝐫)`$ for the $`N`$-component mixture

$$n_\alpha (𝐫)=\underset{i=1}{\overset{N}{}}\mathrm{d}^3r^{}\rho _i(𝐫^{})\omega _i^{(\alpha )}(𝐫𝐫^{}),\alpha =0,\mathrm{},3,$$
(2)
with $`4N`$ scalar weight functions $`\omega _i^{(\alpha )}`$ and two three-component vector weighted densities $`𝐧_\alpha (𝐫)`$

$$𝐧_\alpha (𝐫)=\underset{i=1}{\overset{N}{}}\mathrm{d}^3r^{}\rho _i(𝐫^{})𝝎_i^{(\alpha )}(𝐫𝐫^{}),\alpha =1,2,$$
(3)
with $`2N`$ vector weight functions $`𝝎`$$`{}_{}{}^{(\alpha )}{}_{i}{}^{}`$. The weight functions contain only information about the fundamental geometrical measures of a single sphere of species $`i`$, namely its volume, surface area, and radius $`R_i`$, i.e., in particular they are independent of the density profiles. The explicit expressions for the weight functions are given in the Appendix. $`\mathrm{\Phi }(\{n_\alpha \})=\mathrm{\Phi }_1+\mathrm{\Phi }_2+\mathrm{\Phi }_3`$ is a function of the weighted densities with

$`\mathrm{\Phi }_1`$ $`=`$ $`n_0\mathrm{log}(1n_3),`$ (4)
$`\mathrm{\Phi }_2`$ $`=`$ $`{\displaystyle \frac{n_1n_2𝐧_1𝐧_2}{1n_3}},`$ (5)
and

$`\mathrm{\Phi }_3`$ $`=`$ $`{\displaystyle \frac{\frac{1}{3}n_2^3n_2𝐧_2𝐧_2}{8\pi (1n_3)^2}}={\displaystyle \frac{n_2^3}{24\pi (1n_3)^2}}(13\xi ^2),`$ (6)
where $`\xi (𝐫)|𝐧_2(𝐫)|/n_2(𝐫)`$, which is the ratio of the modulus of the vector weighted density $`𝐧_2(𝐫)`$ and the scalar weighted density $`n_2(𝐫)`$. We want to note that $`\xi (𝐫)0`$ in the bulk and is small for small inhomogeneities. While this original Rosenfeld functional describes very successfully the fluid phase of a one-component hard-sphere system , it fails to predict the freezing transition. This failure has been studied in detail in Refs.  and . For the freezing transition it turns out that the zero-dimensional limit of the functional, in which a small cavity can accommodate only a single sphere, plays a key role. In a crystal the thermal vibrations around a lattice site can be interpreted as the motions in such a cavity formed by the neighboring spheres. Only if the statistical mechanics in such a cavity is described properly by the density functional, the freezing transition is predicted correctly. This is not the case for the original Rosenfeld functional. This problem can be fixed by modifying slightly the contribution $`\mathrm{\Phi }_3`$ \[see Eq. (6)\] such that the freezing transition is predicted by the modified functional while at the same time, for lower packing fractions, the accuracy of the original functional in describing the inhomogeneous fluid is kept. The following modifications have been suggested :

$$\mathrm{\Phi }_{3,q}=\frac{n_2^3}{24\pi }(1\xi ^2)^q$$
(7)
with $`q2`$ and

$$\mathrm{\Phi }_{3,int}=\frac{n_2^3}{24\pi }(13\xi ^2+2\xi ^3).$$
(8)
The first suggestion, $`\mathrm{\Phi }_{3,q}`$, is an antisymmetrized version of $`\mathrm{\Phi }_3`$ in Eq. (6) and the second, $`\mathrm{\Phi }_{3,int}`$, interpolates between $`\mathrm{\Phi }_3`$ of Eq. (6) and $`\mathrm{\Phi }_{3,0D}`$ in the exact zero-dimensional limit

$$\mathrm{\Phi }_{3,0D}=\frac{n_2^3}{24\pi (1n_3)^2}\xi (1\xi )^2.$$
(9)
While the modified Rosenfeld function with $`\mathrm{\Phi }_{3,0D}`$ does successfully predict the freezing transition of the one-component system, it leads to modified bulk properties and hence cannot describe the hard-sphere fluid as accurate as the original Rosenfeld functional. We note that the difference between $`\mathrm{\Phi }_3`$ of the original Rosenfeld functional and both $`\mathrm{\Phi }_{3,q}`$ with $`q=3`$ and $`\mathrm{\Phi }_{3,int}`$ is of the order of $`𝒪(\xi ^3)`$. Therefore we expect the biggest differences between the various versions of the Rosenfeld DFT to occur close to the wall where $`\xi `$ is largest.

Both the original Rosenfeld functional and the modifications corresponding to Eqs. (7) and (8), i.e., the functionals that share common bulk properties, are very successful and accurate for the one-component fluid. But far less is known for binary mixtures. While in this latter respect in Ref.  very good agreement between the density profiles obtained from the Rosenfeld functional and simulation data from Ref.  has been mentioned, in a recent study  significant deviations between the Rosenfeld DFT results and simulations have been found. In this latter study an alternative but equivalent formulation  of the original Rosenfeld functional has been applied.

Here we are interested in the equilibrium density profiles $`\rho _{s,0}(𝐫)`$ and $`\rho _{b,0}(𝐫)`$ of both the small and big components of binary hard-sphere mixtures close to a planar hard wall. To this end we freely minimize the functional

$$\mathrm{\Omega }[\rho _s(𝐫),\rho _b(𝐫)]=[\rho _s(𝐫),\rho _b(𝐫)]+\underset{i=s,b}{}\mathrm{d}^3r^{}\rho _i(𝐫^{})\left(V_{ext}^i(𝐫^{})\mu _i\right)$$
(10)
which is written in terms of the functional

$$[\rho _s(𝐫),\rho _b(𝐫)]=_{id}[\rho _s(𝐫),\rho _b(𝐫)]+_{ex}[\rho _s(𝐫),\rho _b(𝐫)]$$
(11)
with the exactly known ideal gas contribution $`_{id}`$,

$$\beta _{id}=\underset{i=s,b}{}\mathrm{d}^3r^{}\rho _i(𝐫^{})\left(\mathrm{ln}(\lambda _i^3\rho _i(𝐫^{}))1\right),$$
(12)
with $`\lambda _i`$ the thermal wave length of species $`i`$. For the equilibrium density profiles $`\rho _{i,0}(𝐫)`$, $`i=s,b`$, the functionals $``$ and $`\mathrm{\Omega }`$ reduce to the Helmholtz free energy and the grand canonical potential of the mixture, respectively; $`\mu _s`$ and $`\mu _b`$ are the chemical potentials of the two species. The external potentials entering into Eq. (10) model the planar hard wall at $`z=0`$:

$$V_{ext}^i(z)=\{\begin{array}{cc}\hfill \mathrm{},& z<R_i,\hfill \\ \hfill 0,& \text{otherwise},\hfill \end{array}$$
(13)
$`i=s,b`$, with $`z`$ the normal distance from the wall. The external potentials prevent the centers of spheres of species $`i`$ to approach the wall, located at $`z=0`$, closer than $`R_i`$ in which case they are in contact.

In the absence of spontaneous symmetry breaking due to freezing, which we do not consider here, the profiles $`\rho _{i,0}(z)`$, $`i=s,b`$, depend only on the normal coordinate $`z`$ which simplifies the minimization of the functional.

Far away from the wall, i.e., in the bulk system, the vector weighted densities $`𝐧_1`$ and $`𝐧_2`$ and thus $`\xi `$ vanish. In this limit both the original Rosenfeld functional and the two modifications corresponding to Eqs. (7) and (8) reduce to the same bulk expression given by

$$\mathrm{\Phi }_{bulk}=n_0\mathrm{ln}(1n_3)+\frac{n_1n_2}{1n_3}+\frac{n_2^3}{24\pi (1n_3)^2}$$
(14)
and hence they share the same bulk properties. We want to emphasize that as a consequence of this feature all versions of the Rosenfeld functional predict density profiles which show the same asymptotic decay towards the bulk value . The weighted densities in the bulk limit are obtained by inserting the bulk densities $`\rho _{i,bulk}:=\rho _{i,0}(z=\mathrm{})`$ into Eq. (2) yielding

$`n_3`$ $`=`$ $`{\displaystyle \frac{4\pi }{3}}{\displaystyle \underset{i=s,b}{}}R_i^3\rho _{i,bulk}{\displaystyle \underset{i=s,b}{}}\eta _i,`$ (15)
$`n_2`$ $`=`$ $`4\pi {\displaystyle \underset{i=s,b}{}}R_i^2\rho _{i,bulk},`$ (16)
$`n_1`$ $`=`$ $`{\displaystyle \underset{i=s,b}{}}R_i\rho _{i,bulk},`$ (17)
and

$`n_0`$ $`=`$ $`{\displaystyle \underset{i=s,b}{}}\rho _{i,bulk}.`$ (18)
The equation of state following from Eq. (14),

$$\beta p=\frac{n_0}{1n_3}+\frac{n_1n_2}{(1n_3)^2}+\frac{1}{12\pi }\frac{n_2^3}{(1n_3)^3},$$
(19)
is the Percus-Yevick compressibility equation of state of the mixture . This expression is related to the contact values of the density profiles according to the sum rule

$$\beta p=\underset{i=s,b}{}\rho _i(z=R_i+0).$$
(20)
This sum rule is respected by the Rosenfeld functional as by any weighted-density DFT  and therefore provides a test for the numerical accuracy of the calculations. In the following we suppress the subscript 0 which indicates equilibrium profiles as opposed to variational functions.

## III Structural and thermodynamic properties

### A Density profiles

The number density profiles of both components of binary hard-sphere mixtures close to a planar hard wall are obtained by a free minimization of the functional given in Eq. (10). We use the original Rosenfeld functional as well as the modified versions corresponding to Eq. (8) and to Eq. (7) with $`q=2`$ and $`q=3`$. The systems considered here have two different size ratios, $`R_b:R_s=5:3`$ and $`R_b:R_s=3:1`$, and various packing fractions $`\eta _s`$ and $`\eta _b`$ of the small and big spheres , respectively. The resulting density profiles are compared with simulation data published in Ref. . In addition we calculate the local concentrations $`\mathrm{\Phi }_s(z)`$ and $`\mathrm{\Phi }_b(z)`$ of the small and big spheres, respectively, defined as

$$\mathrm{\Phi }_i(z)=\frac{\rho _i(z)}{\rho _s(z)+\rho _b(z)},i=s,b.$$
(21)
We find excellent agreement between the density profiles of both components obtained by density functional theory and the simulation data for all systems under consideration. This holds for all versions of the Rosenfeld functional analyzed here. While at low total packing fractions $`\eta =\eta _s+\eta _b`$ the density functional theory results for all versions of the Rosenfeld functional are practically equivalent, small deviations among the results from different versions of the functional appear at larger values of $`\eta `$, i.e., for $`\eta 0.3`$. We quantify the degree of agreement between our DFT density profiles $`\rho _i(z)`$ and the simulation data from Ref. , available as data points $`(z_j,\rho _i^{sim}(z_j))`$, $`j=1\mathrm{}N_i^{sim}`$ and $`i=s,b`$, by determining the mean square deviations $`\overline{E}_i`$, $`i=s,b`$, defined as

$$\overline{E}_i=\frac{1}{N_i^{sim}}\underset{j=1}{\overset{N_i^{sim}}{}}\left(\frac{\rho _i^{sim}(z_j)\rho _i(z_j)}{\rho _{i,bulk}}\right)^2.$$
(22)
We find that $`\overline{E}_s`$ and $`\overline{E}_b`$ are at most $`5\times 10^4`$ and $`6\times 10^3`$, respectively, for all versions of the Rosenfeld functional. However, since the statistical errors in the simulation data are comparable with or even larger than the differences between the density profiles obtained by different version of the Rosenfeld DFT this approach does not enable us to determine which of the various versions is the most accurate one.

To illustrate the agreement between the DFT and the simulation data, in Fig. 1 we show the density profiles of the small spheres (a) and of the big spheres (b) for $`\eta _s=0.0607`$ and $`\eta _b=0.3105`$ and a size ratio $`R_b:R_s=5:3`$. The symbols ($`\mathrm{}`$) denote the simulation data from Ref.  and the solid line is obtained by the original Rosenfeld functional. The dotted lines denote coarse grained densities $`\overline{\rho }_i^{(j)}`$, $`i=s,b`$ and $`j=0,1`$, defined as

$$\overline{\rho }_i^{(j)}=\frac{1}{z_i^{(j+1)}z_i^{(j)}}\underset{z_i^{(j)}}{\overset{z_i^{(j+1)}}{}}dz\frac{\rho _i(z)}{\rho _{i,bulk}},$$
(23)
with $`z_i^{(0)}0`$ and $`z_i^{(j>0)}`$ the position of the $`j`$th minimum . All details of the density profiles found in the simulations are reproduced very accurately by the density functional theory. The oscillatory behavior, i.e., the amplitudes, the phases, and the decay of the oscillations as obtained by DFT agree excellently with the simulations. The total packing fraction of the system $`\eta =\eta _s+\eta _b=0.3712`$ is already rather high, giving rise to the pronounced structure of the density profiles. The corresponding concentration profiles $`\mathrm{\Phi }_s(z)`$ and $`\mathrm{\Phi }_b(z)`$ of the small and big spheres, respectively, are shown in Fig. 2. These concentration profiles demonstrate that, apart from the purely geometrical constraints, near the wall the big particles are enriched and the small particles depleted. This is in line with the expectation based on the attractive depletion potential near a hard wall of a single big sphere immersed in a fluid of small spheres . Correlation effects reverse this relative distribution in the second layer and restore it in the third.

As mentioned above, small differences between the DFT results corresponding to the various versions of the Rosenfeld functional can be found for these values of $`\eta `$. In order to be able to resolve these small differences magnified parts of the density profiles from Fig. 1 are shown in Fig. 3 together with the simulation data from Ref. ($`\mathrm{}`$). The solid lines in Fig. 3 correspond to the original Rosenfeld functional, the dotted lines correspond to the interpolated version \[Eq. (8)\] whereas the dashed and dashed-dotted lines correspond to the antisymmetrized version \[Eq. (7)\] with $`q=2`$ and $`q=3`$, respectively. All DFT results are very close to the simulation data. However, the deviations between the simulations and the antisymmetrized functional with $`q=2`$ seem to be systematically the largest.

In Figs. 4 and 5 we show the density profiles of a binary mixture with size ratio $`R_b:R_s=3:1`$. The packing fraction of the small spheres is $`\eta _s=0.0047`$ and that of the big spheres is $`\eta _b=0.3859`$ so that the total packing fraction $`\eta =\eta _s+\eta _b=0.3906`$ is again rather high. Therefore there is a strong spatial variation of the density profiles. The agreement between DFT (solid line) and simulations ($`\mathrm{}`$) is again found to be excellent for both the density profile of the small spheres (a) and of the big spheres (b). The dotted lines denote the coarse grained densities as defined in Eq. (23). The concentrations profiles of the small and big spheres, corresponding to these density profiles are shown in Fig. 6. For this larger size ratio the anti-correlated behavior of $`\mathrm{\Phi }_s(z)`$ and $`\mathrm{\Phi }_b(z)`$ is even more pronounced than for the smaller ratio discussed in Fig. 2 and strongly locked in without additional fine structure such as the double peak appearing in Fig. 2.

In addition we test the numerical accuracy of our calculations by means of the sum rule given in Eq. (20). In Table I the sum of the contact values of the binary mixture with size ratio $`R_b:R_s=5:3`$ for various packing fractions is compared with two equations of state. $`\beta p_{PY}^c`$ denotes the Percus-Yevick compressibility equation of state \[Eq. (19)\], to which the Rosenfeld functional reduces in the bulk limit, and $`\beta p_{MCSL}`$ corresponds to the more accurate Mansoori-Carnahan-Starling-Leland equation of state , which represents a generalization to a mixture of the very accurate Carnahan-Starling equation of state  for the one-component fluid. The very good agreement between the contact values and $`\beta p_{PY}^c`$ demonstrates the high accuracy of our numerical procedure. However, at higher packing fractions, $`\beta p_{PY}^c`$ deviates from the more accurate equation of state $`\beta p_{MCSL}`$. The same analysis of our results for a binary mixture with size ratio $`R_b:R_s=3:1`$ is summarized in Table II.

Equation (20) represents a sum rule which must be fulfilled by the density profiles as obtained by any of the density functionals considered here. However, no corresponding rules are available for the individual contact values. We find that for all versions of the Rosenfeld functional under consideration here, the sum rule is respected equally well. However, the individual contact values may differ. This statement is in line with the expectation that the biggest differences between the various versions of the Rosenfeld functional occur in a region where $`\xi `$ is large, i.e., close to the wall, and is substantiated in Table III for the binary mixture with size ratio $`R_b:R_s=5:3`$ and in Table IV for the size ratio $`R_b:R_s=3:1`$.

Vested with this confidence in our numerical procedures we are now able to comment on the comparison of our DFT results for the original Rosenfeld functional with those obtained by earlier DFT studies . We find that DFT approaches other than the Rosenfeld DFT predict the structure of the density profiles of a binary hard-sphere mixture near a planar hard wall only qualitatively , whereas the Rosenfeld functional yields quantitatively reliable predictions. This is demonstrated by Figs. 1 and 35 and the very small values of the mean square deviations $`\overline{E}_s`$ and $`\overline{E}_b`$ \[Eq. (22)\]. The predictions of the Rosenfeld functional agree excellently in all details with the simulation results from Ref. . The deviations between the DFT results of Ref. , calculated with the Rosenfeld functional, and their own simulations originate most likely from numerically problems of the iteration procedure used in Ref. . This suspicion is further supported by the fact that the DFT density profiles shown in Ref.  seem to violate the sum rule Eq. (20). Thus we are led to the conclusion that the deviations between the Rosenfeld DFT results and the simulation data reported in Ref.  most likely are artifacts generated by numerical problems in implementing the iteration procedure used in Ref.  for solving the Euler-Lagrange equations corresponding to the Rosenfeld density functional. Therefore the doubts raised in Ref.  about the performance of the Rosenfeld DFT for binary hard-sphere mixtures are not justified. We conclude that the Rosenfeld DFT exhibits the same high accuracy in predicting density distributions of binary hard-sphere mixtures as for the one-component hard-sphere fluid.

### B Excess adsorption and surface tension

One of the virtues of DFT is that based on the knowledge of the local structural properties $`\rho _i(z)`$, $`i=s,b`$, it is straightforward to calculate also thermodynamic properties such as the excess adsorptions $`\mathrm{\Gamma }_i`$ and the surface tension $`\gamma `$. Here we determine these quantities near a hard wall for a binary hard-sphere fluid whose components exhibits a size ratio $`R_b:R_s=3:1`$. Our analysis is confined to the fluid phase of the mixture; the phase boundary for freezing is estimated from the bulk phase diagrams presented in Refs. .

The excess adsorption of species $`i`$, $`i=s,b`$, is defined as

$`\mathrm{\Gamma }_i`$ $`=`$ $`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}dz(\rho _i(z)\rho _{i,bulk})`$ (24)
$`=`$ $`{\displaystyle \underset{\sigma _i/2}{\overset{\mathrm{}}{}}}dz(\rho _i(z)\rho _{i,bulk}){\displaystyle \frac{\sigma _i}{2}}\rho _{i,bulk}.`$ (25)
This definition of the excess adsorption differs from the definition used in Ref.  as well as from that used in Ref. . To recover the results for the excess adsorption of a one-component hard-sphere fluid in Ref.  and Ref.  one has to subtract from and add to our results, respectively, the constant $`\frac{\sigma _i}{2}\rho _{i,bulk}`$. These differences originate from different choices for the position of the wall. While for the one-component fluid there is no preference for any of these definitions, our choice used here appears to be particularly suited for a mixture because independent of the diameter $`\sigma _i`$ of species $`i`$ the integrals in Eq. (24) start at the same lower bound for all species, namely at the position of the physical wall. We use the same definition for the position of the wall to determine the surface tension. Because the excess adsorptions follow from integrating over oscillatory functions, they depend very sensitively on the precise structure of the density profiles and require very accurate calculations. Moreover, near the phase boundary for freezing the original Rosenfeld functional yields values for $`\mathrm{\Gamma }_s`$ and $`\mathrm{\Gamma }_b`$ which differ significantly from those obtained from the modifications of the Rosenfeld functional. Thus, unless stated otherwise, we have determined the excess adsorption by using the modified Rosenfeld functional corresponding to Eq. (8), which is known to capture the freezing transition of the one-component hard-sphere fluid.

In Fig. 7 we show the excess adsorptions $`\mathrm{\Gamma }_s`$ of the small spheres as function of the packing fractions $`\eta _s`$ and $`\eta _b`$; note that $`\mathrm{\Gamma }_s(\eta _s=0,\eta _b)0`$. For a fixed packing fraction of the small spheres $`\eta _s`$, the excess adsorption of the small spheres increases upon increasing $`\eta _b`$. The reason for this is, that the increasing packing effects of the big spheres associated with large values of $`\eta _b`$ enforces also the packing of the small spheres, giving rise to a strongly enhanced contact value and very pronounced structures in the density profile of the small spheres. For very small $`\eta _s`$, i.e., $`\eta _s0.1`$, and large $`\eta _b`$ we find that $`\mathrm{\Gamma }_s`$ can become positive. As function of $`\eta _b`$, $`\mathrm{\Gamma }_s`$ exhibits a turning point for any fixed value of $`\eta _s`$.

The excess adsorption $`\mathrm{\Gamma }_b`$ of the big spheres is shown in Fig. 8. For constant $`\eta _s`$ and increasing $`\eta _b`$, $`\mathrm{\Gamma }_b`$ decreases. But due to the same mechanism as described above, $`\mathrm{\Gamma }_b`$ increases upon increasing $`\eta _s`$ for constant packing fractions of the big spheres $`\eta _b`$. The square symbols ($`\mathrm{}`$) in Fig. 8 denote simulation data for the excess adsorption of a one-component hard-sphere fluid near a hard wall at packing fractions $`\eta _s`$=0.3680, 0.4103, and 0.4364, respectively, taken from Ref. . Whereas the simulation data for the two smaller packing fractions agree very well with our DFT results, that for the highest value of $`\eta _s`$ clearly deviates from the DFT prediction. Thus from this comparison it remains unclear to which extent the simulation data and the DFT results are reliable close to freezing.

In order to illustrate the large differences between the excess adsorptions calculated from different functionals, in Fig. 9 we show the excess adsorption $`\mathrm{\Gamma }_s`$ and $`\mathrm{\Gamma }_b`$ calculated from the original Rosenfeld functional together with those calculated by using the modified functional corresponding to Eq. (8). The packing fraction of the small spheres is $`\eta _s=0.15`$. Only for small packing fractions of the big spheres, i.e., far away from freezing, we find good agreement between the results obtained from the different functionals. Close to the phase boundaries  the differences become large, which is indeed surprising at first glance. However, the main contribution of the excess adsorption stems from the vicinity of the wall where $`\xi `$ is large and the different versions of the Rosenfeld functional are expected to differ. These differences were already indicated in the previous subsection by the behavior of the individual contact values $`\sigma _i^3\rho _i(R_i+0)`$, $`i=s,b`$, of the density profiles of the small and big spheres, respectively. For large distances from the wall the density profiles exhibit decaying oscillations and therefore are not expected to contribute essentially to the excess adsorption. Moreover, all versions of the Rosenfeld functional display a common characteristic decay because they share the same bulk properties  so that the contributions to the excess adsorption far away from the wall are very similar for the original Rosenfeld functional and its modifications.

The grand potential $`\mathrm{\Omega }`$ of a system in contact with a wall

$$\mathrm{\Omega }=\mathrm{\Omega }_{bulk}+\mathrm{\Omega }_{surf}$$
(26)
decomposes into a bulk contribution, $`\mathrm{\Omega }_{bulk}=pV`$, given by the bulk pressure $`p`$ in the system times the volume $`V`$ occupied by fluid particles, and a surface contribution, $`\mathrm{\Omega }_{surf}=\gamma A`$, which is the surface tension $`\gamma `$ times the surface area $`A`$ of the wall. Scaled particle theory (SPT) provides an approximate expression for $`\gamma `$ for a one-component hard-sphere fluid  as well as a generalization to hard-sphere mixtures  close to a planar hard wall. The surface tension of a one-component hard-sphere fluid within SPT is well tested and turns out to provide reliable results as compared with both DFT calculations  and simulations . In Ref.  a fit to simulation results of the surface tension of a one-component hard-sphere fluid at a planar hard wall is presented, which gives quasi-exact results and closely resembles the SPT expression.

In terms of the weighted bulk densities $`n_0,\mathrm{},n_3`$, defined in Eqs. (15)-(18), the SPT approximation for the surface tension of a hard-sphere mixture close to a planar hard wall  can be written as

$$\beta \gamma _{SPT}=\frac{n_1}{1n_3}+\frac{1}{8\pi }\frac{n_2^2}{(1n_3)^2}.$$
(27)
This expression reduces to the one-component SPT approximation of the surface tension when it is evaluated for the one-component bulk weighted densities. In this latter case the surface tension can be expressed solely in terms of the packing fraction $`\eta `$ of this single component.

Within the Rosenfeld functional the surface tension $`\gamma `$ of a binary hard-sphere mixture at a planar hard wall follows from the equilibrium density profiles $`\rho _i(z)`$, $`i=s,b`$, as obtained in the previous subsection:

$`\beta \gamma `$ $`=`$ $`{\displaystyle \frac{\beta \mathrm{\Omega }+\beta pV}{A}}`$ (28)
$`=`$ $`\underset{L\mathrm{}}{lim}\left[\beta pL+{\displaystyle \underset{0}{\overset{L}{}}}dz\left\{\mathrm{\Phi }[\rho _s(z),\rho _b(z)]+{\displaystyle \underset{i=s,b}{}}\rho _i(z)(V_{ext}^i(z)\mu _i)\right\}\right].`$ (29)
In contrast to the strong dependence of the results of the excess adsorption on the choice of functional, we find that the original Rosenfeld functional as well as its modifications predict very similar results for the surface tension in the whole range of packing fractions studied here. Our results for the surface tension of a binary hard-sphere mixture with size ratio $`R_b:R_s=3:1`$, calculated within the original Rosenfeld functional, are shown in Fig. 10. The deviation between these results and those obtained by the modifications of the Rosenfeld functional are at most 3% and deviate from the predictions of the SPT at most by 10%.

## IV Summary and conclusions

Based on the Rosenfeld density functional we have analyzed structural and thermodynamic properties of binary hard-sphere mixtures near a hard wall with the following main results:

1. Figures 1 and 4 demonstrate that the structure of both density profiles $`\rho _i(z)`$, $`i=s,b`$, of a binary hard-sphere mixture close to a planar hard wall as obtained by the original Rosenfeld functional is in excellent agreement with simulation results for size ratios $`R_b:R_s=5:3`$ and $`R_b:R_s=3:1`$, respectively. The high level of agreement between our DFT results and the simulation data from Ref.  is confirmed quantitatively by small mean square deviations $`\overline{E}_s`$ and $`\overline{E}_b`$ defined in Eq. (22). In terms of these quantities all versions of the Rosenfeld functional are practically of the same quality. Only at high packing fractions rather small differences between the results obtained by different versions of the Rosenfeld DFT become visible like those shown in Figs. 3 and 5.
2. The concentration profiles (see Fig. 2 and Fig. 6) calculated from the density profiles confirm the depletion picture: the small spheres are depleted from regions close to the wall while the big spheres are enriched.
3. The numerical accuracy of our calculations is demonstrated in Tables I and II by the high degree at which the sum rule Eq. (20), which relates the sum of the contact values of both density profiles with the equation of state, is respected. The sum rule, however, makes no prediction for the individual contact values and we find in Tables III and IV that each version of the Rosenfeld functional takes a different route to satisfy the sum rule.
4. Using the modified Rosenfeld functional corresponding to Eq.(8) we have calculated the excess adsorption of the small spheres $`\mathrm{\Gamma }_s(\eta _s,\eta _b)`$ (see Fig. 7) and of the big spheres $`\mathrm{\Gamma }_b(\eta _s,\eta _b)`$ (see Fig. 8) as functions of the packing fractions $`\eta _s`$ and $`\eta _b`$ for a binary hard-sphere mixture with size ratio $`R_b:R_s=3:1`$. We find that these quantities depend very sensitively on the accuracy of the numerical calculations and, as can be seen in Fig. 9, they differ significantly from the excess adsorption calculated by the original Rosenfeld functional.
5. The surface tension of a binary hard-sphere mixture with size ratio $`R_b:R_s=3:1`$ close to a planar hard wall is shown in Fig. 10. We find that all versions of the Rosenfeld functional give results which are in good agreement with the prediction of scaled particle theory \[Eq. (27)\].

From these results we conclude that the class of Rosenfeld functionals yields quantitatively reliable descriptions of interfacial structures in binary hard-sphere fluids. We expect that the same level of reliability also holds for multi-component hard-sphere fluids.

The excess adsorptions $`\mathrm{\Gamma }_s`$ and $`\mathrm{\Gamma }_b`$ of the small and big spheres emphasize the differences between the various versions of the Rosenfeld functional most. In order to decide whether the original Rosenfeld functional or whether its modifications predict these quantities more accurately, additional simulation data of the excess adsorption in a binary hard-sphere fluid are needed.

###### Acknowledgements.

It is a pleasure to thank Bob Evans for many stimulating discussions. We want to thank Dr. Noworyta for providing us with his simulation data.

## A Calculation of the weighted densities

Within the minimization procedure of the Rosenfeld functional the weighted densities $`n_\alpha `$ and $`𝐧_\alpha `$ have to be calculated repeatedly. Therefore it is necessary to optimize these calculations with respect to both computational speed and numerical accuracy.

The weight functions of the Rosenfeld functional are given by

$`\omega _i^{(3)}(𝐫)`$ $`=`$ $`\mathrm{\Theta }(|𝐫|R_i),`$ (A1)
$`\omega _i^{(2)}(𝐫)`$ $`=`$ $`\delta (|𝐫|R_i),`$ (A2)
and

$`𝝎_i^{(2)}(𝐫)`$ $`=`$ $`{\displaystyle \frac{𝐫}{|𝐫|}}\delta (|𝐫|R_i)`$ (A3)
with the Heaviside function $`\mathrm{\Theta }`$ and the Dirac delta function $`\delta `$. The remaining scalar weight functions are proportional to $`\omega _i^{(2)}`$: $`\omega _i^{(1)}=\omega _i^{(2)}/(4\pi R_i)`$ and $`\omega _i^{(0)}=\omega _i^{(2)}/(4\pi R_i^2)`$. The first vector weight function is collinear with $`𝝎`$$`{}_{}{}^{(2)}{}_{i}{}^{}`$: $`𝝎`$$`{}_{}{}^{(1)}{}_{i}{}^{}`$=$`𝝎`$$`{}_{i}{}^{(2)}/(4\pi R_i)`$.

In order to calculate the weighted densities integrals $`I_i^{(\alpha )}`$ of the type

$$I_i^{(\alpha )}=\mathrm{d}^3r^{}\rho _i(𝐫)\omega _i^\alpha (𝐫𝐫^{})$$
(A4)
have to be evaluated. For these convolution type integrals one can exploit the symmetry properties of the density profiles. For the present geometry the weighted densities can be written as

$$n_\alpha (z)=\underset{i=s,b}{}_{R_i}^{R_i}dz^{}\rho _i(z+z^{})\overline{\omega }_i^{(\alpha )}(z^{})$$
(A5)
with $`s`$ and $`b`$ for small and big, respectively, and with reduced weight functions $`\overline{\omega }_i^{(\alpha )}`$ which are functions of $`z`$ only:

$`\overline{\omega }_i^{(3)}(z)`$ $`=`$ $`\pi (R_i^2z^2),`$ (A6)
$`\overline{\omega }_i^{(2)}(z)`$ $`=`$ $`2\pi R_i,`$ (A7)
and

$`\overline{𝝎}_i^{(2)}(z)`$ $`=`$ $`2\pi z𝐞_z`$ (A8)
with the unit vector $`𝐞_z`$ in $`z`$-direction. The relations between these and the remaining weight functions are the same as for the original weight functions. The integrals in Eq. (A5) are one-dimensional convolutions which can be calculated faster and more accurate in Fourier space than in real space. By introducing the Fourier transforms of the density profiles,

$$\widehat{\rho }_i(k)=𝒯(\rho _i(z)),$$
(A9)
and those of the weight functions,

$$\widehat{\omega }_i^{(\alpha )}(k)=𝒯(\overline{\omega }_i^{(\alpha )}(z)),$$
(A10)
the weighted densities can be expressed as

$$n_\alpha (z)=𝒯^1\left(\underset{i=s,b}{}\widehat{\rho }_i(k)\widehat{\omega }_i^\alpha (k)\right).$$
(A11)
This route of calculation offers the important advantage that the numerical calculation of these convolutions can be speed up significantly by applying Fast-Fourier-Transform (FFT) methods. Moreover it turns out that calculations of convolutions in real space depend more sensitively on the grid size $`\mathrm{\Delta }z`$ to be used for discretization than those in Fourier space. We expect that the reason for this is that the FFT algorithm interpolates between data points with trigonometrical functions. To overcome this problem in real space a sophisticated integration scheme would have be applied or a very small grid size would have to be chosen. Both remedies additionally slow down the numerical calculation in real space.

The results presented in this appendix are applicable if the density profiles depend on the $`z`$ coordinate only. However, similar results can be obtained if the density profiles have radial symmetry.


format: text

text:
# Mechanism of spin-triplet superconductivity in Sr2RuO4

\[

## Abstract

The unique Fermi surfaces and their nesting properties of Sr2RuO4 are considered. The existence of unconventional superconductivity is shown microscopically, for the first time, from the magnetic interactions (due to nesting) and the phonon-mediated interactions. The odd-parity superconductivity is favored in the $`\alpha `$ and $`\beta `$ sheets of the Fermi surface, and the various superconductivities are possible in the $`\gamma `$ sheet. There are a number of possible odd-parity gaps, which include the gaps with nodes, the breaking of time-reversal symmetry and $`\stackrel{}{d}\widehat{z}`$.

\]

The nature of superconductivity discovered in Sr2RuO4 has been the subject of intense theoretical and experimental activity. Although Sr2RuO4 has the same layered perovskite structure as La2CuO4, the prototype of the high $`T_\mathrm{c}`$ cuprates superconductors, the electronic structures are very different and the nature of superconductivity seems to be totally different.

The normal state in Sr2RuO4 is characterized as essentially a Fermi liquid below $`50`$K. The resistivities in all directions show $`T^2`$ behavior for $`T50`$K. The effective mass is about $`34m_{\mathrm{electron}}`$ and the susceptibility is also about $`34\chi _0`$ where $`\chi _0`$ is the Pauli spin susceptibility. In contrast to the conventional normal state (below $`50`$K), there are considerable experimental evidences that the superconducting state (below about $`1.5`$K) is unconventional. The nuclear quadrupole resonance(NQR) does not show the Hebel-Slichter peak. The transition temperature is very sensitive to non-magnetic impurities. The 17O NMR Knight shift shows that the spin susceptibility has no change across $`T_\mathrm{c}`$ but stays just the same as in the normal state for the magnetic field parallel to the Sr2RuO4 plane . In addition, spontaneous appearance of an internal magnetic field below the transition temperature is reported by muon spin rotation measurements ($`\mu SR`$) .

Shortly after the discovery of the superconductivity in Sr2RuO4, it was proposed that the odd-parity(spin -triplet) Cooper pairs are formed in the superconducting state in analogy with 3He . The existence of ferromagnetic interaction is crucial in this proposal. In stead, incommensurate antiferromagnetic(AF) fluctuations were found by inelastic neutron scattering experiment . Ferromagnetic interaction is almost negligible compared with AF one.

Earlier specific heat and NQR measurements show a large residual density of states (DOS), $`5060`$% of DOS of the normal state, in the superconducting phase. A possible explanation, so called, orbital dependent superconductivity was proposed . Since four 4d electrons in Ru4+ partially fill the $`t_{2g}`$ band, the relevant orbitals are $`d_{xy}`$, $`d_{xz}`$ and $`d_{yz}`$ which determine the electronic bands. The gap of of one class of bands is substantially smaller than that of the other class of band. The presence of gapless excitations for temperatures greater than the smaller gap would account for the residual DOS. Recent specific measurements on high quality compounds, however, suggest the absence of residual DOS.

Sigrist et al.  proposed the following order parameter which is claimed to be compatible with most of the present experimental data,

$`\stackrel{}{d}=\widehat{z}(k_x\pm ik_y),`$ (1)
where $`\widehat{z}`$ is parallel to the $`\widehat{c}`$ axis and the gap is described as the tensor represented by $`\stackrel{}{d}`$ as $`\mathrm{\Delta }(k)=i(\stackrel{}{d}(k)\stackrel{}{\sigma })\sigma _y,`$ where $`\stackrel{}{\sigma }`$ is the Pauli matrix . Notice that the direction of the order parameter is frozen along the $`\widehat{c}`$ direction due to the crystal field and there is a full gap on the whole Fermi surface. The experiment of Josephson coupling between In and Sr2RuO4 also supports the gap directing $`\widehat{c}`$ axis . Based on the gap (1), the effects of impurity scattering, spin wave excitations and collective modes and sound propagation are studied theoretically.

However, the recent experiments on high quality compounds of Ru NQR and specific heat strongly suggest the existence of nodes and the absence of residual DOS. It is proposed a number of gap order parameters consistent with the existence of nodes phenomenologically . See also Ref..

Details of the Fermi surface have been observed by quantum oscillations . The Fermi surface consists of three sheets, which is consistent with the electronic band calculations . The Fermi sheets are labeled by $`\alpha `$, $`\beta `$, and $`\gamma `$. See Fig. 1. While the $`\gamma `$ sheet of the Fermi surface can be attributed solely to $`d_{xy}`$ Wannier function, the $`\alpha `$ and $`\beta `$ sheets are due to the hybridization of the $`d_{xz}`$ and $`d_{yz}`$ Wannier functions. The $`\gamma `$ band is quasi-isotropic two-dimensional, on the other hand the $`\alpha `$ and $`\beta `$ sheets are quasi-one dimensional which can be visualized as a set of parallel planes separated by $`Q=4\pi /3`$ running both $`k_x`$ and $`k_y`$ directions. Therefore it is natural to expect a sizable nesting effects at the wave vectors $`(\pm Q,k_y,k_z)`$ and ($`k_x,\pm Q,k_z`$). The nesting vectors $`(\pm Q,\pm Q,k_z)`$ and $`(Q,\pm Q,k_z)`$ have the maximum effects since they connect the one dimensional Fermi surfaces in both directions. Collective modes in the spin dynamics are studied based on these nesting effects . In fact the neutron scattering experiment shows peaks at $`(\pm 0.6\pi ,\pm 0.6\pi )`$ close to the nesting vectors (up to $`(\pm 2\pi ,\pm 2\pi )`$) .

– Band structure and the pairing interactions
We denote the annihilation operators for electrons in the three $`4d`$-$`t_{2g}`$ orbitals $`d_{xz}`$, $`d_{yz}`$ and $`d_{xy}`$ of Ru-ions as $`a_{k,s}`$, $`b_{k,s}`$ and $`c_{k,s}`$. The kinetic term of the $`\alpha `$ and $`\beta `$ bands is given by

$`H_{\mathrm{kin}}^{\alpha \beta }`$ $`=`$ $`{\displaystyle \underset{k,s}{}}(\epsilon ^{\alpha \beta }(k_x)\mu )a_{k,s}^{}a_{k,s}`$ (2)
$`+`$ $`{\displaystyle \underset{k,s}{}}(\epsilon ^{\alpha \beta }(k_y)\mu )b_{k,s}^{}b_{k,s}`$ (3)
$`+`$ $`{\displaystyle \underset{k,s}{}}t(k)a_{k,s}^{}b_{k,s}+{\displaystyle \underset{k,s}{}}t^{}(k)b_{k,s}^{}a_{k,s},`$ (4)
and that of the $`\gamma `$ band is given by

$`H_{\mathrm{kin}}^\gamma ={\displaystyle \underset{k,s}{}}(\epsilon ^\gamma (k_x,k_y)\mu )c_{k,s}^{}c_{k,s}.`$ (5)
The $`\alpha `$ and $`\beta `$ bands are quasi one-dimensional and the Fermi surfaces are well-approximated by four sheets $`k_x\pm Q/2`$ and $`k_y\pm Q/2`$. The mixing coefficient $`t(k)`$ can be neglected except around $`k_x=\pm k_y`$. The $`\gamma `$ band is quasi-isotropic two-dimensional one. The mixing term between the $`\alpha `$ or $`\beta `$ band and the $`\gamma `$ one is suppressed by the reflection symmetry of $`z`$.

The nesting of the $`\alpha `$ and $`\beta `$ bands leads to the following AF fluctuations:

$`H_{\mathrm{AF}}^{\alpha \beta }`$ (6)
$`=`$ $`{\displaystyle \underset{k,q,s_i}{}}J_{}^{\alpha \beta }(q_x)(\sigma _z)_{s_1s_3}(\sigma _z)_{s_2s_4}a_{k,s_1}^{}a_{k,s_2}^{}a_{k+q,s_3}a_{kq,s_4}`$ (7)
$`+`$ $`{\displaystyle \underset{k,q,s_i}{}}J_{}^{\alpha \beta }(q_y)(\sigma _z)_{s_1s_3}(\sigma _z)_{s_2s_4}b_{k,s_1}^{}b_{k,s_2}^{}b_{k+q,s_3}b_{kq,s_4}`$ (8)
$`+`$ $`{\displaystyle \underset{k,q,s_i}{}}J_{}^{\alpha \beta }(q_x)\left[(\sigma _x)_{s_1s_3}(\sigma _x)_{s_2s_4}+(\sigma _y)_{s_1s_3}(\sigma _y)_{s_2s_4}\right]`$ (10)
$`\times a_{k,s_1}^{}a_{k,s_2}^{}a_{k+q,s_3}a_{kq,s_4}`$
$`+`$ $`{\displaystyle \underset{k,q,s_i}{}}J_{}^{\alpha \beta }(q_y)\left[(\sigma _x)_{s_1s_3}(\sigma _x)_{s_2s_4}+(\sigma _y)_{s_1s_3}(\sigma _y)_{s_2s_4}\right]`$ (12)
$`\times b_{k,s_1}^{}b_{k,s_2}^{}b_{k+q,s_3}b_{kq,s_4},`$
where $`J_{}^{\alpha \beta }(q_i)>0`$ and $`J_{}^{\alpha \beta }(q_i)>0`$ have peaks at $`q_i\pm Q`$. In general, $`J_{}^{\alpha \beta }(q_i)J_{}^{\alpha \beta }(q_i)`$. The magnetic field generated by the AF fluctuations above induces the interaction $`H_{\mathrm{AF}}^\gamma `$ in the $`\gamma `$ band:

$`H_{\mathrm{AF}}^\gamma `$ (13)
$`=`$ $`{\displaystyle \underset{k,q,s_i}{}}J_{}^\gamma (q_x,q_y)(\sigma _z)_{s_1s_3}(\sigma _z)_{s_2s_4}c_{k,s_1}^{}c_{k,s_2}^{}c_{k+q,s_3}c_{kq,s_3}`$ (14)
$`+`$ $`{\displaystyle \underset{k,q,s_i}{}}J_{}^\gamma (q_x,q_y)\left[(\sigma _x)_{s_1s_3}(\sigma _x)_{s_2s_4}+(\sigma _y)_{s_1s_3}(\sigma _y)_{s_2s_4}\right]`$ (16)
$`\times c_{k,s_1}^{}c_{k,s_2}^{}c_{k+q,s_3}c_{kq,s_3},`$
where $`J_{}^\gamma (q_x,q_y)>0`$ and $`J_{}^\gamma (q_x,q_y)>0`$ have peaks at $`(q_x,q_y)=(\pm Q,\pm Q)`$ and $`(q_x,q_y)=(Q,\pm Q)`$.

In addition to the AF fluctuations, we consider the electron-phonon interaction. Due to the low dimensionality of the system, this interaction is weakly screened. We assume that the phonon-mediated interaction for the $`\alpha `$ and $`\beta `$ bands $`H_{\mathrm{ph}}^{\alpha \beta }`$ is quasi-one dimensional,

$`H_{\mathrm{ph}}^{\alpha \beta }`$ $`=`$ $`{\displaystyle \underset{k,k^{},s,s^{}}{}}f^{\alpha \beta }(q_x)a_{k,s}^{}a_{k,s^{}}^{}a_{k+q,s}a_{kq,s^{}}`$ (17)
$`+`$ $`{\displaystyle \underset{k,k^{},s,s^{}}{}}f^{\alpha \beta }(q_y)b_{k,s}^{}b_{k,s^{}}^{}b_{k+q,s}b_{kq,s^{}},`$ (18)
where $`f^{\alpha \beta }(q_i)>0`$ has a peak at $`q_i=0`$. For the $`\gamma `$ band, the phonon-mediated interaction we consider is

$`H_{\mathrm{ph}}^\gamma ={\displaystyle \underset{k,k^{},s,s^{}}{}}f^\gamma (q_x,q_y)c_{k,s}^{}c_{k,s^{}}^{}c_{k+q,s}c_{kq,s^{}}`$ (19)
where $`f^\gamma (q_x,q_y)>0`$ has a peak at $`(q_x,q_y)=(0,0)`$.

– Superconductivity in the $`\alpha `$ and $`\beta `$ bands
First, we examine the superconductivity in the $`\alpha `$ and $`\beta `$ bands. As we will show immediately, the odd-parity superconductivity is realized due to the quasi-one dimensionality . For simplicity, we put $`f^{\alpha \beta }(q_i)f^{\alpha \beta }(0)\delta _{q_i,0}`$, $`J_{}^{\alpha \beta }(q_i)J_{}^{\alpha \beta }(Q)\delta _{q_i,\pm Q}`$ and $`J_{}^{\alpha \beta }(q_i)J_{}^{\alpha \beta }(Q)\delta _{q_i,\pm Q}`$ in the following.

In the lowest order approximation, we neglect the mixing term $`t(k)`$, so the gap equation is separated for $`a_{k,s}`$ and $`b_{k,s}`$. For $`a_{k,s}`$ electrons, the gap $`\mathrm{\Delta }^{(a)}`$ is defined by

$`\mathrm{\Delta }_{s_2,s_1}^{(a)}(k)={\displaystyle \underset{k^{},s_3,s_4}{}}V_{s_1,s_2,s_3,s_4}(k_x,k_x^{})a_{k^{},s_3}a_{k^{},s_3},`$ (20)
where

$`V_{s_1,s_2,s_3,s_4}(k_x,k_x^{})`$ (21)
$`=f^{\alpha \beta }(k_x+k_x^{})\delta _{s_1s_3}\delta _{s_2s_4}f^{\alpha \beta }(k_xk_x^{})\delta _{s_1s_4}\delta _{s_2s_3}`$ (22)
$`+J_{}^{\alpha \beta }(k_x+k_x^{})(\sigma _z)_{s_1s_3}(\sigma _z)_{s_2s_4}`$ (23)
$`J_{}^{\alpha \beta }(k_xk_x^{})(\sigma _z)_{s_1s_4}(\sigma _z)_{s_2s_3}.`$ (24)
$`+J_{}^{\alpha \beta }(k_x+k_x^{})\left[(\sigma _x)_{s_1s_3}(\sigma _x)_{s_2s_4}+(\sigma _y)_{s_1s_3}(\sigma _y)_{s_2s_4}\right]`$ (25)
$`J_{}^{\alpha \beta }(k_xk_x^{})\left[(\sigma _x)_{s_1s_4}(\sigma _x)_{s_2s_3}+(\sigma _y)_{s_1s_4}(\sigma _y)_{s_2s_3}\right].`$ (26)
If the gap is unitary (as we assume), the gap equation becomes

$`\mathrm{\Delta }_{s_2,s_1}^{(a)}(k)`$ $`={\displaystyle \underset{k^{},s_3,s_4}{}}V_{s_1,s_2,s_3,s_4}(k_x,k_x^{})`$ (28)
$`\times {\displaystyle \frac{\mathrm{\Delta }_{s_3,s_4}^{(a)}(k^{})}{2E_k^{}^{(a)}}}\mathrm{tanh}\left({\displaystyle \frac{\beta E_k^{}^{(a)}}{2}}\right),`$
where $`E_k^{(a)}=\sqrt{ϵ^{\alpha \beta }(k_x)+\mathrm{tr}(\mathrm{\Delta }^{(a)}\mathrm{\Delta }^{(a)})/2}`$. The gap $`\mathrm{\Delta }^{(a)}`$ is parameterized as

$`\mathrm{\Delta }^{(a)}(k)=\{\begin{array}{cc}i\psi ^{(a)}(k)\sigma _y\hfill & \text{for even-parity gap}\hfill \\ i(\stackrel{}{d}^{(a)}(k)\stackrel{}{\sigma })\sigma _y\hfill & \text{for odd-parity gap}\hfill \end{array},`$ (31)
and we obtain the solutions as $`\psi ^{(a)}|_{k_x\pm Q/2}=\mathrm{const}.`$ for even-parity gap and $`\stackrel{}{d}^{(a)}|_{k_x\pm Q/2}=\pm \mathrm{const}.`$ for odd-parity gap. For these solutions, the critical temperature is $`k_\mathrm{B}T_\mathrm{c}=1.13\mathrm{}\omega _\mathrm{D}e^{1/N(\mu )\lambda },`$ where $`\omega _\mathrm{D}`$ is the cut-off for the interactions, $`N(\mu )`$ is the DOS at the Fermi surface for $`a_{k,s}`$ electrons and $`\lambda `$ is

$`\lambda =\{\begin{array}{cc}f^{\alpha \beta }(0)J_{}^{\alpha \beta }(Q)2J_{}^{\alpha \beta }(Q)\hfill & \text{for even-parity gap}\hfill \\ f^{\alpha \beta }(0)+J_{}^{\alpha \beta }(Q)2J_{}^{\alpha \beta }(Q)\hfill & \text{for }\stackrel{}{d}^{(a)}\stackrel{}{z}\hfill \\ f^{\alpha \beta }(0)J_{}^{\alpha \beta }(Q)\hfill & \text{for }\stackrel{}{d}^{(a)}\stackrel{}{z}\hfill \end{array}.`$ (35)
Therefore, the odd-parity superconductivity is realized in $`a_{k,s}`$, and $`\stackrel{}{d}^{(a)}`$ vector becomes

$`\stackrel{}{d}^{(a)}|_{kyQ/2}=\stackrel{}{d}^{(a)}|_{kyQ/2}=\mathrm{const}.,`$ (36)
and

$`\begin{array}{cc}\stackrel{}{d}^{(a)}\widehat{z}\hfill & \text{if }J_{}^{\alpha \beta }(Q)>J_{}^{\alpha \beta }(Q)\text{ }\hfill \\ \stackrel{}{d}^{(a)}\widehat{z}\hfill & \text{if }J_{}^{\alpha \beta }(Q)<J_{}^{\alpha \beta }(Q)\text{ }\hfill \end{array}.`$ (39)
Note that if $`J_{}^{\alpha \beta }(Q)`$ is large enough, the superconductivity is realized even when $`f^{\alpha \beta }(0)=0`$.

In a similar manner, the superconductivity in $`b_{k,s}`$ is odd-parity and characterized by (36) and (39) when every $`a`$’s, $`k_x`$ and $`k_y`$ are replaced by $`b`$’s, $`k_y`$ and $`k_x`$. The critical temperature is same as $`a_{k,s}`$. Due to $`D_{4h}`$ symmetry, it holds that $`|\stackrel{}{d}^{(a)}|=|\stackrel{}{d}^{(b)}|`$.

The mixing term $`t(k)`$ affects the gaps $`\stackrel{}{d}^{(a)}`$ and $`\stackrel{}{d}^{(b)}`$ at $`k=\pm K(\pm Q/2,\pm Q/2,k_z)`$ and $`k=\pm \stackrel{~}{K}(Q/2,\pm Q/2,k_z)`$. If $`t(k)`$ is large enough, the Fermi surface at $`k=\pm K`$ and $`k=\pm \stackrel{~}{K}`$ is largely deformed. Because of the weakness of the screening, the attractive force for the electron at those points decreases. This causes nodes at those points. The relative phase of $`\stackrel{}{d}^{(a)}`$ and $`\stackrel{}{d}^{(b)}`$ is not determined in this case, so the time-reversal symmetry is broken in general.

If $`t(k)`$ is small, $`\stackrel{}{d}^{(a)}`$ and $`\stackrel{}{d}^{(b)}`$ do not disappear at $`k=\pm K`$ and $`k=\pm \stackrel{~}{K}`$. In this case, the relative phase of $`\stackrel{}{d}^{(a)}`$ and $`\stackrel{}{d}^{(b)}`$ is determined. Since $`t(k)`$ at $`k_x=\pm k_y`$ is real, it can be easily shown that $`\mathrm{arg}\stackrel{}{d}^{(a)}=\mathrm{arg}\stackrel{}{d}^{(b)}`$ or $`\mathrm{arg}\stackrel{}{d}^{(a)}=\mathrm{arg}\stackrel{}{d}^{(b)}+\pi `$ at $`k=\pm K`$ and $`k=\pm \stackrel{~}{K}`$. The superconductivity in this case is illustrated in Fig. 2. The existence of nodes in the small $`t(k)`$ case depends on the symmetry of the system. If the system is invariant under the following transformation,

$`\begin{array}{cc}a_{k,s}b_{k,s},b_{k,s}a_{k,s}\hfill & \text{for }k_x=\pm k_y\hfill \\ a_{k,s}a_{k,s},b_{k,s}b_{k,s}\hfill & \text{for }k_x\pm k_y\hfill \end{array},`$ (42)
the following pairing interaction is added.

$`\mathrm{\Delta }H_{\mathrm{int}}^{\alpha \beta }=`$ $`{\displaystyle \underset{s,s^{}}{}}\mathrm{\Delta }_{ss^{}}^{(a)}(K)b_{K,s}b_{K,s}+(\text{h.c.})`$ (45)
$`{\displaystyle \underset{s,s^{}}{}}\mathrm{\Delta }_{ss^{}}^{(a)}(\stackrel{~}{K})b_{\stackrel{~}{K},s}b_{\stackrel{~}{K},s}+(\text{h.c.})`$
$`+(ab).`$
This interaction changes the spectrum at $`k=\pm K`$ and $`k=\pm \stackrel{~}{K}`$. If $`\stackrel{}{d}^{(a)}\stackrel{}{d}^{(b)}`$ holds, nodes exists at either $`k=\pm K`$ or $`k=\pm \stackrel{~}{K}`$. If the system does not have the symmetry above or $`\stackrel{}{d}^{(a)}\parallel ̸\stackrel{}{d}^{(b)}`$, no node exists in the gap.

– Superconductivity in the $`\gamma `$ band
Next, we examine the superconductivity in the $`\gamma `$ band based on the phonon-mediated interaction and AF fluctuations. In contrast to the $`\alpha `$ and $`\beta `$ bands, both even-parity and odd-parity superconductivity can be realized. If $`J_{}^\gamma (q_x,q_y)`$ is large enough, the odd-parity superconductivity is also realized in the $`\gamma `$ band, but if not, the various superconductivity can be realized. To show this, we introduce the following AF interaction and phonon-mediated interaction:

$`J_{\left\{\genfrac{}{}{0pt}{}{}{}\right\}}^\gamma (q_x,q_y)`$ (46)
$`={\displaystyle \underset{Q_{\left\{\genfrac{}{}{0pt}{}{x}{y}\right\}}=\pm 4\pi /3}{}}{\displaystyle \frac{g_{\left\{\genfrac{}{}{0pt}{}{}{}\right\}}}{(q_xQ_x)^2+(q_yQ_y)^2+Q_{\mathrm{AF}}^2}},`$ (47)
and

$`f^\gamma (q_x,q_y)={\displaystyle \frac{g_{\mathrm{ph}}}{q_x^2+q_y^2+Q_{\mathrm{ph}}^2}},`$ (48)
where $`g_{\left\{\genfrac{}{}{0pt}{}{}{}\right\}}`$, $`g_{\mathrm{ph}}`$, $`Q_{\mathrm{ph}}`$ are constants. From the experimental data in , we take $`Q_{\mathrm{AF}}=0.0797\times 2\pi `$.

To solve the gap equation, we approximate the Fermi surface of the $`\gamma `$ band as a cylinder with $`k_\mathrm{F}=2\pi /3`$ and use the mean-field approximation. We have solved the gap equation analytically when the phonon-mediated interaction dominates. The gap in this case is given by

$`\psi (k)\{\begin{array}{cc}\mathrm{const}.\hfill & \text{for }s\text{ gap}\hfill \\ k_x^2k_y^2\hfill & \text{for }d\text{ gap}\hfill \end{array},`$ (51)
$`\stackrel{}{d}(k)(k_x\pm ik_y)\text{for }p\text{ gap}.`$ (52)
In Fig.3, we show the phase diagram in the case where the phonon-mediated interaction dominates and $`g_{}=g_{}`$.

– Discussions
To conclude, we have studied the superconductive properties of Sr2RuO4based on AF fluctuations and phonon-mediated interactions. Due to the quasi-one dimensionality, the odd-parity superconductivity is realized in the $`\alpha `$ and $`\beta `$ bands. The existence of nodes in the gaps of these bands depends on the symmetry of the system (see Eq.(42)) or the strength of the hybridization of $`d_{xz}`$ and $`d_{yz}`$. The direction of $`\stackrel{}{d}`$ vector is determined by the anisotropy of the AF fluctuations.

The $`\gamma `$ band has many possibilities of superconductivity. If the AF fluctuations direct to the $`\widehat{c}`$ axis and are strong enough, the odd-parity superconductivity with $`\stackrel{}{d}\widehat{z}`$ is realized, however, if not, $`s`$\- or $`d`$-wave superconductivity is also possible.

There are various possibilities to explain the experimental results of Sr2RuO4 from our theory. Here we give one of them which is not in the previous theories. The NMR data  and the experiment of the Josephson coupling between In and Sr2RuO4  support the odd-parity superconductivity with $`\stackrel{}{d}\widehat{z}`$. This is easily realized in our theory with large $`J_{}^{\alpha \beta }(Q)`$ and/or $`J_{}^\gamma (Q,Q)`$. The nodes suggested by NQR and specific heat are likely to be those in the $`\alpha `$ and $`\beta `$ bands. The square vortex lattice observed by the neutron scattering  is consistent with the odd-parity superconductivity in the $`\alpha `$ and $`\beta `$ bands since the superconductivity occurs in the orthogonal quasi-one dimensional systems. The $`\mu `$SR data  can be explained by the superconductivity in the $`\alpha `$ and $`\beta `$ bands with the large hybridization of $`d_{xz}`$ and $`d_{yz}`$ or/and the superconductivity in $`\gamma `$ bands with the gap (52). Since larger $`J_{}^{\alpha \beta }(Q)`$ induces larger $`J_{}^\gamma (Q,Q)`$, it is likely that the odd-parity superconductivity in the $`\alpha `$ and $`\beta `$ band is followed by that in $`\gamma `$ band, so the absence of residual DOS also can be explained. The transition temperatures could be different between $`\gamma `$ and $`\alpha `$, $`\beta `$ Fermi surfaces, but there will be a single $`T_\mathrm{c}`$ if the hybridization between $`d_{xy}`$ and $`d_{yz}`$,$`d_{xz}`$ on the Ru atoms is large enough.

Acknowledgments

It is a pleasure to thank P. Bourges, Y. Hasegawa, K. Ishida, H.-Y. Kee, Y. Maeno, K. Maki, Y. Matsuda, Y. Sidis, and M. Sigrist, for useful discussions.