text:
# Dependence of hadronic properties on Quark Masses and Constraints on their Cosmological Variation

## I Introduction

The understanding of how exactly the fundamental parameters of the Standard Model enter any observable is certainly one of the most important aims of hadronic/nuclear physics. Two old deep questions drives its discussion: (i) Can there be “alternative Universes” with different set of parameters, and what are the boundaries of the World we know in the parameter space?; (ii) How to observe cosmological variations of weak and strong scales?

Discussion of both issues has been significantly revived recently. We will not go into questions (i) (see e.g. ) and only mention the latter (ii). The issue of cosmological time variation of major constants of physics has been recently revived by astronomical data which seem to suggest a variation of electromagnetic $`\alpha `$ at the $`10^5`$ level for the time scale 10 bn years, see . The statistical significance of the effect at the moment obviously exclude any random fluctuations, so the effect definitely exists. Whether it may or may not have a conventional explanation is not yet clear: more experimental work is clearly needed to reach any conclusions.

Nevertheless, it is quite timely to have another look at existing limits on time variation of all the fundamental constants. In particular, since the electromagnetic and weak forces are mixed together in the Standard Model, one may expect a similar modification of the weak couplings, the weak scale in general and quark masses in particular. In fact, one can measure only variation of dimensionless parameters. Therefore, we obtain limits on variation of $`m_s/\mathrm{\Lambda }_{QCD}`$ where $`m_s`$ is the strange quark mass and $`\mathrm{\Lambda }_{QCD}`$ is the QCD scale defined as a position of Landau pole in logarithm for running coupling constant. It is convenient to put $`\mathrm{\Lambda }_{QCD}=const`$.

A generic further argument goes as follows. The masses of 3 heavy quarks – c,b,t – are too large to be important in hadronic and nuclear physics. The masses of two light ones – u,d – are important, in particular via the pion contribution to nuclear forces we studied in our previous paper  and a subsequent one . Much more extensive discussion of this issue in a contexts of chiral perturbation theory can be found in literature, see e.g.  and references therein. The conclusion reached in those studies is that the $`m_u,m_d`$ are too small to be really important.

Thus, we focus on the dependence on the $`strange`$ quark mass $`m_s`$, the only one which has the right magnitude to generate a maximal sensitivity of the hadronic/nuclear physics to the weak scale. Indeed, as follows from QCD phenomenology, its variation from 0 to experimental value influence vacuum parameters such as the quark condensates $`\overline{q}q`$ at the factor-2 level. It also affects the masses of even $`nonstrange`$ hadrons such as the nucleon at about 20 percent level, etc. The fundamental explanation of why such unexpected “strangeness” is in fact present is related with the important role of the instanton-induced effects in QCD. As it follows from the chiral anomaly relation and is explained in details e.g. in a review , the multi-fermion ’t Hooft interaction necessarily involves all 3 flavors at all steps, even in the interactions of light u,d quarks.

The main new element of the present work is identification of the specific hadron, the $`\sigma `$, now known as $`f_0(600)`$ meson in the Review of Particle Properties , as the important ingredient of the nuclear physics most sensitive to $`m_s`$. This resonance is known for a very long time as a structure in $`\pi \pi `$ scattering, and is widely used in nuclear physics, see e.g. a review on applications of Walecka model . And nevertheless, its acceptance as a real hadronic state was difficult and had a very complicated history. Non-relativistic quark models of hadrons intensify it is as a l=1 or p-state, and tend to predict its mass to be around 1.4 GeV or so, in striking contradiction to experiment which finds a mass of only about $`m_\sigma 500MeV`$. Also large width of $`\sigma `$ makes this state to be easily deformed by all kind of effects, see e.g. recent work by one of us predicting drastic modification on the sigma shape in pp and heavy ion collisions . On the other hand, recent data have elucidated sigma production in a set of much simpler situations, from heavy quark hadrons, $`\mathrm{{\rm Y}}`$ transitions and $`D`$ decays. Those consistently point toward the smaller width $`\mathrm{\Gamma }_\sigma 250MeV`$. This reduced the controversy and put this particle back into the Review of Particle Properties .

In the next section we will argue that the sigma mass has stronger sensitivity to $`m_s`$, than that of ordinary non-strange hadrons like $`N`$ or $`\omega `$. This happens because of a valence strange part plus the repulsion from the nearby $`\overline{K}K,\eta \eta `$ continuum thresholds. We then estimate the derivative of the deuteron binding and neutron resonance energies to $`m_s`$.

## II Sensitivity of hadronic masses to variation of quark masses

### A Why sigma?

The very first appearance of the sigma mesons was as a two-pion scalar-isoscalar resonance. It has been gradually learned that the corresponding channel for $`\overline{q}q`$ interaction is the “maximally attractive channel”, with attraction so strong that it breaks spontaneously chiral symmetry and produce the non-zero quark condensate. The mechanism of that attraction is attributed mostly to instanton-induced ’t Hooft interaction, for a review see .

Sigma meson is an excitation on top of the scalar condensate, a kind of a Higgs boson of strong interactions. If indeed one naively assumes that it underlines all hadronic masses, e.g. that of the nucleon, the corresponding coupling can be estimated as

$`g_s=M_n/f_\pi 10`$ (1)
. This large value in turn implies that the perturbation theory can only be used as a qualitative guide, at best.

In passing, let us mention that arguments about development of the most optimum effective description of hadronic/nuclear physics in terms of mesonic degrees of freedom, known also as Quantum Hadro Dynamics, are still going on. Using some field variables can be better than others: in this respect let us mention the paper  which emphasized the instead of the traditional $`\sigma `$ of the linear sigma-model, a chiral partner of the pion, one better use the radial field $`\sqrt{\sigma ^2+\stackrel{}{\pi }^2}`$ which has normal derivative coupling to pions. There is extensive literature on loop corrections and related observables, such as resonance mass and shape modification in nuclear and excited hadronic matter, see e.g.  as a recent example.

For the purposes of this work it would be sufficient to use simple and widely used Walecka model, which keeps only the sigma and the omega exchanges in the effective nuclear forces

$`V={\displaystyle \frac{g_s^2}{4\pi }}{\displaystyle \frac{e^{rm_\sigma }}{r}}+{\displaystyle \frac{g_v^2}{4\pi }}{\displaystyle \frac{e^{rm_\omega }}{r}}`$ (2)
The very important lesson about nuclear forces this model had emphasized is that the nuclear potential is in fact a highly tuned small difference of two large terms.

We will argue in the next section that there are reasons to think that the sensitivity of these two terms to the fundamental weak scale is quite different. Sigma (scalars) involve all quark flavors while omega (vectors) do not, forbidden by Zweig rule. As a result, there is no fine tuning in the derivative over $`m_s`$, which significantly enhances the effect to be derived.

The values of the two coupling constants used in the nuclear matter applications of this model  are

$`g_s^2=357.4m_\sigma ^2/m_N^2100`$ (3)
$`g_v^2=273.8m_\omega ^2/m_N^2190`$ (4)
for $`m_\sigma =500MeV`$. Note that the effective scalar coupling is close to naive value (1) mentioned above.

### B  Strange valence and strange sea of the $`\sigma `$ meson

Scalar and pseudoscalar mesons are different from more familiar vector and axial ones in terms of their flavor composition. In the latter case the so called Zweig rule applies, forbidding flavor mixing: so for example the $`\omega `$ meson we will discuss in his work has a truly negligible mixing with a strange counterpart, $`\varphi `$ meson. Scalar and pseudoscalar mesons are on the contrary nearly ideally $`SU(3)`$ octets and singlets, so different flavors are very strongly mixed together. The pseudoscalar channel has been studied extensively, and we know that in this case $`\eta ^{}`$ is much heavier than (twice more strange) $`\eta `$ meson. This is the famous Weinberg U(1) problem, resolved by existence of the instanton-induced repulsive interaction pushing the singlet state upward.

The same instanton-induced interaction has the opposite sign in the $`SU(3)`$ singlet scalar channel, pushing $`m_\sigma `$ downward, see  for details. The magnitude of flavor mixing matrix element in the scalar channel has been also evaluated on the lattice, and the results agree with the instanton-based predictions in sign and magnitude.

Apart of theoretical motivation, there is a simple phenomenological fact that no purely strange $`\overline{s}s`$ counterpart to sigma resonance $`f_0(600)`$ is seen. There are strong evidence that the pair of states $`f_0(980),a_0(980)`$ are the $`\overline{K}K`$ molecule, so the next $`f_0`$ are at 1300 and 1500 MeV.

For those reasons, we think that a description of $`\sigma `$ as a $`SU(3)`$ singlet state

$`\sigma ={\displaystyle \frac{1}{\sqrt{3}}}(\overline{u}u+\overline{d}d+\overline{s}s)`$ (5)
strongly split from the $`SU(3)`$ octet one is a reasonable approximation. This means that the valence contribution to the derivative is

$`{\displaystyle \frac{m_\sigma }{m_s}}|_{val}=<\sigma |\overline{s}s|\sigma >=2/3`$ (6)
since with probability 1/3 there is a strange $`pair`$. A mixing between different scalar mesons $`f_0`$ ($`\sigma f_0(600)`$, $`f_0(980)`$, $`f_0(1370)`$) would further change the valence contributions downward. However, based on large gap between sigma and other states, and also based on better studied $`\eta \eta ^{}`$ mixing, one might think that the relative change due to mixing with next $`f_0`$ states is not significant.

Let us now consider the contribution of the so called strange sea, virtual $`\overline{s}s`$ pairs, which are always present even in a completely non-strange hadrons like a nucleon. The relatively well studied case is the nucleon mass sensitivity to $`m_s`$, already discussed in . Let briefly remind that $`KN`$ scattering data imply that

$`{\displaystyle \frac{m_N}{m_s}}=<N|\overline{s}s|N>1.5`$ (7)
and thus about 1/5 of the nucleon mass comes from the strange sea ($`m_s`$ 120 MeV).

Similar matrix elements for $`\sigma ,\omega `$ mesons is not possible to obtain experimentally, although it can be done on the lattice. To estimate it we will adopt a simple constituent quark picture, assuming additivity of the strange sea. If so, the derivatives analogous to (7) for all mesons should be $`2/3`$ of it, or

$`{\displaystyle \frac{m_{mesons}}{m_s}}|_{sea}1`$ (8)
As we will see later, the exact value of the common sea contribution is not actually important, since its contributions to $`\sigma `$ and $`\omega `$ mesons tend to cancel each other nearly exactly when we calculate variation of $`NN`$ interaction. What matters is the $`difference`$ between their strange seas, to be discussed in the next section.

It is convenient to present the effect of the possible quark mass variation on the $`\sigma `$ mass in the following form

$`{\displaystyle \frac{\delta m_\sigma }{m_\sigma }}=0.4({\displaystyle \frac{\delta m_s}{m_s}}+{\displaystyle \frac{m_u+m_d}{m_s}}{\displaystyle \frac{\delta m_q}{m_q}})=0.4{\displaystyle \frac{\delta m_s}{m_s}}+0.04{\displaystyle \frac{\delta m_q}{m_q}}`$ (9)
where we have used $`m_s=120MeV`$. We see that the relative change of the strange quark mass produces much larger effect than the relative change of the light quark mass. This is similar to the case of the nucleon mass variation

$`{\displaystyle \frac{\delta m_N}{m_N}}=0.19{\displaystyle \frac{\delta m_s}{m_s}}+0.045{\displaystyle \frac{\delta m_q}{m_q}}`$ (10)
### C $`\overline{K}K,\eta \eta `$ mixing with $`\sigma `$

As the second approximation, we will discuss loop effects, or mixing with 2-meson states. Those are also completely different for $`\sigma `$ and $`\omega `$ mesons. As we mentioned earlier, the latter practically does not mix with $`\overline{K}K`$ states, while $`\sigma `$ does mix with them strongly. An admixture of virtual $`\overline{K}K,\eta \eta `$ pairs can be viewed as an additional contribution to the strange sea, on top of the strange content of the non-strange constituent quarks (8).

The $`\sigma `$ mixing with continuum of pseudoscalars is described by the standard mass operator given by the usual loop diagram

$`\mathrm{\Sigma }(Q^2)={\displaystyle \frac{d^4k}{2\pi ^4}\frac{\lambda ^2(k,Q)}{[(k+Q/2)^2m^2+iϵ][(kQ/2)^2m^2+iϵ]}}`$ (11)
where $`\lambda `$ is the $`KK\sigma ,\eta \eta \sigma `$ or $`\pi \pi \sigma `$ couplings. Its real part describes the shift of $`m_\sigma `$ due to repulsion from these states: the imaginary part (for pions only) gives the width. The sign of the shift due to $`\overline{K}K`$ is obviously negative, since $`m_\sigma <2m_K`$. The effect of $`\pi \pi `$ has contributions of both signs. For constant coupling the shift is logarithmically divergent, in reality it has to be regulated by formfactors in the vertexes. The total shift is negative and large, of the order of very large width $`\mathrm{\Gamma }_\sigma 300MeV`$. Note that this large negative shift is partly the reason why the $`\sigma `$ mass is so small. However, in this paper we focus on the dependence on quark masses. Assuming that the main dependence comes from masses of Goldstone bosons, $`m_\pi ,m_K`$, we differentiate the mass operator over these masses and obtain the convergent result. Thus one can ignore formfactors and extract the effective coupling constant out of the integral.

For the derivative at $`Q=m_\sigma =0.5GeV`$ we get the following numerical value for the shifts

$`{\displaystyle \frac{\mathrm{\Sigma }}{m_K^2}}=0.0229GeV^2\lambda _{\sigma KK}^2`$ (12)
$`{\displaystyle \frac{\mathrm{\Sigma }}{m_\eta ^2}}=0.019GeV^2\lambda _{\sigma \eta \eta }^2`$ (13)
The couplings are not experimentally known, so we rely on the SU(3) symmetry and relate them to $`\lambda _{\sigma \pi \pi }`$, which is in turn related to sigma meson width

$`\mathrm{\Gamma }_\sigma ={\displaystyle \frac{3}{2}}{\displaystyle \frac{\lambda _{\sigma \pi \pi }^2}{16\pi m_\sigma }}\sqrt{1{\displaystyle \frac{4m_\pi ^2}{m_\sigma ^2}}}`$ (14)
Taking $`\mathrm{\Gamma }_\sigma `$=250 MeV we obtain $`\lambda _{\sigma \pi \pi }^2`$= 5 Gev2. The factor 3/2 account for $`\pi ^+\pi ^{},\pi ^0\pi ^0`$ modes. However in the mass shift there are contribution of $`K^+K^{},\overline{K}^0K^0,\eta \eta `$ channels which we would count as 5/2. Substituting numbers and using standard Gell-Mann-Oaks-Renner expressions for $`m_\eta ,m_K`$ ($`m^2m_s\mathrm{\Lambda }_{QCD}`$), we obtain the additional sensitivity of the sigma mass shift arising from the mixing effects

$`{\displaystyle \frac{\delta m_\sigma }{m_\sigma }}={\displaystyle \frac{\delta m_s}{m_s}}\left[2{\displaystyle \frac{m_K^2}{2m_\sigma ^2}}{\displaystyle \frac{\mathrm{\Sigma }}{m_K^2}}+0.5{\displaystyle \frac{m_\eta ^2}{2m_\sigma ^2}}{\displaystyle \frac{\mathrm{\Sigma }}{m_\eta ^2}}\right]{\displaystyle \frac{\delta m_s}{m_s}}0.14`$ (15)
We included $`K^+K^{},\overline{K}^0K^0`$ modes with the coefficient 1 and $`\eta \eta `$ with 1/2: the latter only contribute about 1/5 of the final answer. We ignored even smaller contribution of the $`\eta ^{}\eta ^{}`$ loop.

### D The total sensitivity of $`m_\sigma `$ as compared to $`m_N`$

Together with the one estimated in the previous section, it leads to total

$`{\displaystyle \frac{\delta m_\sigma }{m_\sigma }}(0.24+0.16+0.14){\displaystyle \frac{\delta m_s}{m_s}}=0.54{\displaystyle \frac{\delta m_s}{m_s}}`$ (16)
where 3 terms are the contributions of the common strange sea (8), valence strangeness (6) and the loop mixing (15), respectively.

In the same units the sensitivity of the nucleon mass is

$`{\displaystyle \frac{\delta m_N}{m_N}}0.19{\displaystyle \frac{\delta m_s}{m_s}}`$ (17)
We conclude that the sigma mass is about 3 times more sensitive to the variation of the strange quark mass than the nucleon mass.

We also need the sensitivity of $`\omega `$ meson to the strange mass variation. This meson does not have valence strange quarks and practically does not mix with $`\varphi `$, $`K`$ and $`\eta `$ mesons. Therefore, only the strange see contributes:

$`{\displaystyle \frac{\delta m_\omega }{m_\omega }}0.15{\displaystyle \frac{\delta m_s}{m_s}}`$ (18)
## III The modification of the deuteron binding

### A Preliminary analytic estimates

Simple analytic estimate for sensitivity of the deuteron binding to sigma and omega mass modification is obtained by the differentiation of the potential over the mass and averaging the resulting expression $`exp(mr)`$ over the radial wave function.

The simplest short range approximation leads to the free motion wave function

$`\psi (r)={\displaystyle \frac{\sqrt{2\kappa }}{r}}exp(\kappa r)`$ (19)
This wave function tends to infinity at small distances. However, the real deuteron wave function should be small there because of the repulsive core in the potential $`V(r)`$. Therefore, we introduced a small cut-off radius $`b`$ in the integration over $`r`$. We estimated $`b`$ from the condition $`V(b)=0`$ which gives $`b=0.45`$ fm. We get the following shift of the deuteron binding

$`{\displaystyle \frac{\delta Q_d}{Q_d}}={\displaystyle \frac{m_\sigma }{Q_d}}{\displaystyle \frac{\delta m_\sigma }{m_\sigma }}{\displaystyle \frac{g_s^2}{4\pi }}{\displaystyle \frac{2\kappa }{2\kappa +m_\sigma }}e^{(2\kappa +m_\sigma )b}75{\displaystyle \frac{\delta m_\sigma }{m_\sigma }}`$ (20)
A variation of $`m_\omega `$ gives

$`{\displaystyle \frac{\delta Q_d}{Q_d}}={\displaystyle \frac{m_\omega }{Q_d}}{\displaystyle \frac{\delta m_\omega }{m_\omega }}{\displaystyle \frac{g_v^2}{4\pi }}{\displaystyle \frac{2\kappa }{2\kappa +m_\omega }}e^{(2\kappa +m_\omega )b}80{\displaystyle \frac{\delta m_\omega }{m_\omega }}`$ (21)
The sign difference between these two derivatives and very large derivative value is due to fine tuning between omega and sigma terms, which are separately much larger than the sum.

The next step one can do analytically is to add a simplest square potential to the hard core. The energy of a shallow level in such potential is equal to :

$$E_d=Q_d=\frac{\pi ^2}{16}\frac{(UU_0)^2}{U_0},$$
(22)
$$U_0=\frac{\pi ^2}{8ma^2}.$$
(23)
Here $`U`$ and $`a`$ are the depth and width of the potential well ($`a=cb`$, where $`c`$ and $`b`$ are outer and inner radii), $`m=m_N/2`$ is the reduced mass. Selecting the width and depth of the well to be $`a=1.6fm,U_0=40.2MeV,U=52.6MeV,`$ we got

$`{\displaystyle \frac{\delta Q_d}{Q_d}}81.6e^{m_\sigma b}{\displaystyle \frac{\delta m_\sigma }{m_\sigma }}`$ (24)
$`{\displaystyle \frac{\delta Q_d}{Q_d}}87.4e^{m_\omega b}{\displaystyle \frac{\delta m_\omega }{m_\omega }}`$ (25)
By changing the core radius from $`b=0`$ to $`b=0.4`$ fm one can vary the answer by about factor 3. Simple exponential dependence on the core radius $`b`$ appears because of translational invariance of 1d Scredinger equation for $`r\psi (r)`$.

Another effect one should consider is the modification of the nucleon mass: its contribution to modification of the deuteron binding is

$`\delta Q_d={\displaystyle \frac{\delta M_N}{M_N}}<d|{\displaystyle \frac{p^2}{2M_N}}|d>`$ (26)
which leads to

$$\frac{\delta Q}{Q}=\frac{U+U_0}{UU_0}\frac{\delta m_N}{m_N}7.7\frac{\delta m_N}{m_N}.$$
(27)
Although the sensitivity to the nucleon mass is much weaker than that for mesons, it is still quite strong: we attribute it to the fact that the small deuteron binding energy is in turn the delicate balance between larger kinetic and the potential energies.

### B Using the Walecka potential

The estimates of the preceding section are given for orientation only, and in fact one of course have to solve numerically the radial Schreodinger equation and obtain the correct wave function. The one can either average the potential derivative over it or simply vary all masses involved explicitly. We did the latter and determined the sensitivity of $`Q_d`$ to the sigma, omega and nucleon mass.

Strictly speaking, at this point it is no longer possible to limit ourselves to Walecka model, with the coupling constants (3), since it does not describe correctly the deuteron binding. In fact, by ignoring all spin-dependent forces one cannot even separate the spin-singlet and the spin-triplet states. The tensor forces, attributed to pion and rho exchanges, are needed for this task. Instead of doing so, we have chosen to modify a bit the strength of the omega term, reducing $`g_v^2`$ by a factor 0.953 as compared to (3) and obtaining the correct deuteron binding.

Our results are***In the derivative over the omega mass we have divided back by the factor 0.953, restoring the original strength of the vector term.

$`{\displaystyle \frac{\delta Q_d}{Q_d}}48{\displaystyle \frac{\delta m_\sigma }{m_\sigma }}26{\displaystyle \frac{\delta m_s}{m_s}}`$ (28)
$`{\displaystyle \frac{\delta Q_d}{Q_d}}50{\displaystyle \frac{\delta m_\omega }{m_\omega }}7.5{\displaystyle \frac{\delta m_s}{m_s}}`$ (29)
$`{\displaystyle \frac{\delta Q_d}{Q_d}}6{\displaystyle \frac{\delta m_N}{m_N}}1.1{\displaystyle \frac{\delta m_s}{m_s}}`$ (30)
One can see that the first two derivatives are more sensitive to exact shape of the wave function: they agree qualitatively but not quantitatively with the analytic estimates above. We will not show here such details as exact and approximate wave functions, but just comment that the difference between those explain the difference in the integrals.

Summing all the contributions we find

$`{\displaystyle \frac{\delta Q_d}{Q_d}}17{\displaystyle \frac{\delta m_s}{m_s}}`$ (31)
Using limits on Big Bang Nucleosynthesis from

$`|{\displaystyle \frac{\delta Q_d}{Q_d}}|<0.1`$ (32)
one gets the final limit on the $`m_s`$ variation to be

$`\left|{\displaystyle \frac{\delta (m_s/\mathrm{\Lambda }_{QCD})}{(m_s/\mathrm{\Lambda }_{QCD})}}\right|<0.006`$ (33)
## IV  OKLO

In this section we extract limits on $`\delta m_s`$ following from data on natural nuclear reactor in Oklo active about 2 bn years ago. The most sensitive phenomenon (used previously for limits on the variations of the electromagnetic $`\alpha `$) is disappearance of certain isotopes (especially $`Sm^{149}`$) possessing a neutron resonance close to zero . Today the lowest resonance energy is only $`E_0=0.0973\pm 0.0002eV`$ is large compared to its width, so the neutron capture cross section $`\sigma 1/E_0^2`$. The data constrain the ratio of this cross section to a non-resonance one (which was used to measure number of neutrons emitted by the reactor). It therefore implies Of course, under assumption that the $`same`$ resonance was the lowest one at the time of Oklo reactor. that these data constrain the variation of the following ratio $`\delta (E_0/E_1)`$ where $`E_11MeV`$ is a typical single-particle energy scale, which may be viewed as the energy of some 1-body “doorway” state.

A generic expression for the level energy in terms of fundamental parameters of QCD can be written as follows

$`\delta E_0=A\delta \mathrm{\Lambda }_{QCD}+B_q\delta m_q+B_s\delta m_s+C\delta \alpha \mathrm{\Lambda }_{QCD}`$ (34)
where $`A,B,C`$ are some coefficients. The first term is the basic QCD term, while others are corrections due to modification of the quark masses and the electromagnetic $`\alpha `$.

In this section we provide new estimate of the $`B_s`$. More specifically, we estimate the variation of the resonance energy resulting from a modification of the sigma mass. The energy of the resonance $`E_0=E_{excitation}S_n`$ consists of excitation energy of a compound nucleus, minus the neutron separation energy $`S_n`$. This, in turn, is a depth of the potential well $`V`$ minus the neutron Fermi energy $`ϵ_F`$, $`S_n=Vϵ_F`$. The latter scales like $`1/R^2`$ if the radius of the well is changed. The kinetic part of the excitation energy $`E_{excitation}`$ scales in the same way. So,

$`E_0=E_{excitation}S_n=E_{excitation}+ϵ_FV=K{\displaystyle \frac{\mathrm{}^2}{MR}}V`$ (35)
where K is a numerical constant which can be found from the present time condition $`E_00`$. The shift of the resonance then is

$`\delta E_0=K{\displaystyle \frac{\mathrm{}^2}{MR}}\left({\displaystyle \frac{\delta M}{M}}+{\displaystyle \frac{2\delta R}{R}}\right)\delta V=V\left({\displaystyle \frac{\delta M}{M}}+{\displaystyle \frac{2\delta R}{R}}+{\displaystyle \frac{\delta V}{V}}\right)`$ (36)
Using eq. (2) we can find the depth of the potential well Note that the suppression of N-N wave function at small separation due to the repulsive core reduces the depth of the effective potential $`V`$. However, this effect is not so important in the ratio $`\frac{\delta V}{V}`$.

$`V={\displaystyle \frac{3}{4\pi r_0^3}}\left({\displaystyle \frac{g_s^2}{m_\sigma ^2}}{\displaystyle \frac{g_v^2}{m_\omega ^2}}\right)`$ (37)
Here $`r_0`$= 1.2 fm is an inter-nucleon distance. Numerical estimates shows that the contribution of the variation of $`r_0`$ (and the variation of $`R=A^{1/3}r_0`$ in eq. (36)) is not as important as the direct contribution of the $`m_\sigma `$ variation in the equation above. This gives us

$`{\displaystyle \frac{\delta V}{V}}=8.6{\displaystyle \frac{\delta m_\sigma }{m_\sigma }}+6.6{\displaystyle \frac{\delta m_\omega }{m_\omega }}=3.5{\displaystyle \frac{\delta m_s}{m_s}}`$ (38)
$`\delta E_0=1.7\times 10^8eV\times {\displaystyle \frac{\delta m_s}{m_s}}`$ (39)
We used $`V`$= 50 MeV in eq. (36). Comparison of this result with the the observational limits claimed in  $`|\delta E_0|<0.02`$ eV gives a very strong limit

$`|{\displaystyle \frac{\delta (m_s/\mathrm{\Lambda }_{QCD})}{(m_s\mathrm{\Lambda }_{QCD})}}|<1.2\times 10^{10}`$ (40)
at $`time1.8`$ billion years ago.

Note that the authors of the last work in  found also the non-zero solution $`\delta E_0=0.097\pm 0.008`$ eV. This solution corresponds to the same resonance moved below thermal neutron energy. In this case

$`{\displaystyle \frac{\delta (m_s/\mathrm{\Lambda }_{QCD})}{(m_s\mathrm{\Lambda }_{QCD})}}=(0.56\pm 0.05)\times 10^9`$ (41)
The error here does not include the theoretical uncertainty.

The production of nuclei with $`A>5`$ during BBN is strongly suppressed because of the absence of stable nuclei with $`A=5`$. 5He is unstable nucleus which is seen as a resonance in n-$`\alpha `$ elastic scattering. The ground state lies at 0.89 MeV above neutron threshold. The variation of the strange quark mass may influence the position of the resonance making, for example, 5He stable. Stable 5He at the time of BBN would change strongly the primordial abundances of light elements. The estimate similar to that we made for Sm nucleus gives us a limit

$`{\displaystyle \frac{\delta (m_s/\mathrm{\Lambda }_{QCD})}{(m_s\mathrm{\Lambda }_{QCD})}}>0.006`$ (42)
This limit corresponds to $`\delta E>`$ \- 0.89 MeV at the time of BBN.

We obtained limits on variation of $`m_s/\mathrm{\Lambda }_{QCD}`$ during the interval between the Big Bang and present time and on shorter time scale from Oklo natural nuclear reactor which was active 1.8 billion years ago. It is also possible to obtain limits on the intermediate time scale. One possibility is related to position of the resonance in $`{}_{}{}^{12}C`$ during production of this element in stars. This famous resonance at E=380 KeV is needed to produce enough carbon and create life. According to Ref.  the position of this resonance can not shift by more than 60 KeV ( one can also find in Ref.  the limits on the strong interactions and other relevant references). We have made a very rough estimate of the limit on the strange quark mass variation which can be obtained via $`m_\sigma `$-mechanism:

$`|{\displaystyle \frac{\delta (m_s/\mathrm{\Lambda }_{QCD})}{(m_s\mathrm{\Lambda }_{QCD})}}|<0.001`$ (43)
This limit can be improved after an accurate calculation.

Acknowledgments We thank J.Bjorken, G.Brown and V.Dmitriev for helpful discussions. One of the authors (VF) is thankful to the State University of New York and Institute for Theoretical Atomic and Molecular Physics, Harvard-Smithsonian Center for Astrophysics for hospitality and support. This work is also supported by the Australian Research Council. ES work is partially supported by the US-DOE grant No. DE-FG02-88ER40388.


format: text

text:
# Testing the stability of fundamental constants with the 199Hg+ single-ion optical clock

Present address: \] BNM-SYRTE, Observatoire de Paris, France Permanent address: \] Department of Physics, University of Notre Dame, also supported by DOE and NSF.

## Abstract

Over a two-year duration, we have compared the frequency of the 199Hg+ $`5d^{10}6s{}_{}{}^{2}S_{1/2}^{}(F=0)5d^96s^2{}_{}{}^{2}D_{5/2}^{}(F=2)`$ electric-quadrupole transition at $`282`$ nm with the frequency of the ground-state hyperfine splitting in neutral 133Cs. These measurements show that any fractional time variation of the ratio $`\nu _{\mathrm{Cs}}/\nu _{\mathrm{Hg}}`$ between the two frequencies is smaller than $`\pm 7\times 10^{15}`$yr-1 ($`1\sigma `$ uncertainty). According to recent atomic structure calculations, this sets an upper limit to a possible fractional time variation of $`g_{\mathrm{Cs}}(m_e/m_p)\alpha ^{6.0}`$ at the same level.

Atomic Frequency Standards,Fundamental Constants,Precision Spectroscopy

The development of string theory models aiming at a unified description of gravity and quantum mechanics has renewed interest for improved experimental tests of Einstein’s Equivalence Principle (EEP). Indeed, a common feature of these models is that they allow for, or even predict, violations of EEP Damour et al. (2002). These include violation of the universality of free-fall as well as variation of fundamental constants with time and space. Interestingly, a recent analysis of the spectrum of quasars Webb et al. (2001) suggests that the fine-structure constant $`\alpha `$ may have changed over the cosmological timescale ($`10^{10}`$ yr). The Oklo reactor analysis Shlyakhter (1976); Damour and Dyson (1996) on the other hand puts a stringent limit to possible variation of $`\alpha `$ on the geological timescale ($`10^9`$ yr). Owing to their high accuracy, comparisons between atomic frequency standards based on different atomic species and/or types of transitions provide one of the best ways to perform laboratory tests of the stability of fundamental constants. Present and future efforts to improve atomic frequency standards, in both the optical and the microwave domains, will improve these tests, leading to significant constraints on theoretical work aimed at a unified theory.

In this letter we describe frequency comparisons conducted over a two-year period between a 199Hg+ single-ion optical clock and a 133Cs fountain atomic clock that set a new stringent limit to a possible variation of fundamental constants. The theoretical background for such a test is given first. We then describe the experiment and conclude with the results of the test.

In the experiment, the frequency $`\nu _{\mathrm{Hg}}`$ of the 199Hg+ $`5d^{10}6s{}_{}{}^{2}S_{1/2}^{}(F=0)5d^96s^2{}_{}{}^{2}D_{5/2}^{}(F=2,m_F=0)`$ electric-quadrupole transition at $`\lambda =282`$ nm is compared to the frequency $`\nu _{\mathrm{Cs}}`$ of the ground-state hyperfine transition $`6S_{1/2}(F=3,m_F=0)6S_{1/2}(F=4,m_F=0)`$ in neutral 133Cs. Including relativistic and many-body effects, $`\nu _{\mathrm{Hg}}`$ can be expressed as $`\nu _{\mathrm{Hg}}R_yF_{\mathrm{Hg}}(\alpha )`$, where $`R_y=R_{\mathrm{}}c`$ is the Rydberg constant expressed as a frequency, and $`F_{\mathrm{Hg}}(\alpha )`$ is a dimensionless function of the fine-structure constant $`\alpha =e^2/4\pi \epsilon _0\mathrm{}c`$. Similarly, the hyperfine frequency of cesium can be approximated by $`\nu _{\mathrm{Cs}}g_{\mathrm{Cs}}(m_e/m_p)\alpha ^2R_yF_{\mathrm{Cs}}(\alpha )`$, where $`g_{\mathrm{Cs}}`$ is the gyro-magnetic ratio of the 133Cs nucleus and $`m_e/m_p`$ the electron-to-proton mass ratio. $`F_{\mathrm{Hg}}(\alpha )`$ and $`F_{\mathrm{Cs}}(\alpha )`$ are calculated in Prestage et al. (1995); Dzuba et al. (1999). With numerical values from Dzuba et al. (1999), we find

$`\alpha {\displaystyle \frac{}{\alpha }}\mathrm{ln}F_{\mathrm{Hg}}(\alpha )3.2`$ (1)
$`\alpha {\displaystyle \frac{}{\alpha }}\mathrm{ln}F_{\mathrm{Cs}}(\alpha )+0.8.`$ (2)
Therefore, sequential measurements of the ratio $`\nu _{\mathrm{Cs}}/\nu _{\mathrm{Hg}}`$ actually test the stability of the product of fundamental constants $`g_{\mathrm{Cs}}(m_e/m_p)\alpha ^{6.0}`$. Note that the sensitivity to variation of the individual constants needs to be known to only 1-10 $`\%`$ in order to adequately describe how this comparison constrains possible variations of these fundamental constants. It is therefore justified to omit higher-order terms in the expression of the hyperfine splitting of 133Cs, such as finite nuclear size. The high sensitivity of the Hg+ vs Cs frequency comparison to a change of $`\alpha `$ arises from the large relativistic effects encountered in heavy atoms combined with the negative sign of the relativistic effects in 199Hg+. These factors make the 199Hg+ $`{}_{}{}^{2}S_{1/2}^{}{}_{}{}^{2}D_{5/2}^{}`$ optical transition one of the best choices for a search of a variation of $`\alpha `$ Dzuba et al. (1999).

The 199Hg+ single-ion optical frequency standard has been described previously Young et al. (1999a); Rafac et al. (2000) and only the main features are outlined here. A single 199Hg+ ion is stored in a radio-frequency (rf) spherical Paul trap held at a cryogenic temperature ($`4.2`$ K). It is laser-cooled to near the $`1.7`$ mK Doppler limit using the strongly allowed $`{}_{}{}^{2}S_{1/2}^{}(F=1){}_{}{}^{2}P_{1/2}^{}(F=0)`$ transition at 194 nm, as shown in Fig. 1. After a cooling and state preparation phase ($`16`$ ms), the ion is left in the $`{}_{}{}^{2}S_{1/2}^{}(F=0)`$ lower state of the clock transition. The cooling laser is then switched off and the probe laser light at $`282\mathrm{nm}`$ is directed onto the ion for a typical duration of $`T=50`$ ms. The cooling laser is turned on again to determine the ion’s internal state using the technique of electron shelving Dehmelt (1975); Bergquist et al. (1987). The $`194\mathrm{nm}`$ fluorescence photons emitted by the ion are counted for $`15`$ ms. The absence of scattered photons indicates that the ion has been excited to the $`{}_{}{}^{2}D_{5/2}^{}(F=2)`$ state by the $`282\mathrm{nm}`$ probe laser. Similarly, a typical count rate of $`6000`$ s-1 from the ion indicates that it has remained in the $`{}_{}{}^{2}S_{1/2}^{}(F=0)`$ state after the probe period. When the ion is detected in the $`{}_{}{}^{2}D_{5/2}^{}`$ excited state, the $`194\mathrm{nm}`$ is left on until scattered photons are detected again, indicating that the ion has spontaneously decayed to the $`{}_{}{}^{2}S_{1/2}^{}`$ state. A new interrogation cycle is then started.

The 282 nm radiation used to probe the clock transition is obtained by frequency-doubling 563 nm light from a dye laser in a deuterated ammonium dihydrogen phosphate (ADP) crystal, as shown in Fig. 2. The light from the dye laser is stabilized to a resonance of a stable high-finesse ($`200000`$) Fabry-Pérot cavity. Two acousto-optic modulators (AOM1,AOM2) driven by precision rf sources shift and fine-tune the frequency of the 563 nm radiation so as to match the frequency of the 282 nm radiation to the $`SD`$ resonance frequency. By compensating the predictable linear drift of the stable reference cavity ($`500`$ mHz s-1) using AOM1, we realize an interrogation oscillator with a fractional frequency instability of $`3\times 10^{16}`$ between 1 s and 10 s, corresponding to a linewidth of $`640`$ mHz at 282 nm Young et al. (1999a). For a typical probe time $`T=50`$ ms, which is shorter than the $`90`$ ms spontaneous decay time of the $`{}_{}{}^{2}D_{5/2}^{}`$ state, the observed linewidth is Fourier-transform limited with a full width at half maximum of $`16`$ Hz, much larger than the frequency fluctuations of the probe laser during the probe pulse time. Under these conditions, the fluctuations of the measured transition probability, and hence the frequency stability of the Hg+ frequency standard, are limited mainly by atomic quantum projection noise Itano et al. (1993).

In our trap, the radial and axial secular frequencies are 1.1 MHz and 2.2 MHz respectively. Both frequencies are much higher than the recoil frequency $`\nu =h/2m\lambda ^212.6\mathrm{kHz}`$, where $`m`$ is the ion mass and $`h`$ is the Planck constant. The observed spectrum therefore consists of a central feature at the $`{}_{}{}^{2}S_{1/2}^{}(F=0){}_{}{}^{2}D_{5/2}^{}(F=2,m_F=0)`$ transition resonance frequency, which is free of the recoil frequency shift, together with vibrational sidebands at integer multiples of the trap secular frequencies Bergquist et al. (1987). The spectrum also exhibits features corresponding to the carrier and vibrational sidebands of the $`{}_{}{}^{2}S_{1/2}^{}(F=0){}_{}{}^{2}D_{5/2}^{}(F=2,m_F=\pm 1,\pm 2)`$ “Zeeman components.” Typically, the frequency difference between the clock transition and the first Zeeman components is $`6`$ MHz, which corresponds to a quantization magnetic field $`B0.3`$ mT.

In order to lock the frequency of the probe laser to the atomic resonance, a square-wave frequency modulation is applied to the 563 nm beam using AOM2 so as to probe the carrier of the clock transition alternately on each side of the resonance. Typically, 24 measurements are averaged on each side of the resonance to find the transition probability. The difference between these two transition probabilities reveals the detuning between the center frequency of the probe laser and the center of the atomic resonance. This error signal is used in a digital servo loop that steers the average frequency of the probe laser to the center of the atomic resonance, applying frequency corrections to the synthesizer driving AOM1. The time constant of this servo loop is on the order of $`\tau _{\mathrm{loop}}15`$ s. The analysis of the frequency corrections indicates that for $`\tau >\tau _{\mathrm{loop}}`$, the fractional frequency instability of the probe light stabilized to the atomic resonance is $`57\times 10^{15}\tau ^{1/2}`$, in good agreement with the quantum-limited instability expected for our experimental conditions (trap secular frequencies, ion temperature and measurement cycle time).

In order to compare the Hg+ optical frequency standard to other frequency standards, some fraction of the 563 nm light at the output of AOM1 is launched into a 180 m long optical fiber and delivered into a separate room where it is frequency shifted by AOM3. AOM3 is used to actively cancel the optical path length fluctuations of the fiber in order to preserve the high degree of coherence of the light Young et al. (1999b); Ma et al. (1996).

The frequency of the 563 nm light is measured with an optical frequency comb generated by a mode-locked Ti:Sapphire laser whose femtosecond pulses are spectrally broadened in a microstructure fiber Holzwarth et al. (2000); Diddams et al. (2000); Udem et al. (2001); Diddams et al. (2001). The optical spectrum at the output of the fiber consists of equally spaced, phase-coherent modes with frequencies $`f_n=nf_r+f_0`$, where $`f_r`$ is the repetition rate of the mode-locked laser, $`n`$ is an integer and $`f_0`$ is a frequency offset caused by the difference between the phase and group velocities in the laser cavity. In our setup, the repetition rate is typically $`f_r1`$ GHz and the spectrum spans more than one octave from $`520`$ nm to $`1170`$ nm. The repetition rate is detected directly with a fast photodiode. The offset frequency $`f_0`$ is detected by the self-referencing method, where the frequency-doubled red part of the comb $`2f_n=2nf_r+2f_0`$ is heterodyned with the blue part of the comb $`f_{2n}=2nf_r+f_0`$. We also detect the beatnote $`f_b`$ between the 563 nm light from the Hg+ optical frequency standard and the closest mode of the optical frequency comb $`f_m=mf_r+f_0`$ ($`m`$ is known from previous coarse measurements of the Hg+ $`SD`$ transition frequency). In practice, $`f_0`$ and $`f_b`$ are phase-locked to precision rf sources by acting on the cavity length and the pump power of the femtosecond laser, respectively. Finally, $`f_r`$ is measured by counting the low-frequency beatnote between $`f_r`$ and a third rf synthesizer.

In order to perform absolute optical frequency measurements in terms of the SI second, all synthesizers and frequency counters are referenced to the $`5`$ MHz output of a hydrogen maser with a typical frequency instability $`2\times 10^{13}`$ at 1 s. The maser itself is part of an ensemble of 5 masers and 3 commercial cesium clocks used to realize the local NIST atomic time scale, which is in turn periodically calibrated using the NIST-F1 cesium fountain primary standard Jefferts et al. (2002), as well as international cesium standards. The frequency of the reference hydrogen maser is known within $`1.8`$ parts in $`10^{15}`$ with respect to the ground-state hyperfine splitting of 133Cs. As shown in recent investigations Udem et al. (1999); Diddams et al. (2002), the additional noise and inaccuracy from the optical frequency comb itself is negligible at this level. Using this setup, it is thus possible to perform absolute measurements of the frequency of the Hg+ optical standard (together with the frequency of each component of the optical comb) with a fractional frequency uncertainty of $`1.8\times 10^{15}`$. Typically, measurements are performed for 2 hours, leading to a statistical (type A) fractional frequency uncertainty of 2.4 parts in $`10^{15}`$, which corresponds to a $`2.5`$ Hz uncertainty on the frequency of the 282 nm stabilized light.

The full evaluation of all systematic effects of the Hg+ optical standard is still under way. At the present time, $`10`$ Hz is a conservative upper bound for the total systematic (type B) uncertainty, as shown by the following preliminary analysis of systematic shifts. The second-order Zeeman frequency shift of the clock transition is given by $`\delta \nu _Z=K_Z^{(2)}B^2,K_Z^{(2)}=18.925(28)`$ kHz mT-2 Itano (2000). The typical field $`B0.3`$ mT corresponds to a second-order Zeeman frequency shift $`\delta \nu _Z1.7`$ kHz. Therefore, the $`1.5\times 10^3`$ fractional uncertainty on coefficient $`K_Z^{(2)}`$ leads to a $`2.6`$ Hz uncertainty on the clock frequency. For the same value of the bias magnetic field, the sensitivity of the clock frequency to field fluctuations is $`11`$ Hz $`\mu `$T-1. In our unshielded environment, $`B`$-field fluctuations up to $`\pm 0.2\mu \mathrm{T}`$ have been observed, leading to a $`2.2`$ Hz uncertainty. The most troublesome frequency shift is the electric-quadrupole shift $`\delta \nu _Q`$ due to the coupling between the atomic electric-quadrupole moment in the upper $`{}_{}{}^{2}D_{5/2}^{}(F=2,m_F=0)`$ clock state with electric field gradients due, for example, to stray charges on the trap electrodes. From our trap geometry and the bias voltages applied to compensate for stray electric fields, it is possible to calculate an upper bound for an electric field gradient at the location of the ion and to estimate that $`\delta \nu _Q<1`$ Hz Itano (2000). In addition, two successive versions of the trap electrodes (gold and molybdenum surfaces) have been used, leading to no detectable change of the clock frequency. Similarly, the background pressure of helium (the only species remaining with significant pressure at 4.2 K) has been changed by more than one order of magnitude (estimated by measuring the ion heating rate due to background collisions) without producing any detectable frequency shift. Empirically, this implies that the He background pressure shift is smaller than $`1`$ Hz. At a temperature of $`300`$ K, the black-body radiation shift is $`0.08`$ Hz, and it is considerably lower in the $`4.2`$ K cryogenic environment Itano (2000). At the Doppler cooling limit of $`1.7`$ mK, the second-order Doppler shift due to thermal motion is $`0.003`$ Hz. Finally, by suitably compensating the stray static electric field, systematic shifts related to the trapping oscillating electric field (second-order Doppler shift due to micromotion, AC Stark shift) are made smaller than $`0.1`$ Hz Berkeland et al. (1998).

Over a period of two years, we have performed 20 measurements of the 199Hg+ $`SD`$ transition frequency with respect to the 133Cs ground state hyperfine splitting. Figure 3 shows these measurements, corrected for the second-order Zeeman frequency shift, with their statistical (type A) error bars. The absolute frequency of the Hg+ optical standard is given by the weighted average of these data: $`\nu _{\mathrm{Hg}}=1064721609899143.4`$ Hz. The total statistical uncertainty is only $`1.0`$ Hz. Measurements in Fig. 3 clearly show a reproducibility better than $`10`$ Hz at $`1.06\times 10^{15}`$ Hz, the most stringent comparison of optical and microwave frequencies to date. The total uncertainty of the measurement is dominated by the $`10`$ Hz systematic uncertainty of the Hg+ clock. With a $`10^{14}`$ fractional uncertainty and a measurement period of two years, our measurement constrains a possible fractional variation of $`\nu _{\mathrm{Cs}}/\nu _{\mathrm{Hg}}`$ at the level of $`\pm 7\times 10^{15}`$ yr-1 ($`1\sigma `$ uncertainty).

This result can be interpreted as constraining a possible fractional variation of $`g_{\mathrm{Cs}}(m_e/m_p)\alpha ^{6.0}`$ at the same level: $`\pm 7\times 10^{15}`$ yr-1. Assuming that any change in this quantity is due to the $`\alpha ^{6.0}`$ factor, we derive an upper bound for a possible linear variation of the fine-structure constant: $`|\dot{\alpha }/\alpha |<1.2\times 10^{15}`$ yr-1, a factor of $`30`$ improvement over Prestage et al. (1995). However, there may be a significant change of $`g_{\mathrm{Cs}}`$ or $`m_e/m_p`$ due to possible variation of the strength of the strong and electroweak interactions. In fact, recent theoretical work argues that within the framework of a grand unified theory a fractional variation of $`\alpha `$ is necessarily accompanied by a fractional variation of $`m_e/m_p`$ that is $`38`$ times larger Calmet and Fritzsch (2002). Fortunately, comparisons between two optical clocks will test the stability of $`\alpha `$ alone Dzuba et al. (1999). One interesting possibility is to compare the 199Hg+ $`SD`$ transition to the 657 nm $`{}_{}{}^{1}S_{0}^{}(m=0){}_{}{}^{3}P_{1}^{}(m=0)`$ transition in neutral 40Ca Udem et al. (2001); Helmcke et al. (2002). With an independent constraint to the stability of $`\alpha `$, comparisons involving hyperfine transitions, such as the present $`\nu _{\mathrm{Hg}}`$ (optical) vs $`\nu _{\mathrm{Cs}}`$ comparison or the $`\nu _{\mathrm{Rb}}`$ (microwave) vs $`\nu _{\mathrm{Cs}}`$ comparison Bize et al. (2001); H. Marion et al. (2001), will test the stability of the strong and electroweak interactions. Reference Karshenboïm (2000) investigates in more detail possible laboratory tests of the stability of fundamental constants.

We thank Robert Windeler of OFS laboratories for providing the microstructure fiber and Thomas Udem for his contributions in the early measurements. This work was partially supported by the Office of Naval Research and through a cooperative research and development agreement with Timing Solutions, Inc., Boulder, CO. This work of an agency of the U.S. government is not subject to U.S. copyright.


text:
# Kohn-Luttinger instability of the 𝑡-𝑡' Hubbard model in two dimensions: variational approach

## I Introduction

In a seminal paper, Kohn and Luttinger Kohn65  (KL) presented a perturbative argument showing that the generic ground state of a Fermi liquid at weak coupling is superconducting, even if the bare electron-electron interaction is not attractive in the Cooper channel in any angular momentum sector. Recently, this result has been rederived utilizing various versions of the renormalization group method.Shankar94 ; Zanchi97 ; Honerkamp01

However, no attempt has been made so far to construct a variational wavefunction for a KL superconductor, which could be (at least in principle) tested in a direct numerical simulation. The purpose of this paper is to construct such a wavefunction for the repulsive Hubbard model in the limit of weak coupling $`UW`$, where $`U`$ is the local repulsive interaction and $`W`$ is the electron bandwidth. Our strategy follows quite closely the canonical transformation approach to the electron-phonon problem introduced by Fröhlich Frohlich52  and sketched in the Appendix. In some sense it is also similar to the customary analysis of the Hubbard model close to half filling in the limit of strong coupling $`UW`$. Namely, in the latter case a canonical transformation can be found Harris67 ; Gros87  which transforms the purely repulsive Hubbard model to the so-called $`t`$-$`J`$ model which explicitly contains an attractive interaction (of the order $`JW^2/U`$) favoring the formation of singlet superconductivity. Unfortunately, the $`t`$-$`J`$ model contains a constraint on the local number of electrons and as such is very difficult to study. Here, on the other hand, we construct a canonical transformation which eliminates the scattering of quasiparticles to first order in $`U/W`$ and generates a weak-coupling model with an attractive interaction of the order $`U^2/W`$.

The outline of the paper is as follows. The construction of the effective model is described in Section II and our method is applied to the two-dimensional $`t`$-$`t^{}`$ Hubbard model in Section III.

## II Effective Hamiltonian

We consider the Hubbard model on a square lattice with periodic boundary conditions and $`L=l\times l`$ sites, described by the Hamiltonian $`H=H_0+H_1+H_2`$, where $`H_0=_{𝐤,\sigma }\epsilon _𝐤n_{𝐤,\sigma }`$ is the kinetic energy operator and the interaction term has been split into two parts:

$`H_1`$ $`=`$ $`{\displaystyle \frac{U}{L}}\left({\displaystyle \frac{N^2}{4}}𝐒^2\right)+{\displaystyle \frac{U}{L}}{\displaystyle \underset{𝐤,𝐩}{}}c_𝐤^{}c_𝐤^{}c_𝐩c_𝐩,`$ (1)
$`H_2`$ $`=`$ $`{\displaystyle \frac{U}{L}}{\displaystyle \underset{\{𝐤\}}{\overset{}{}}}c_{𝐤_3}^{}c_{𝐤_1}c_{𝐤_4}^{}c_{𝐤_2}\delta _{𝐤_1+𝐤_2,𝐤_3+𝐤_4},`$ (2)
where $`N=N_{}+N_{}`$ and $`𝐒`$ are the total electron number and the total spin, respectively. Terms which are not extensive in the thermodynamic limit $`L\mathrm{}`$ have been neglected. $`H_2`$ is the generic interaction term which scatters electrons from $`𝐤_1`$ to $`𝐤_3`$ and from $`𝐤_2`$ to $`𝐤_4`$. The prime on the summation means that terms with $`𝐤_1=𝐤_3`$ (forward scattering channel), $`𝐤_1=𝐤_4`$ (exchange channel), and $`𝐤_1+𝐤_2=0`$ (Cooper channel) are excluded from $`H_2`$ and are singled out into $`H_1`$. The operator $`H_0+H_1`$ can be thought of as a reduced BCS Hamiltonian of a Landau Fermi liquid, while $`H_2`$ contains the scattering processes leading to a finite lifetime of the Landau quasiparticles.

Now we look for a canonical transformation from the bare electrons to the dressed ones, $`\stackrel{~}{H}=e^{iS}He^{iS}`$, such that the scattering of the quasiparticles vanishes to first order in $`U`$. This happens if $`H_2+i[S,H_0]=0`$ in which case we can write $`\stackrel{~}{H}=H_0+H_1+i[S,H_1]+i[S,H_2]/2`$, where terms of order $`O(U^3)`$ have been neglected. Making use of the identity

$$[c_P^{}c_Q,c_R^{}c_S]=\delta _{QR}c_P^{}c_S\delta _{PS}c_R^{}c_Q$$
(3)
where the indices $`P,Q,R,S`$ include also the spin label, one verifies readily that

$$S=\frac{iU}{L}\underset{\{𝐤\}}{\overset{}{}}\frac{\delta _{𝐤_1+𝐤_2,𝐤_3+𝐤_4}}{\epsilon _{𝐤_1}+\epsilon _{𝐤_2}\epsilon _{𝐤_3}\epsilon _{𝐤_4}}c_{𝐤_3}^{}c_{𝐤_1}c_{𝐤_4}^{}c_{𝐤_2}$$
(4)
is the Hermitian generator we are looking for. Note that the forward, exchange, and Cooper channel processes do not contribute to $`S`$. In order to calculate the commutators $`[S,H_1]`$ and $`[S,H_2]`$ we first note that the operators $`A=c_3^{}c_1`$ and $`B=c_\gamma ^{}c_\alpha `$ commute with the operators $`X=c_4^{}c_2`$ and $`Y=c_\delta ^{}c_\beta `$ and therefore

$$[AX,BY]=[A,B]XY+BA[X,Y].$$
(5)
Combining Eqs. (3,5) we thus find the identity

$`[c_3^{}c_1c_4^{}c_2,c_\gamma ^{}c_\alpha c_\delta ^{}c_\beta ]=(\delta _{1\gamma }c_3^{}c_\alpha \delta _{3\alpha }c_\gamma ^{}c_1)c_4^{}c_2c_\delta ^{}c_\beta +(\delta _{2\delta }c_4^{}c_\beta \delta _{4\beta }c_\delta ^{}c_2)c_\gamma ^{}c_\alpha c_3^{}c_1,`$
making use of which one can straightforwardly calculate the effective Hamiltonian $`\stackrel{~}{H}`$. The result is fairly involved, including also three-body forces among the electrons. Nevertheless, the expectation value $`E=\stackrel{~}{\psi }|\stackrel{~}{H}|\stackrel{~}{\psi }`$ in a BCS state $`|\stackrel{~}{\psi }=\mathrm{\Pi }_𝐤(u_𝐤+v_𝐤c_𝐤^{}c_𝐤^{})|0`$ takes a simple and physically transparent form,

$$E=\underset{𝐤\sigma }{}\epsilon _𝐤f_{𝐤\sigma }+\frac{U}{L}\frac{N^2}{4}+\frac{1}{L}\underset{𝐤,𝐩}{}V_{\mathrm{𝐤𝐩}}b_𝐤^{}b_𝐩+\frac{U^2}{L^2}\underset{\{𝐤\}}{\overset{}{}}\frac{f_{𝐤_1}f_{𝐤_2}(1f_{𝐤_3})(1f_{𝐤_4})}{\epsilon _{𝐤_1}+\epsilon _{𝐤_2}\epsilon _{𝐤_3}\epsilon _{𝐤_4}}\delta _{𝐤_1+𝐤_2,𝐤_3+𝐤_4},$$
(6)
where we have introduced a momentum-resolved BCS order parameter $`b_𝐩=\stackrel{~}{\psi }|c_𝐩c_𝐩|\stackrel{~}{\psi }`$ and we assumed $`\stackrel{~}{\psi }|𝐒^2|\stackrel{~}{\psi }=0`$. In the Cooper channel the effective interaction reads $`V_{\mathrm{𝐤𝐩}}=U+U^2\chi _1(𝐤+𝐩,\epsilon _𝐩\epsilon _𝐤)`$, where $`\chi _1(𝐪,\omega )`$ is the real part of the particle-hole susceptibility $`\chi (𝐪,\omega )=L^1_𝐊(f_𝐊f_{𝐊+𝐪})/(\epsilon _{𝐊+𝐪}\epsilon _𝐊\omega i0)`$. The physical content of the KL argument is that the effective interaction between the dressed electrons can be attractive in the Cooper channel in some angular momentum sector. If this is the case, then the energy $`E=\stackrel{~}{\psi }|e^{iS}He^{iS}|\stackrel{~}{\psi }`$ is minimized by nonvanishing $`b_𝐩`$ and the wavefunction $`|\psi =e^{iS}|\stackrel{~}{\psi }`$ provides a variational ansatz for the KL superconductor.

Before proceeding we would like to point out that our variational approach bears some similarity to a very recent formulation of the flow equation method.Grote02  The comparison of both methods is left for future.

## III BCS variational solution of the effective Hamiltonian

We wish to apply the above formalism to the study of the $`t`$-$`t^{}`$ Hubbard model in two dimensions whose single particle dispersion is $`\epsilon _𝐤=2t(\mathrm{cos}k_x+\mathrm{cos}k_y)+4t^{}\mathrm{cos}k_x\mathrm{cos}k_y`$, where we have set the lattice constant $`a=1`$. The superconducting phase diagram of the weak-coupling $`t`$-$`t^{}`$ Hubbard model has been studied previously by considering the effective KL interaction $`V_{\mathrm{𝐤𝐩}}^{KL}=U+U^2\chi _1(𝐤+𝐩,0)`$ with both $`𝐤`$ and $`𝐩`$ lying at the Fermi surface.Hlubina99  Roughly speaking, such an approach provides us with the coupling constant $`g`$, but not with the prefactor $`\mathrm{\Omega }`$ in the BCS formula for the transition temperature, $`T_c=\mathrm{\Omega }\mathrm{exp}(1/g)`$.

In what follows we perform the standard BCS minimization of an effective free energy whose $`T=0`$ limit reduces to the variational energy Eq. (6). We assume that the interaction $`U`$ is sufficiently weak so that no particle-hole instabilities can develop. (In particular, this means that we can’t study systems close to the van Hove filling.) Concentrating on lowest order effects in $`U/W`$, we neglect the last term in Eq. (6) which represents the standard second-order perturbation theory effects, and treat it as an additive constant. A more detailed study of Eq. (6) taking into account also the last term is postponed to future publications. Moreover, since we expect that $`T_ct`$, we make use of the $`T=0`$ Fermi functions in the calculation of the particle-hole susceptibility $`\chi (𝐪,\omega )`$.

A standard BCS calculation leads to the self-consistent equation for the gap function $`\mathrm{\Delta }_𝐤=L^1_𝐩V_{\mathrm{𝐤𝐩}}b_𝐩`$,

$$\mathrm{\Delta }_𝐤=\frac{1}{L}\underset{𝐩}{}V_{\mathrm{𝐤𝐩}}\mathrm{\Delta }_𝐩\frac{\mathrm{tanh}(E_𝐩/2T)}{2E_𝐩},$$
(7)
where $`E_𝐩=(\xi _𝐩^2+\mathrm{\Delta }_𝐩^2)^{1/2}`$, $`\xi _𝐩=\epsilon _𝐩\mu `$, and $`\mu `$ is the chemical potential.

For temperatures infinitesimally below the mean field transition temperature, the gap equation can be linearized with respect to $`\mathrm{\Delta }_𝐩`$ and can be written as $`D_𝐤=_𝐩I_{\mathrm{𝐤𝐩}}(T)D_𝐩`$, where $`D_𝐤=\mathrm{\Delta }_𝐤\varphi _𝐤T^{1/2}`$ is a dimensionless gap function, $`\varphi _𝐤=[\mathrm{tanh}(\xi _𝐤/2T)/2\xi _𝐤]^{1/2}`$, and $`I_{\mathrm{𝐤𝐩}}(T)=L^1V_{\mathrm{𝐤𝐩}}\varphi _𝐤\varphi _𝐩`$ is a real symmetric matrix corresponding to a dimensionless pair scattering function. At high temperatures $`I_{\mathrm{𝐤𝐩}}(T)T^1`$ and the gap equation does not have solutions. With decreasing $`T`$ the maximal eigenvalue $`\lambda `$ of $`I_{\mathrm{𝐤𝐩}}(T)`$ grows and the mean field transition temperature $`T_c`$ can be calculated from the condition $`\lambda (T_c)=1`$.

We have studied the eigenvalues of $`I_{\mathrm{𝐤𝐩}}(T)`$ numerically. In order to minimize the cost of the calculation, the matrix $`I_{\mathrm{𝐤𝐩}}`$ was calculated on a special lattice, see Fig. 1. The lattice consists of a sequence $`i=1,2,\mathrm{},n`$ of $`n`$ regular $`l_i\times l_i`$ grids, which become progressively more dense as the distance to the Fermi level diminishes, $`l_{i+1}=2l_i`$. The borders between the subsequent levels are constructed as follows: Let the maximum and minimum band energies be $`\xi _{\mathrm{max}}`$ and $`\xi _{\mathrm{min}}`$ and let $`\xi _0=\mathrm{max}(\xi _{\mathrm{max}},|\xi _{\mathrm{min}}|)`$. We construct a sequence of energies $`\xi _i=q^i\xi _0`$ with a quotient $`q<1`$. For $`i<n`$, the $`i`$-th level grid $`l_i\times l_i`$ is realized for all $`𝐤`$ points satisfying $`\xi _i<|\xi _𝐤|<\xi _{i1}`$, whereas the $`n`$-th grid applies for $`|\xi _𝐤|<\xi _{n1}`$. The lattice is characterized by three parameters: $`q`$, $`l_1`$, and the number of levels $`n`$. In what follows we always take $`q=4^1`$ and $`l_1=32`$, and instead of the number of levels we describe the lattice by specifying the finest $`𝐤`$-point grid $`l_n\times l_n`$.

The susceptibility $`\chi (𝐪,\omega )`$ was calculated on usual lattices $`l_n\times l_n`$ by a straightforward modification of the method described in Ref. Zlatic00,  which makes use of the fast Fourier transform algorithm. We have used $`4l_n`$ time points and the Nyquist frequency was chosen to be $`8W`$.

The matrix $`I_{\mathrm{𝐤𝐩}}`$ was diagonalized by the modified Lanczos method.Dagotto94  The initial vector was generated randomly in one octant (quadrant) of the Brillouin zone and subsequently continued to the whole zone so as to transform according to the one-dimensional $`s`$, $`d`$, $`d_{xy}`$ and $`g`$ (two-dimensional $`p`$) irreducible representations of the point group of the square.notation

As an illustrative example, the temperature dependence of the maximal eigenvalue in the $`d_{xy}`$ sector for $`t^{}/t=0.78`$ and electron density $`\rho =1`$ is shown in Fig. 2 for various $`l_n\times l_n`$. Several points are worth mentioning. For a finite lattice size $`l_n`$, the divergence of the particle-particle susceptibility is cut off at temperatures $`T^{}v_F/l_n`$, and therefore $`\lambda (T)`$ saturates at low temperatures $`TT^{}`$. Thus, even making use of a moderate interaction strengthinteraction  $`U=W/2`$ (this value of $`U`$ is used in most numerical examples) the largest eigenvalue $`\lambda `$ of $`I_{\mathrm{𝐤𝐩}}`$ is typically less than unity down to the lowest temperatures. Therefore we can’t determine $`T_c`$ directly from $`\lambda (T_c)=1`$. Instead, the $`T`$-dependence of the eigenvalue is fitted in the low temperature limit by a BCS-like expression $`\lambda (T)=g\mathrm{log}(\mathrm{\Omega }/T)`$ and afterwards $`T_c`$ is calculated from $`T_c=\mathrm{\Omega }\mathrm{exp}(1/g)`$. This procedure is shown explicitly in Fig. 2. It should be pointed out that usually we can determine the transition temperature with a reasonable degree of confidence only if $`T_c/t>10^6`$. In this sense the parameter values used in Fig. 2 are somewhat special.

Let us also stress that the large value of $`\mathrm{\Omega }`$ estimated from Fig. 2 should not be interpreted as a large pairing energy scale. This is evident from an equivalent but physically more meaningful expression $`\lambda (T)=\lambda (T_0)+g\mathrm{log}(T_0/T)`$ which makes explicit use of the pairing energy $`T_0`$ (its value is arbitrary, but can be estimated from the maximal temperature where the logarithmic scaling of $`\lambda `$ applies, $`T_00.005t`$ in our example). The large value of $`\mathrm{\Omega }=T_0\mathrm{exp}[\lambda (T_0)/g]`$ is a consequence of the large $`\lambda (T_0)`$, which in turn measures the contribution of the high energy processes to $`\lambda `$.

### III.1 Half-filled $`t`$-$`t^{}`$ Hubbard model

Now let us apply the abovementioned method to the study of the superconducting phase diagram of the $`t`$-$`t^{}`$ Hubbard model at half filling ($`\rho =1`$). In the strong-coupling limit $`UW`$, this would correspond to the Heisenberg model $`H=J_1_{𝐢,𝐣}𝐒_𝐢𝐒_𝐣+J_2_{𝐢,𝐣}𝐒_𝐢𝐒_𝐣,`$ where $`𝐢,𝐣`$ and $`𝐢,𝐣`$ are pairs of nearest neighbor and next-nearest neighbor sites, respectively, $`J_1=4t^2/U`$, and $`J_2=4(t^{})^2/U`$. There is general consensus that, due to frustration, for $`0.4<J_2/J_1<0.6`$ the Heisenberg model does not exhibit the Néel order (see, e.g. Ref. Capriotti01,  and references therein). The nature of the emergent state has been discussed intensively, but very recently it has been shownCapriotti01  that a projected BCS wavefunction with a mixed $`d`$ and $`d_{xy}`$ symmetry provides a very good approximation for the $`J_1`$-$`J_2`$ model with $`J_2/J_1`$ close to 1/2.

In this paper we look for a weak-coupling analogue of the phase proposed by Capriotti et al.Capriotti01  To this end, in Fig. 3 we plot $`T_c`$ in the $`d`$-wave and $`d_{xy}`$-wave sectors as a function of $`\alpha =2t^{}/(t+2t^{})`$ for $`U=W/2`$ and $`\rho =1`$. The dominant pairing instability is seen to change from the $`d`$-wave at small $`t^{}/t`$ to the $`d_{xy}`$ symmetry at large $`t^{}/t`$, since the susceptibility $`\chi (𝐪,0)`$ peaks at $`(\pi ,\pi )`$ and $`(\pi ,0)`$ in these two limits. Surprisingly, in the region corresponding to $`J_2/J_11/2`$, namely $`t^{}/t2^{1/2}`$ or $`\alpha 0.59`$, we find that neither $`d`$ nor $`d_{xy}`$-wave pairing instabilities are dominant, and a small island of $`s`$-wave pairing is realized. For $`\alpha [0.45,0.57]`$ the transition temperatures are too low to be accessible by the present method. Therefore we show in Fig. 3 the coupling constant $`g`$ in the (dominant) $`d`$ and $`s`$-wave symmetry sectors calculated according to Ref. Hlubina99, . From the qualitative agreement of these two different methods we conclude that in the vicinity of $`\alpha 0.55`$ the pairing symmetry changes from $`d`$ to $`s`$-wave.

We expect that for $`\alpha 0.59`$, upon increase of the interaction strength there is presumably a true Mott-Hubbard superconductor-insulator transition at some critical value of $`U`$. We hypothesize that before this happens, with increasing $`U`$ the $`d`$ and $`d_{xy}`$ phases grow at the expense of the $`s`$ phase, until at some stage the $`s`$ phase disappears from the phase diagram. On the other hand, underneath the crossing point of the $`d`$-wave and $`d_{xy}`$-wave $`T_c`$’s, we expect a region with mixed $`d+id_{xy}`$ pairing. Such mixed-symmetry pairing can be justified by a simple model calculation for a two dimensional isotropic system with a separable interaction $`V(\varphi ,\varphi ^{})=V_1(\varphi ,\varphi ^{})+V_2(\varphi ,\varphi ^{})`$ in the Cooper channel. Here $`V_1(\varphi ,\varphi ^{})\mathrm{cos}2\varphi \mathrm{cos}2\varphi ^{}`$ with coupling constant $`g_1`$ and cutoff $`\omega _1`$ leads to $`d`$-wave pairing and $`V_2(\varphi ,\varphi ^{})\mathrm{sin}2\varphi \mathrm{sin}2\varphi ^{}`$ with coupling constant $`g_2`$ and cutoff $`\omega _2`$ yields $`d_{xy}`$-wave pairing. For simplicity we further assume that by changing a control parameter ($`t^{}/t`$ in our case), only the coupling constant $`g_2`$ changes. Then, if we denote that value of $`g_2`$ for which the $`T_c`$ crossing happens as $`g_{20}`$, it can be shown readily that a mixed $`d+id_{xy}`$ pairing state is in fact stabilized at zero temperature in a finite window of the control parameter for which $`g_{20}^11/2<g_2^1<g_{20}^1+1/2`$. Note that in the weak-coupling limit this window is quite narrow.

Before proceeding we should mention that the concept of mixed-symmetry superconductivity has been introduced previously in the context of strong coupling models, for the $`t`$-$`J`$ model by KotliarKotliar88  and for the $`t`$-$`t^{}`$-$`J_1`$-$`J_2`$ model in a recent paper by Sachdev.Sachdev01  Our discussion shows that (as pointed out already in Ref.Hlubina99, ), mixed-symmetry superconductivity is in fact a generic consequence of phase diagrams like the one shown in Fig. 3, where superconducting states of different symmetry cross.

### III.2 $`d`$-wave region of the $`t`$-$`t^{}`$ Hubbard model

As another application of the present formalism, in the rest of this paper we concentrate on the $`d`$-wave region of the superconducting phase diagram of the $`t`$-$`t^{}`$ Hubbard modelHlubina99  which is relevant to the cuprates. Close to the van Hove density, superconductivity competes with the SDW stateHlubina97 ; Hlubina99 ; Honerkamp01  and the question whether one of these two states or even more exotic statesFurukawa98  wins is still open. Here we focus on a different question, namely that of the full $`𝐤`$-dependence of the superconducting gap $`\mathrm{\Delta }_𝐤`$ away from the van Hove density. Fig. 4 shows the self-consistent solution $`\mathrm{\Delta }_𝐤`$ of Eq. (7) for $`t^{}/t=0.3`$, $`\rho =0.8`$, and $`T=0`$. The maximal gap on the Fermi surface $`\mathrm{\Delta }_{\mathrm{max}}=0.0023t`$ and the critical temperature $`T_c=0.0013t`$ imply $`\mathrm{\Delta }_{\mathrm{max}}/T_c1.8`$.

The overall shape of $`\mathrm{\Delta }_𝐤`$ is well described by the simplest formula consistent with $`d`$-wave symmetry, $`\mathrm{\Delta }_𝐤\mathrm{cos}k_x\mathrm{cos}k_y`$. Nevertheless, fine structure on top of this overall shape is clearly visible. This is shown more quantitatively in Fig. 5, where $`\mathrm{\Delta }_𝐤`$ is plotted along two $`k_y=`$const cuts of the Brillouin zone. A nontrivial structure is seen to develop in the direction perpendicular to the Fermi surface which is analogous to the strong-coupling effects in conventional superconductors.

It is well known that in the case of phonon mediated superconductivity the energy dependence of $`\mathrm{\Delta }_𝐤`$ leads to a structure in the density of states. In the Appendix we have shown that our method is capable to describe such effects. Therefore we have looked for features in the density of states caused by the structure in $`\mathrm{\Delta }_𝐤`$ in our purely electronic model as well. The electron spectral function $`A(𝐤,\omega )`$ and the density of states $`N(\omega )`$ have been calculated from

$`A(𝐤,\omega )`$ $`=`$ $`u_𝐤^2\delta (\omega E_𝐤)+v_𝐤^2\delta (\omega +E_𝐤),`$ (8)
$`N(\omega )`$ $`=`$ $`L^1{\displaystyle \underset{𝐤}{}}A(𝐤,\omega ),`$ (9)
where $`u_𝐤^2,v_𝐤^2=2^1(1\pm \xi _𝐤/E_𝐤)`$. Fig. 6 shows that at weak coupling (e.g. for $`U=W/2`$ when $`\mathrm{\Delta },T_ct`$) no structures associated with the energy dependence of $`\mathrm{\Delta }_𝐤`$ can be observed. This is because features in the energy dependence of $`\mathrm{\Delta }_𝐤`$ can be visible also in the density of states only if they appear in a $`𝐤`$-point for which the single-particle energy $`\xi _𝐤`$ is not much larger than $`\mathrm{\Delta }_𝐤`$. In other words, in order to be observable in $`N(\omega )`$, the structure has to appear at a distance of at most $`\delta k\mathrm{\Delta }/v_F`$ from the Fermi line. At weak coupling this criterion is not satisfied and therefore Fig. 6 does not exhibit any deviations from the textbook form of the density of states for a $`d`$-wave superconductor.

We have solved the self-consistent equation Eq. (7) also for a larger interaction strength $`U=0.75W`$, when the gap increases by more than an order of magnitude with respect to $`U=W/2`$. In this case we did find a nontrivial structure in the density of states (see Fig. 7), which could however be attributed to the van Hove singularity in the density of states of the noninteracting spectrum rather than to a many body effect.

Now let us turn to the angular dependence of $`\mathrm{\Delta }_𝐤`$. In Fig. 8 we plot $`\mathrm{\Delta }(\phi )/\mathrm{\Delta }(0)`$ along the Fermi line for $`T=0`$ and $`T/t=1`$. The variable $`\phi `$ is the angle between the radius vector of the Fermi surface point under study \[measured with respect to $`M=(\pi ,\pi )`$\] and the $`x`$ axis. The overall shape of the high temperature data is close to the simple $`\mathrm{cos}k_x\mathrm{cos}k_y`$ form of the gap. With decreasing temperature $`\mathrm{\Delta }_𝐤`$ is seen to be locally more and more enhanced close to the crossing point of the Fermi line with the magnetic zone (hot spot). Nevertheless, this so-called hot spot effect is quite weak down to $`T=0`$. Fig. 8 shows explicitly that the hot spot effect is in qualitative agreement with an earlier formulation Hlubina99  in which the energy dependence of the Cooper channel interaction is neglected and $`\mathrm{\Delta }_𝐤`$ lives only on the Fermi line. On the other hand, a recent paperGuinea02  reports (within the formalism of Ref. Hlubina99, ) a much stronger hot spot effect than found here.

## IV Conclusions

In conclusion, we have constructed a canonical transformation of the Hubbard model which eliminates quasiparticle scattering processes to first order in the interaction $`U`$. The resulting effective Hamiltonian contains terms of order $`U^2/W`$ where $`W`$ is the bandwidth, which are attractive in the Cooper channel. This allows us to construct a variational wavefunction for the Kohn-Luttinger superconductor and to show explicitly that even a purely repulsive model gains energy by developing superconducting correlations.

As an application of the method, the effective Hamiltonian has been solved in the BCS approximation and the superconducting phase diagram of the $`t`$-$`t^{}`$ Hubbard model at half filling has been found. The superconducting state at $`t^{}/t0.7`$ is a weak-coupling analog of the wavefunction proposed for the $`J_1`$-$`J_2`$ model by Capriotti et al.Capriotti01

As another application, in the $`d`$-wave region of the $`t`$-$`t^{}`$ Hubbard model Hlubina99  we have found a nontrivial energy dependence of the gap function reminiscent of strong-coupling effects in conventional superconductors. There is a hot spot effect in the angular dependence of $`\mathrm{\Delta }_𝐤`$, which is quite consistent with the estimates from calculationsHlubina99  in which the gap lives only on the Fermi line.

###### Acknowledgements.

This work was supported by the Slovak Scientific Grant Agency under Grant No. VEGA-1/9177/02 and by the Slovak Science and Technology Assistance Agency under Grant No. APVT-51-021602. *

## Appendix A

As a test of the applicability of the variational method to an effective electronic Hamiltonian, we apply it to the case of phonon mediated superconductivity. We consider the simplest model of electrons with a local coupling to a phonon mode. The Hamiltonian of the system can be written $`H=H_0+H^{}`$, where

$`H_0`$ $`=`$ $`{\displaystyle \underset{𝐤,\sigma }{}}\epsilon _𝐤c_{𝐤,\sigma }^{}c_{𝐤,\sigma }+{\displaystyle \underset{𝐪0}{}}\omega _𝐪a_𝐪^{}a_𝐪,`$
$`H^{}`$ $`=`$ $`{\displaystyle \frac{D}{\sqrt{L}}}{\displaystyle \underset{𝐤,\sigma ,𝐪0}{}}c_{𝐤+𝐪,\sigma }^{}c_{𝐤,\sigma }\left(a_𝐪+a_𝐪^{}\right).`$ (10)
Following Fröhlich,Frohlich52  we seek a canonical transformation $`\stackrel{~}{H}=e^{iS}He^{iS}`$ which eliminates electron-phonon coupling to first order in $`D`$. This is accomplished by

$$iS=\frac{D}{\sqrt{L}}\underset{𝐤,\sigma ,𝐪0}{}\left(\frac{c_{𝐤+𝐪,\sigma }^{}c_{𝐤,\sigma }a_𝐪}{\epsilon _{𝐤+𝐪}\epsilon _𝐤\omega _𝐪}+\frac{c_{𝐤+𝐪,\sigma }^{}c_{𝐤,\sigma }a_𝐪^{}}{\epsilon _{𝐤+𝐪}\epsilon _𝐤+\omega _𝐪}\right),$$
(11)
since, as one verifies easily, $`[iS,H_0]=H^{}`$. The effective Hamiltonian can therefore be written as $`\stackrel{~}{H}=H_0+\stackrel{~}{H}^{}`$ where, to second order in $`D`$,

$`\stackrel{~}{H}^{}={\displaystyle \frac{1}{2}}[iS,H^{}]`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{𝐤,𝐤^{},\sigma ,\sigma ^{}}{}}{\displaystyle \underset{𝐪0}{}}{\displaystyle \frac{D^2\omega _𝐪}{(\epsilon _{𝐤+𝐪}\epsilon _𝐤)^2\omega _𝐪^2}}c_{𝐤^{}𝐪,\sigma ^{}}^{}c_{𝐤^{},\sigma ^{}}c_{𝐤+𝐪,\sigma }^{}c_{𝐤,\sigma }`$
$`+`$ $`{\displaystyle \frac{D^2}{2L}}{\displaystyle \underset{𝐪0}{}}{\displaystyle \underset{𝐪^{}0}{}}{\displaystyle \underset{𝐤,\sigma }{}}(c_{𝐤+𝐪,\sigma }^{}c_{𝐤𝐪^{},\sigma }c_{𝐤+𝐪+𝐪^{},\sigma }^{}c_{𝐤,\sigma })\left[{\displaystyle \frac{a_𝐪(a_𝐪^{}+a_𝐪^{}^{})}{\epsilon _{𝐤+𝐪}\epsilon _𝐤\omega _𝐪}}+{\displaystyle \frac{a_𝐪^{}^{}(a_𝐪^{}+a_𝐪^{}^{})}{\epsilon _{𝐤+𝐪}\epsilon _𝐤\omega _𝐪}}\right].`$
The expectation value of $`\stackrel{~}{H}^{}`$ in the product of a BCS state with a phonon state diagonalizing $`H_0`$, $`|\stackrel{~}{\psi }`$, simplifies to

$$E=\stackrel{~}{\psi }|\stackrel{~}{H}^{}|\stackrel{~}{\psi }=\frac{1}{L}\underset{𝐤,𝐤^{}}{}V_{𝐤,𝐤^{}}b_𝐤^{}^{}b_𝐤+\frac{D^2}{L}\underset{𝐤,𝐤^{},\sigma }{}\frac{f_{𝐤\sigma }(1f_{𝐤^{}\sigma })}{\epsilon _𝐤\epsilon _𝐤^{}\omega _{𝐤𝐤^{}}}+\underset{𝐪0}{}\delta \omega _𝐪n(\omega _𝐪),$$

where $`n(\omega )`$ is the Bose distribution function. The first term describes the effective phonon mediated electron-electron interaction in the Cooper channel with $`V_{𝐤,𝐤^{}}=2D^2\omega _{𝐤𝐤^{}}/[(\epsilon _𝐤\epsilon _𝐤^{})^2\omega _{𝐤𝐤^{}}^2]`$, the second term is the standard second-order perturbation theory correction to the ground state energy of the electron system, and the last term with $`\delta \omega _𝐪=2D^2\chi _1(𝐪,\omega _𝐪)`$ is the phonon energy renormalization.

In order to facilitate a numerical solution of the gap equation Eq. (7) for phonon mediated pairing, we have made use of a Cooper channel interaction

$$V_{𝐤,𝐤^{}}=D^2\mathrm{Re}\underset{\sigma =\pm }{}\frac{\sigma }{\epsilon _𝐤\epsilon _𝐤^{}+\sigma (\omega _{𝐤𝐤^{}}+i\mathrm{\Gamma })},$$

regularized by a finite phonon lifetime $`\mathrm{\Gamma }`$. The BCS gap equation has been solved numerically on the same $`t`$-$`t^{}`$ square lattice as for the Hubbard model. The superconducting gap obtained numerically for an Einstein mode with frequency $`\omega _0/t=0.2`$ and $`\mathrm{\Gamma }/t=0.02`$ for electronic parameters $`t^{}/t=0.45`$ and $`\rho =1.0`$ and electron-phonon coupling $`D/t=0.51`$ is shown in Fig. 9. Note that the gap changes sign in the vicinity of the phonon energy $`\omega _0`$. This nontrivial energy dependence of $`\mathrm{\Delta }_𝐤`$ translates into a feature of the tunneling density of states in the same energy region, see Fig. 10. It should be pointed out that although both of the above results are in qualitative agreement with the Eliashberg theory,Scalapino69  there is one important difference of the present approach with respect to Ref. Scalapino69, . Namely, within the Eliashberg theory the features of $`\mathrm{\Delta }(\omega )`$ and $`N(\omega )`$ occur at energy $`\omega _0+\mathrm{\Delta }(0)`$ and not $`\omega _0`$ as in our case. This difference is presumably due to the fact that within the canonical transformation approach, it is the bare electron energies which enter the energy denominators in Eq. (11) and not the renormalized energies as in the Eliashberg theory. However, in a truly weak coupling theory this difference is marginal.


text:
# Non-Gaussian PDFs from Maximum-Entropy-principle considerations

## Abstract

In this work we develop on the recently suggested concept of superstatistics \[C. Beck and E.G.D. Cohen, Physica A 322, 267 (2003)\], face the problem of devising a viable way for estimating the correct statistics for a system in absence of sufficient knowledge of its microscopical dynamics, and suggest to solve it through the Maximum Entropy Principle. As an example, we deduce the Probability Distribution Function for velocity fluctuations in turbulent fluids, which is slightly different from the form suggested in \[C. Beck, Phys. Rev. Lett. 87, 180601 (2001)\].

Nonextensive statistical mechanics (NESM) has been raising considerable interest along these last 15 years. However, while there is an universal consensus about its most famous prediction, i.e. the ubiquitous existence in nature of power-law Probability Distribution Functions (PDFs), the soundness of the theoretical foundation on which is based-the generalized definition of entropy tsallis -has been questioned by several authors. We quote here, e.g., the paper nauenberg  which dealt in particular with contradictions between theory’s predictions and thermodynamical constraints.
The paper sattinepl  started from these critiques but focussed more on the information–theoretic aspect, suggesting that Tsallis’ entropy should be regarded merely as a practical tool for doing predictions in presence of a reduced amount of information about the system. Indeed, the need of resorting to modified definitions of the informational entropy when the knowledge about the states of the system is insufficient, is well known in statistics and has been extensively pointed out in luzzi1  and references therein.
Quite recently, Sattin and Salasnich pre –starting from an earlier work by Beck beckprl –and, in a more formal and rigorous way, Beck and Cohen superstatistics  demonstrated, without direct reference to any definition of entropy, that Tsallis’ statistics is just a particular case of an infinite class of statistics (hence the name “superstatistics”). All the elements of the class are characterized by one or more parameters, and reduce to ordinary Maxwell-Boltzmann statistics for particular values of the parameters. Beck and Cohen started from a model of dynamical system, a Brownian particle moving according to a Langevin equation, where the noise ($`\sigma `$) and friction ($`\gamma `$) terms are allowed to fluctuate, to arrive to the famous equation

$$P(E)=K_0^{\mathrm{}}e^{\beta E}f(\beta )𝑑\beta $$
(1)
Here, $`P(E)`$ is the PDF for the system of being in the state of “energy” $`E`$, $`\beta `$ is a fluctuating parameter (generalized inverse temperature), which in the present formulation is function of $`\sigma `$, $`\gamma `$, and $`f(\beta )`$ is the PDF for the realization of the particular value $`\beta `$; $`e^{\beta E}`$ is the usual Boltzmann factor. Eq. (1) is more general than suggested by the Brownian particle model. Indeed, it can be written for any system $`\mathrm{\Sigma }`$ interacting with a fluctuating environment $`B`$, hence $`P(E)`$ is a weighted average of the standard Boltzmann statistics over different realizations of the interaction between the system and its environment, quantified by $`f(\beta )`$. Notice that the existence of finite fluctuations in the environment is strictly correlated with finite-size-effects of the environment itself; indeed, the result that Tsallis statistics can arise within the framework of a system interacting with a finite thermal bath was recently reported by Aringazin and Mazhitov aringazin , and it had already been suggested much earlier by Plastino and Plastino plastino .
The physical content of the theory, thus, shifts from the entropic index $`q`$ to the PDF $`f(\beta )`$. Of course, the explicit expression for $`f(\beta )`$ must vary for any single problem, constrained just by some rather intuitive criteria (normalizability, etc …); Beck and Cohen give several possible examples of functions which are potential candidate for $`f(\beta )`$. However, a simple criterion able to guide the user towards a plausible functional form for $`f`$, lacking a more detailed knowledge of the underlying microscopical details, would be very satisfying. But such a criterion is readily available: it is the well known Jaynes’ Maximum Entropy Principle. In this case, the system about which we have not a proper knowledge is no longer $`\mathrm{\Sigma }`$ but the environment $`B`$ itself, and $`f(\beta )`$ is a measure of the probability of $`B`$ of occupying a state in an abstract one-dimensional space parameterized by $`\beta `$.
Lacking any further information, the most probable realization of $`f(\beta )`$ will be the one that maximizes (Shannon) entropy $`S(f)=f\mathrm{ln}fd\beta `$ with suitable constraints.
In this paper we present a straightforward application of this principle to an important example, namely the PDF of velocity fluctuations in turbulent flows bls . This accurate experimental measurement is thought to be one of the strongest evidences in support of Tsallis’ theory, since the empirical PDF appear very well matched by power-laws. In this work we suggest instead that the true curve can be very close numerically, but rather different in its analytical expression, from a power-law. Indeed, further experimental investigations of turbulent flows are now suggesting the existence of small deviation from pure power-laws (experiments of Jung and Swinney, cited in superstatistics ). Recent work by Aringanzin and Mazhitov ar1  deal with the attempt of theoretically recovering the new accurate experimental distributions for fluid particle accelerations. In this work, instead, we will consider, only fluid velocity differences.
As a first test case, let us suppose of knowing anything about $`f`$ but the average value of the inverse temperature

$$<\beta >=\beta f(\beta )𝑑\beta .$$
(2)
The most probable $`f`$ extremizes the functional

$$F=f\mathrm{ln}fd\beta \lambda f\beta 𝑑\beta .$$
(3)
and the solution is the Boltzmann-like function

$$f(\beta )=f_0\mathrm{exp}(\lambda \beta )$$
(4)
($`\lambda `$ is the Lagrange multiplier). Such an expression was obtained, by example, in batanov , from which the original idea of entropy maximization was taken. It reads slightly different from the Gamma (or $`\chi ^2`$) distribution originally devised by Beck beckprl  to recover power-law PDF. If we replace $`f`$ of (4) into (1) we get

$$P(E)\frac{1}{E+\lambda },$$
(5)
which is not a physically acceptable solution, since $`P(E)`$ must be normalizable.
In order to make a step beyond this rough scheme it is necessary to input more information within the model. We do it starting from the same premises as Beck’s beckprl  but diverging at just the next step: infact, he defines $`\beta ϵ_r\tau `$, where $`ϵ_r`$ is the energy dissipation rate of the fluid on microscopical scale, and $`\tau `$ the typical energy transfer time. We note that $`ϵ_r\tau `$ has units of \[energy\], not of \[energy\]-1, and a more intuitive way of writing this relation should be $`\beta (ϵ_r\tau )^1`$. However, the inverse of a sum of squares of random variables does not yield a $`\chi ^2`$-distributed random variable, hence Beck is forced to re-establish the correct dimensions by multiplying by the constant $`\mathrm{\Lambda }`$, with the units of \[speed\]4. Thus, a characteristic speed $`\mathrm{\Lambda }^{1/4}`$ has entered the calculations, whose physical interpretation remains obscure. Of course, Beck was forced to do this assumption because of the sought agreement with Tsallis’ theory. Now, we are free from this constraint, and can allow for more natural choices, although we must agree that a certain arbitrariness in the choice of the physically meaningful variables is unavoidable within this context.
To start with, we shift from $`\beta `$ parameter to $`T=1/\beta `$. $`T`$, which has units of \[energy\] is a more convenient variable, as explained above. We can write the equivalent of Eq. (1) in terms of the new parameter; furthermore, in order to adhere to existing literature, it is convenient using a generalized velocity instead of energy: $`E=u^2/2`$. Therefore Eq. (1) becomes

$$p(u)=𝑑T\sqrt{\frac{1}{2\pi T}}\mathrm{exp}\left(\frac{u^2}{2T}\right)g(T).$$
(6)
Again, we follow Beck’s recipe: the parameter $`T`$ is written in terms of fluctuating Kolgomorov velocities $`u_i`$

$$T=\frac{T_0}{3}\underset{i=1,2,3}{}u_i^2$$
(7)
where $`T_0`$ is a constant. We stress that the r.h.s. of the previous equation and that of Eq. (28) in beckprl  are the same, although the l.h.s.’s are each the inverse of the other. This is due to our discarding of the artificial constant $`\mathrm{\Lambda }`$. $`T_0`$, conversely, has a straightforward physical interpretation as average thermal energy.
The PDF $`g(T)`$ can be rewritten as a function of $`u`$’s: $`g^{}(u_1,u_2,u_3)g^{}(𝐮)g(T)`$. The Maximum-Entropy principle imposes of extremizing

$$F=g^{}(𝐮)\mathrm{ln}(g^{}(𝐮))\frac{1}{u_M^2}(u_1^2+u_2^2+u_3^2)g^{}(𝐮)$$
(8)
where $`1/u_M^2`$ is the Lagrange multiplier corresponding to the fact that we fix the average kinetic energy. The solution is, obviously, $`g^{}(𝐮)=C\mathrm{exp}((u_1^2+u_2^2+u_3^2)/u_0^2)`$. From here, reversing to $`T`$ variable,

$$g(T)=KT^{1/2}\mathrm{exp}\left(\frac{3T}{2T_0}\right).$$
(9)
The term $`T^{1/2}`$ comes from the volume element, $`C,K`$ are normalization constants and we have rescaled velocities such that $`u_M^22/3`$. Notice that Eq. (9) is a Gamma distribution and could be obtained straightforwardly from Eq. (7) by assuming from the start that the $`u_i`$ were normal random variables, just as done by Beck.
Let us now replace Eq. (9) into (6): we get

$$p(u)=K^{}_0^{\mathrm{}}\mathrm{exp}\left(\frac{u^2}{2T}\frac{3T}{2T_0}\right)𝑑T.$$
(10)
The explicitly normalized solution reads

$$p(u)_{PW}=\frac{1}{\widehat{u}\pi }\frac{u}{\widehat{u}}K_1\left(\frac{u}{\widehat{u}}\right)\left(\widehat{u}=\sqrt{\frac{T_0}{3}}\right)$$
(11)
and $`K_1`$ is the Bessel $`K`$ function of order one. This result appears rather different from usual power-laws. Indeed, we chose to plot in Fig. (1) this curve together with the best fitting curve found by Beck for the velocity PDF:

$$p(u)_{\text{Beck}}=\frac{1}{Z_q}\frac{1}{\left(1+(q1)\stackrel{~}{\beta }C|u|^{2\alpha }\right)^{1/(q1)}}(q1.1,\alpha 0.9).$$
(12)
(See dashed line in Fig. 1 of beckprl ). On the whole, the two curves match rather closely. Some differences appear at low $`u`$’s.

Since we have not available the original experimental data $`p_{\text{expt}}`$, we cannot directly compare on them the goodness of our fit. A good insight comes, however, by plotting the relative difference $`\mathrm{\Delta }=(p_{\text{Beck}}p_{PW})/p_{\text{Beck}}`$ (Fig. 2, $`p_{PW}`$ is our solution, given in Eq. 11). The quantity $`\mathrm{\Delta }_1=(p_{\text{Beck}}p_{\text{expt}})/p_{\text{Beck}}`$ appears plotted in Fig. 2 of bls . If $`\mathrm{\Delta }`$ and $`\mathrm{\Delta }_1`$ agree, $`p_{PW}=p_{\text{expt}}`$. It appears that our fit slightly overestimates experiment at $`u0`$, but for the same amount as $`p_{\text{Beck}}`$ does underestimate it. On the whole, the agreement with experiment is remarkable. We remark that, as only adjustable parameter, we used the hypothesis that $`T`$ is of the form (7) with the index ranging from one to three. The equivalent of $`\alpha `$ parameter used by Beck does not enter our calculations. Thus, we have realized an economy in our way of modelling the data.

We call the reader’s attention to the fact that, in his Fig. 1, Beck studied at once (i) velocity spatial differences, and (ii) accelerations. We, instead, considered only the former quantity, and it is possible to see that the latter cannot be reproduced by the present treatment, even allowing for a varying number of Kolmogorov velocities in Eq. (7). This must be traced back to the fact that now $`g(T)`$ is no longer a good weight function, but we must define an equivalent of $`T`$, defining the “average acceleration” of the system.
Finally, we point to two important issues: I) we have the asymptotic trend $`xK_1(x)x^{1/2}\mathrm{exp}(x),(x\mathrm{})`$, that is, we have not power-law decay. Notwithstanding this, our curve nicely fits data that have been previously considered as stemming from a power-law PDF. The reason lies in the finite $`u`$-range sampled and therefore calls for an inherent ambiguity in this kind of studies, related to finite experimental scans: unless one is sure of being investigating the true asymptotic region, one can never be completely confident about the fitting curve used.
II) As reported above, all superstatistics must collapse to the single Boltzmann statistics when $`q1`$, that is, differences between different models are at least order $`𝒪(q1)`$. Since, in this case, we are dealing with a parameter not far from unity, $`q1.1`$, it is to be expected that any two reasonable models would give close results. (The author wishes to thank C. Beck for calling these two points to his attention).

In conclusion, the use of the concept of superstatistics together with Maximum-Entropy principle appears to be an efficient way of estimating statistical properties in general systems using a minimal amount of information.

###### Acknowledgements.

The author wishes to thank A. Rapisarda, C. Beck, and R. Luzzi for useful comments and suggestions.