text:
# Meson Form Factors and Non-Perturbative Gluon Propagators

## I Introduction

The study of the infrared (IR) limit of the constant coupling of the Quantum Chromodynamics (QCD) had attracted much attention in the last years, as well as the infrared finite form of the gluon propagator. The behavior of the running coupling constant claimed for several theoretical and experimental studies indicates a frozen value in this kinematical region amn02 . The question of the infrared form of Green functions of QCD is still a controversial field of research (for a review see avs01 ; rwja ). Due to the different methods employed to analyze the Green functions, for example, the Dyson-Schwinger equations (DSE) and simulations in the lattice field theory (LFT) (as well as the approximations used to avoid several difficulties found in these methods) we have distinct different forms of Green functions. However, as we will show, the combination of the results from LFT and solutions of DSE restricts the form of the gluon propagator.

In a recent paper amn02 , Aguilar, Mihara and Natale perform a study of the running coupling constant in the infrared region, using the dependence of the meson form factors on the coupling constant. To calculate the pion form factor, a mixed framework between perturbative and non-perturbative physics is employed with the following assumptions: freezing of the coupling constant and finite gluon propagator in the IR region, although using the perturbative scheme of calculation of the meson form factor lb79 .

Nevertheless, the use of a frozen running coupling constant to describe QCD exclusive processes, particularly the meson form factor, is not a novelty. Ji and Amiri ja90  and Brodsky et al. bjpr98  use a modified coupling constant to study the form factor and the related reaction $`\gamma \gamma \pi ^+\pi ^{}`$. This modified coupling constant is based on the work by Cornwall jmc82  in which a gauge independent set of Dyson-Schwinger equations is solved, giving a coupling constant finite in the infrared and a gluon propagator with a running mass.

However, there are another approaches to include infrared contributions. In mrt , Maris and collaborators perform a full non-perturbative calculation of the meson form factor. The amplitude for the vertex $`\pi \pi \gamma `$ (or $`KK\gamma `$), which is related with the form factor, is obtained using the following ingredients: the Bethe-Salpeter amplitudes for the scattering kernel $`q\overline{q}`$ (solutions of the Bethe-Salpeter equation); the solution of the DSE for the quark propagator in the rainbow truncation and an ansatz for the quark-photon vertex based in a Bethe-Salpeter equation in the ladder truncation. This approach gives a good description of the experimental data for the low momentum region for both form factors (pion and kaon) as well as another observables (see last paper of mrt ). In another work kcj01 , the authors employed a light-front Bethe-Salpeter model, where the quark, instead of the gluon as the model above, has a dynamical mass, giving a good agreement with the experimental data for the pion form factor.

The meson form factor in turn, has many attempts of description in pure perturbative QCD vlhp . For example, Stefanis ngs  uses an unified factorization scheme that includes logarithmic corrections (which origin is the gluon emission) and power correction whose origin is non-perturbative. Melić et al. mnp99  calculate the pion form factor in next-to-leading order (NLO) in perturbation theory. Yeh yeh02  also made a NLO calculation and found a fittable expression for the pion form factor which includes non-perturbative effects. However, these attempts have only validity where the perturbation theory is valid, in other words, when the momentum transfer is large.

In this article, we extend the results of amn02  for the pion form factor for a different set of coupling constants and propagators and compare them with more recent experimental data. Also the robustness of the model is tested with the kaon, a more massive meson than the pion. As we will show, the model employed describes both, pion and kaon form factors with good agreement with the available data for the kinematical region of exchange of low momenta. The article is organized as follows: in the section II, we will review the model used to calculate the meson form factor. Next, in section III, we present the modified meson form factor calculation scheme with the main features of the non-perturbative gluons and the frozen running coupling constant and their relation. Finally, we present the results and conclusions of the work.

## II The Meson Form Factor Calculation

The meson form factor is given by the factorized expression in perturbative QCD (pQCD) form lb79 ; blk83 ; fld (see figure 1)

$$F_M(Q^2)=_0^1𝑑x_0^1𝑑y\varphi _M^{}(y,\stackrel{~}{Q}_y)T_H(x,y,Q^2)\varphi _M(x,\stackrel{~}{Q}_x),$$
(1)
where $`\stackrel{~}{Q}_z=\mathrm{min}(z,1z)Q`$ ($`z=x,y`$) and $`Q`$ is the 4-momentum transfer by the photon. The equation above can be seen as a convolution of the initial and final states, represented by the quark amplitude distributions $`\varphi _M(z,\stackrel{~}{Q}_z)`$ (obtained in a non-perturbative calculation) with a hard (perturbative) scattering amplitude, $`T_H(x,y,Q^2)`$. The quark amplitude distribution is interpreted as the amplitude to find the quark or antiquark within the meson with fractional momentum $`z`$ or $`1z`$, respectively.

The function $`T_H(x,y,Q^2)`$ was calculated in leading order (LO) in bt87  and in next leading order (NLO) in mnp99 . The LO expression is given by (fig. 1):

$`T_H(x,y,Q^2)`$ $`=`$ $`{\displaystyle \frac{64\pi }{3Q^2}}\{{\displaystyle \frac{2}{3}}{\displaystyle \frac{\alpha _s[(1x)(1y)Q^2]}{(1x)(1y)}}`$ (2)
$`+{\displaystyle \frac{1}{3}}{\displaystyle \frac{\alpha _s(xyQ^2)}{xy}}\},`$
where $`\alpha _s(Q^2)`$ is the running coupling constant.

The quark amplitude distributions $`\varphi _M(x,\stackrel{~}{Q}^2)`$ are usually calculated in a certain energy scale ($`Q_0`$) by QCD sum rules cz84  and then evaluated in an arbitrary energy scale solving a Bethe-Salpeter-type evolution equation lb79 ; fld , which solution is

$$\varphi _M(x,Q^2)=x(1x)\underset{n=0}{\overset{\mathrm{}}{}}C_n^{(3/2)}(2x1)\left[\frac{\alpha _s(Q^2)}{\alpha _s(Q_0^2)}\right]^{d_n}\varphi _n^{(M)}(Q_0^2),$$
(3)
where $`C_n^{(3/2)}(2x1)`$ are the Gegenbauer orthogonal polynomials ast ;

$$d_n=\frac{2A_n^{\mathrm{NS}}}{\beta _0},$$
(4)
where $`A_n^{\mathrm{NS}}`$ is the non-singlet anomalous dimension given by

$$A_n^{\mathrm{NS}}=\frac{4}{3}\left[\frac{1}{2}+\frac{1}{(n+1)(n+2)}2\underset{j=2}{\overset{n+1}{}}\frac{1}{j}\right],$$
(5)
and $`\beta _0=11\frac{2}{3}n_f`$ ($`n_f`$ is the number of flavors). To obtain the factors $`\varphi _n^{(M)}(Q_0^2)`$, given by

$$\varphi _n^{(M)}(Q_0^2)=\frac{4(2n+3)}{(n+1)(n+2)}_0^1𝑑xC_n^{(3/2)}(2x1)\varphi _M(x,Q_0^2),$$
(6)
we use the quark amplitude distributions in a given scale of energy. From cz84 , we have at $`Q_0=500\mathrm{MeV}`$ for the pion and the kaon, respectively,

$$\varphi _\pi (x,Q_0^2)=\frac{30f_\pi }{2\sqrt{3}}x(1x)(2x1)^2,$$
(7a)
$$\varphi _K(x,Q_0^2)=\frac{30f_K}{2\sqrt{3}}\left[0.6(2x1)^2+0.25(2x1)^3+0.08\right],$$
(7b)
where the normalization factor $`f_M`$ is given by

$$_0^1𝑑x\varphi _M(x,Q_0^2)=\frac{f_M}{2\sqrt{3}},$$
(8)
and thereby $`f_\pi =93\mathrm{MeV}`$, $`f_K=112\mathrm{MeV}`$.

Using the orthogonality relation of the Gegenbauer polynomials, it is easy to find that the only non zero $`\varphi _n^{(M)}(Q_0^2)`$’s for the pion are

$$\varphi _0^{(\pi )}(Q_0^2)=6\frac{f_\pi }{2\sqrt{3}},\varphi _2^{(\pi )}(Q_0^2)=4\frac{f_\pi }{2\sqrt{3}},$$

and for the kaon

$`\varphi _0^{(K)}(Q_0^2)=6{\displaystyle \frac{f_K}{2\sqrt{3}}}`$ , $`\varphi _1^{(K)}(Q_0^2)={\displaystyle \frac{30}{28}}{\displaystyle \frac{f_K}{2\sqrt{3}}},`$ (9a)
$`\varphi _2^{(K)}(Q_0^2)={\displaystyle \frac{12}{5}}{\displaystyle \frac{f_K}{2\sqrt{3}}}`$ , $`\varphi _3^{(K)}(Q_0^2)={\displaystyle \frac{3}{7}}{\displaystyle \frac{f_K}{2\sqrt{3}}},`$ (9b)
therefore the exponents $`d_n`$ are

$`d_0=0`$ , $`d_1={\displaystyle \frac{32}{996n_f}},`$ (10a)
$`d_2={\displaystyle \frac{50}{996n_f}}`$ , $`d_3={\displaystyle \frac{314}{49530n_f}}.`$ (10b)
Finally, the quark amplitude distribution is, for the pion,

$$\varphi _\pi (x,Q^2)=\frac{f_\pi }{2\sqrt{3}}x(1x)\left\{6+\left[30(2x1)^26\right]\left(\frac{\alpha _s(Q^2)}{\alpha _s(Q_0^2)}\right)^{d_2}\right\},$$
(11)
and for the kaon,

$`\varphi _K(x,Q^2)={\displaystyle \frac{f_K}{2\sqrt{3}}}x(1x)`$
$`\times \{6+{\displaystyle \frac{45}{14}}(2x1)\left({\displaystyle \frac{\alpha _s(Q^2)}{\alpha _s(Q_0^2)}}\right)^{d_1}+`$
$`+{\displaystyle \frac{6}{5}}\left[15(2x1)^23\right]\left({\displaystyle \frac{\alpha _s(Q^2)}{\alpha _s(Q_0^2)}}\right)^{d_2}+`$
$`+{\displaystyle \frac{3}{14}}[35(2x1)^315(2x1)]\left({\displaystyle \frac{\alpha _s(Q^2)}{\alpha _s(Q_0^2)}}\right)^{d_3}\}.`$ (12)
However, Dziembowski and Mankiewicz dm87  use the constituent-quark model to calculate the quark amplitude distribution for the pion, obtaining

$$\varphi _\pi ^{(\mathrm{DM})}(x)=N\mathrm{exp}\left[\frac{m^2}{8x(1x)\beta ^2}\right]\left\{\frac{(xM+m)\left[(1x)M+m\right]}{4\beta ^4}2x(1x)\right\},$$
(13)
where $`N`$ is determined by the normalization condition Eq.(8)($`N=0.622`$), $`M`$ is the spin-averaged meson mass ($`M=614.4`$ MeV for the pion), $`m`$ is the constituent quark mass ($`m=330`$ MeV) and $`\beta `$ is a Gaussian parameter chosen as $`460`$ MeV. We will test the changes in the pion form factor when this quark distribution is employed.

## III The modified meson form factor

In amn02 , the modifications introduced in the pion form factor by the change of the usual running coupling amplitude by the frozen infrared one have been tested, as well as the case of the change of the gluon propagator. In this case, the hard scattering amplitude can be written as

$$\stackrel{~}{T}_H(x,y,Q^2)=\frac{64\pi }{3}\left\{\frac{2}{3}\stackrel{~}{\alpha }_s(\widehat{k}^2)D(\widehat{k}^2)+\frac{1}{3}\stackrel{~}{\alpha }_s(\widehat{p}^2)D(\widehat{p}^2)\right\},$$
(14)
where $`\widehat{k}^2=(1x)(1y)Q^2`$, $`\widehat{p}^2=xyQ^2`$ and $`D(Q^2)`$ is the gluon propagator, and $`\stackrel{~}{\alpha }_s`$ is the modified running coupling constant which expression is given below.

Therefore, the expression for the meson form factor from Eq.(1) and (14) is

$$\stackrel{~}{F}_M(Q^2)=_0^1𝑑x_0^1𝑑y\stackrel{~}{\varphi }_M^{}(y,\stackrel{~}{Q}_y)\stackrel{~}{T}_H(x,y,Q^2)\stackrel{~}{\varphi }_M(x,\stackrel{~}{Q}_x),$$
(15)
where $`\stackrel{~}{\varphi }_M`$ is the meson wave function with the modified running coupling constant.

The infrared form of the gluon propagator is still a controversial aspect (for a review for solutions and methods see rwja  and references therein). The general formula for the gluon propagator is given by (in the Landau gauge):

$$D_{\mu \nu }^{ab}(q^2)=\delta ^{ab}\left(\delta _{\mu \nu }\frac{q_\mu q_\nu }{q^2}\right)D(q^2),$$
(16)
where the usual perturbative propagator is given by $`D(q^2)=1/q^2`$ which diverges when $`q^20`$. In the literature, there are several different forms for the gluon propagator, with different behaviors for the infrared region: finite, zero and more divergent that $`1/k^2`$, each one with its advantages and inconveniences. The reasons of these different behaviors are the methods and approximations employed to obtain the gluon propagator. The more popular methods are the Schwinger-Dyson equations and the simulations from lattice field theory. The most recent results in the last method fdrb01 ; lrg02  discard the solution of the type $`1/k^4`$ and therefore we will not use it in this work. In spite of the controversy, the non-perturbative gluon propagator was used successfully in many phenomenological applications as, for example, proton-proton scattering hkpr  and elastic production of vector mesons ghnw .

In this work, we use the following gluon propagators:

* calculated by Cornwall jmc82  using the gauge independent pinch technique, given by

$$D_\mathrm{C}(q^2)=\frac{1}{q^2+M_g^2(q^2)},$$
(17)
where $`M_g^2(q^2)`$ is a dynamical gluon mass term,

$$M_g^2(q^2)=m_g^2\left\{\frac{\mathrm{ln}\left(\frac{q^2+4m_g^2}{\mathrm{\Lambda }_{\mathrm{QCD}}^2}\right)}{\mathrm{ln}\left(\frac{4m_g^2}{\mathrm{\Lambda }_{\mathrm{QCD}}^2}\right)}\right\}^{12/11}$$
(18)
with $`m_g^2=500\pm 200\mathrm{MeV}`$. With this propagator, amn02  obtains the best description for the pion form factor. Its features include the correct ultraviolet behavior (according to the renormalization group) and the dynamically generated mass $`M_g^2(q^2)`$.
* from Häbel et al. hkrsw90 , calculated using an approach that employs the same features of the perturbative theory, which implies a simplified set of Schwinger-Dyson equations, which solution is given by

$$D_\mathrm{H}(q^2)=\frac{1}{q^2+\frac{b^2}{q^2}},$$
(19)
We should note that this propagator goes to zero when $`q^20`$, unlike the other propagators employed in this work, which have a non-zero value in this region.
* Alkofer et al. ahs  solve a set of Dyson-Schwinger equations for the QCD Green’s functions and the gluon propagator found by this manner can be fitted amn02  by

$$D_\mathrm{A}(q^2)=\frac{bq^2}{q^4+c^2},$$
(20)
where $`b=3.707`$ and $`c=0.603`$.
* Atkinson and Bloch abw  solve a set of Dyson-Schwinger equations for the gluon and ghost propagator, improving the solution found by ahs . The asymptotic infrared behavior of the gluon propagator is

$$D_{\mathrm{AB}}^{\mathrm{IR}}(q^2)\frac{1}{q^2}\left\{A_0(q^2)^{2\kappa }\left(1+\underset{\lambda =1}{\overset{3}{}}f_\lambda a_\lambda (q^2)^{\lambda \rho }\right)\right\},$$
(21)
where $`A_0=1`$, $`\kappa =0.769475`$, $`a_1=10.27685`$, $`\rho =1.96964`$, $`f_1=1`$, $`f_2=0.408732`$, $`f_3=0.761655`$.
* from Gorbar and Natale gn00 , that use the operator product expansion (OPE) to relate the gluon and quark propagators with their respective condensates, the gluon propagator obtained is

$$D_{\mathrm{GN}}(q^2)=\frac{1}{q^2+\mu _g^2\theta (\chi \mu _g^2q^2)+\frac{\mu _g^4}{q^2}\theta (q^2\chi \mu _g^2)},$$
(22)
where $`\mu _g^2=\mathrm{\hspace{0.17em}0.61149}\mathrm{GeV}`$ is a parameter fixed by the gluon condensate gn00  and $`\chi =`$0.9666797. We should note that this propagator interpolates two different propagators: the pure massive and one similar to the Häbel solution.

The modification of the gluon propagator is closely related with the modification of the running coupling constant in the infrared region. As calculated by Cornwall jmc82 , the coupling constant is frozen in the low momentum by the addition of the massive term, Eq.(18),

$$\alpha _s^{(\mathrm{C})}(q^2)=\frac{4\pi }{\beta _0\mathrm{ln}\left(\frac{q^2+4M_g^2(q^2)}{\mathrm{\Lambda }_{\mathrm{QCD}}^2}\right)},$$
(23)
where $`\mathrm{\Lambda }_{\mathrm{QCD}}^2`$ is the QCD scale parameter and $`\beta _0=11\frac{2}{3}n_f`$ is the first coefficient of the beta function and $`n_f`$ is the number of flavors. In the case of the other propagators, we consider the terms in the denominator as massive terms and we will substitute it in the Eq.(23) as well as in the Cornwall propagator, as shown below

* The Häbel propagator gets

$$\alpha _s^{(\mathrm{H})}(q^2)=\frac{4\pi }{\beta _0\mathrm{ln}\left(\frac{q^2+M_\mathrm{H}^2(q^2)}{\mathrm{\Lambda }_{\mathrm{QCD}}^2}\right)},$$
(24)
with

$$M_\mathrm{H}^2(q^2)=\frac{b^2}{q^2}.$$
* The Gorbar-Natale propagator gets

$$\alpha _s^{(\mathrm{GN})}(q^2)=\frac{4\pi }{\beta _0\mathrm{ln}\left(\frac{q^2+M_{\mathrm{GN}}^2(q^2)}{\mathrm{\Lambda }_{\mathrm{QCD}}^2}\right)},$$
(25)
with

$$M_{\mathrm{GN}}^2(q^2)=\mu _g^2\theta (\chi \mu _g^2q^2)+\frac{\mu _g^4}{q^2}\theta (q^2\chi \mu _g^2).$$

However, in the case of the propagators given by the equations (20) and (21), the coupling constant comes from the solutions of the Dyson-Schwinger equations for the Green’s functions. For the case of the solution of the Alkofer et al., the coupling constant can be fitted by the following formula amn02

$`\alpha _s^\mathrm{A}(q^2)=\{\begin{array}{cc}\alpha _s^\mathrm{a}:\hfill & q^2<0.31\mathrm{GeV}^2\hfill \\ \alpha _s^\mathrm{b}:\hfill & 0.31<q^2<1.3\mathrm{GeV}^2,\hfill \\ \alpha _s^\mathrm{c}:\hfill & q^2>1.3\mathrm{GeV}^2\hfill \end{array}`$ (29)
with

$`\alpha _s^\mathrm{a}(q^2)`$ $`=`$ $`0.2161+9.2621\mathrm{exp}\left(2{\displaystyle \frac{(q^20.0297)^2}{(0.6846)^2}}\right),`$
$`\alpha _s^\mathrm{b}(q^2)`$ $`=`$ $`1.4741+8.6072\mathrm{exp}\left({\displaystyle \frac{q^20.1626}{0.3197}}\right),`$
$`\alpha _s^\mathrm{c}(q^2)`$ $`=`$ $`{\displaystyle \frac{1.4978}{\mathrm{ln}(1.8488q^2)}}.`$ (30)
The solution of Atkinson and Bloch abw  is given analytically in the asymptotic regions

$$\alpha _s^{\mathrm{AB}}(q^2)\stackrel{\mathrm{IR}}{}4\pi \nu \left[1\underset{k=1}{\overset{3}{}}b_k\left(\frac{q^2}{\mathrm{\Omega }^2}\right)^{k\rho }\right],$$
(31)
where $`\nu =0.912771`$, $`\mathrm{\Omega }^2=0.1864754`$, $`b_1=1`$, $`b_2=0.760753`$, $`b_3=0.370785`$, and

$$\alpha _s^{\mathrm{AB}}(q^2)\stackrel{\mathrm{UV}}{}\frac{4\pi }{4\mathrm{log}(\frac{q^2}{\mathrm{\Lambda }_{\mathrm{QCD}}^2})},$$
(32)
where $`\mathrm{\Lambda }_{\mathrm{QCD}}^2=0.06802\mathrm{GeV}^2`$.

The procedure is the following: using the Eq.(15), we substitute in the hard scattering amplitude $`\stackrel{~}{T}_H`$ the different gluon propagators in the order above as well as the running coupling constant, $`\stackrel{~}{\alpha }_s`$, to verify the changes in the form factor. A comparison with the perturbative fit of Yeh yeh02 , as well as, with the full non-perturbative calculation of Maris and collaborators mrt  is also performed.

## IV Results

The results obtained for the pion form factor in the flavorless case ($`n_f=0`$) are displayed in the figure 2. A remarkable feature of the result is the significant difference between the distinct coupling constants and their respective gluon propagators. The origin of this difference is the strong dependence of the form factor in the running coupling constant. The results are also consistent with the infrared behavior of the coupling constant, for example, in the case of the Häbel propagator (Eq. (24)) when $`Q^20`$ the massive term diverges, which implies a divergent behavior. With the other propagators, the behavior is distinct because the coupling constant is finite when $`Q^20`$, giving a finite result.

As already pointed out in amn02 , the best description for the data is given by the Cornwall propagator with $`m_g=300\mathrm{MeV}`$ while the other choices for the propagator listed above do not give a reasonable result in comparison with the available experimental data.

When a more realistic number of flavors, $`n_f=3`$, is used, the result for the form factor does not present a significant change, as expected in amn02  and displayed in fig. (3), however the result with three flavors is better than with no flavors.

As the number of flavors, the change in the quark amplitude distribution, Eq. (11) and (13) does not make a significant modification in the form factor.

We display in the figures (2), (3) and (4), the full non-perturbative solution for the pion form factor obtained by Maris and collaborators mrt , only in the low momentum transfer region ($`Q^20.8\mathrm{GeV}^2`$). As pointed out in the Introduction, this solution gives a very good description of the low momentum data. Our best result, using the Cornwall’s propagator and $`m_g=0.3`$ GeV, is below of this prediction although with the same global behavior. If we diminish the smaller massive parameter $`m_g`$, we can obtain a better fit to the data (for both mesons), but this value is outside of the range of the values for this parameter, which provides good description of another processes hkpr ; ghnw .

We also display a fit with the experimental points, based on the approach of yeh02 . In that paper, the author uses the collinear expansion to analyze the NLO power corrections to the pion form factor, and founds the following expression for the pion form factor

$$F_\pi (Q^2)=\frac{16\pi \alpha _s(Q_{\mathrm{eff}}^2)f_\pi ^2}{Q^2}\left(1\frac{32\pi ^2f_\pi ^2J(Q_{\mathrm{eff}}^2)}{Q^2}\right),$$
(33)
where $`Q_{\mathrm{eff}}^2`$ is the effective energy scale and $`J(Q_{\mathrm{eff}}^2)`$ is called a jet function, introduced to absorb the infrared divergences from the NLO corrections. In the same paper, it was made a phenomenological fit to the data for the pion form factor with the following form

$$F_\pi ^{\mathrm{fit}}(Q^2)=\frac{A}{Q^2}\left(1\frac{B}{Q^2}\right),$$
(34)
In yeh02 , the best values for $`A`$ and $`B`$ to fit the above formula to the experimental data from bb76 ; amen86  are $`A=0.46895`$ and $`B=0.3009`$.

The phenomenological fit has a divergence when the momentum transferred by the photon goes to zero, in opposition to the results that employ the non-perturbative propagators and a finite infrared running coupling constant, and therefore does not give a good description for the data in low $`Q^2`$.

In order to compare with the previous results found in the literature, we plot $`Q^2F_\pi (Q^2)`$ as a function of $`Q^2`$, as shown in the fig. 4. In ja90 , its prediction does not go to zero, as in our result, neither describes the data in the region of low $`Q^2`$. In jflb01 , a comparison of predictions found in the literature for the pion form factor is made and the best result of this study (Cornwall with $`m_g=300`$ MeV) still is a good result.

In the case of the kaon, we used the same procedure, although using only the quark amplitude distribution of the Eq.(12) in the Eq.(14). The result of the flavorless case ($`n_f=0`$) is displayed in the fig.(5) while the case $`n_f=3`$ is shown in the fig.(6). As shown in the figures, the global behavior is the same as in the pion form factor, in other words, the best description for the data is given again by the Cornwall’s choices for the coupling constant and gluon propagator with $`m_g=300`$ MeV.

In the figures (5) and (6) we present a result, based in yeh02 , with the same functional form of the pion, but now applied to the available data for the kaon,

$$F_K^{\mathrm{fit}}(Q^2)=\frac{C}{Q^2}\left(1\frac{D}{Q^2}\right),$$
(35)
since the only dependence of the meson in the Eq.(33) is the normalization factor $`f_M`$. We compare Eq.(35) to the data available for the kaon iv81  and found for the coefficients $`C=0.0594695`$ and $`D=0.0130963`$ with $`\chi ^2=1.93314`$.

As the data available for the kaon restricts to the region of low $`Q^2`$, the predictions of different models cannot be tested in the kinematic region of high $`Q^2`$ for the moment, in opposition to the case of the pion.

For the pion, as well as for the kaon, when the coupling constant and the propagator of Atkinson and Bloch abw  work is employed, due to the absence of a full analytical form for the coupling constant and for the gluon propagator, we use only the infrared analytical forms, exception to the case of the fixed scale of energy (see Eq.(7a)), when we use the ultraviolet analytical result. There is a difference in the full numerical result and the analytical asymptotic expression for the coupling constant (see the first paper of abw ), but it is not significant for our result.

## V Conclusions

In this work, we calculated the pion and kaon form factor following the approach of amn02  obtaining a reasonable agreement with the data available (for the pion bb76 ; amen86 ; jflb01  and kaon iv81 ) in the case of the Cornwall’s gluon propagator jmc82  (with parameter $`m_g=300\mathrm{GeV}`$) and the running coupling constant that is frozen in the small momentum transfer region, avoiding infrared divergences. We point out that a smaller mass gives a better result, but its value is out of the mass interval found in previous results for another processes ghnw . The other propagators (Cornwall with $`m_g=700\mathrm{GeV}`$, Häbel el al. hkrsw90 , Alkofer el al. ahs , Atkinson and Bloch abw  and Gorbar and Natale gn00 ) give results in disagreement with the experimental results. The significant difference between the results has its origin in the well-known meson form factor sensibility to the running coupling constant. Since we modified the IR behavior of the $`\alpha _s`$, this difference is expected as in amn02 , where the same behavior is observed. The case with a non-zero number of flavors, $`n_f=3`$, does not have a significant difference in the results. The exchange of the quark amplitude distribution in the pion case also does not modify significantly the form factor according to ja90 . The model employed also describes the available experimental data in the low $`Q^2`$ region for both (pion, kaon) mesons form factors, in opposition to the results ja90  that employ only a frozen form of the running coupling constant, without the finite infrared form for the gluon propagator, but agrees with the results shown in jflb01  for the pion form factor. In mnp99 ; rad91 , restrictions are made concerning the applicability of the frozen coupling constant and the modified gluon propagator. In rad91  the pion form factor is calculated using the same ideas of this work, although employing a different pion wave function and a gluon propagator with fixed mass, resulting for $`Q^2F_\pi (Q^2)`$ a smaller value than the experimental one by a factor of 10. In comparison with the pure non-perturbative calculation of mrt  and the NLO perturbative calculation of yeh02 , our results for the pion (except by a normalization factor, given by adjusting the massive parameter) is an interpolation of the above results, having the fit characteristics of both calcutations. The good agreement with the experimental data of this work indicates that it is possible to include non-perturbative effects in exclusive processes in QCD without dramatic changes in the perturbative scheme of calculation of observables. Otherwise, the inclusion of non-perturbative effects is still a field in which much research is required.

###### Acknowledgements.

This research work was supported by the Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq). The authors thank A. N. Natale for the fruitful discussions during this work.


format: text

text:
# Strong interaction of correlated electrons with phonons: Exchange of phonon clouds by polarons

## I Introduction

Since the discovery of high-temperature superconductivity by Bednorz and Müller,Bednorz\_86 (1) the Hubbard model and related models such as RVB and $`t`$-$`J`$ have widely been used to discuss the physical properties of the normal and superconducting state. Dagotto\_94 (2, 3, 4, 5, 6) However, a unanimous explanation of the origin of the condensate in high-temperature superconductors has not emerged so far. One of the unsolved questions is in how far can phonons be involved in the formation of the superconducting state. In experimental and theoretical works mostly the change of phonon frequencies and phonon life times associated with the superconducting transition were discussed. For example, the decrease of frequencies of Raman-active phonons at the transition Thomson\_89 (7) and observation of the isotope effect for the case of not optimally doped superconductors Frank\_93 (8), as well as the observation of phonon-induced structure in the tunnel characteristics, Vedneev\_92 (9) speak in favor of strong electron-phonon coupling in the cuprates.

The aim of the present paper is to gain further insight into the mutual influence of strong on-site Coulomb repulsion and strong electron-phonon interaction by using the single-band Hubbard-Holstein model and a recently developed diagrammatic approach. Vladimir\_90 (10, 11, 12, 13, 14) For simplicity we consider coupling to dispersionless phonons only, although this might not be the most interesting case with respect to superconductivity. However, previous investigations Moskalenko\_97a (15, 16, 17) have shown that the Hubbard-Holstein model Hubbard\_63 (18, 19) constitutes a formidable problem of its own. Other authors have also intensively studied this model Hamiltonian. Freericks\_94 (20, 21, 22, 23)

Because the interactions between electrons and electrons and phonons are strong, we include the Coulomb repulsion in the zero-order Hamiltonian and apply the canonical transformation of Lang and Firsov Lang\_63 (24) in order to eliminate the linear electron-phonon interaction. In the strong electron-phonon coupling limit the resulting Hamiltonian of hopping polarons (i.e., hopping electrons surrounded by clouds of phonons) can lead to an attractive interaction among electrons being mediated by the phonons. In this limit the chemical potential, on-site Coulomb energy as well as the frequency of the collective mode of phonon clouds (which is much larger than the bare frequency of the Einstein oscillators) are strongly renormalized Moskalenko\_97c (17, 25, 26) affecting the dynamical properties of the polarons and the character of the superconducting transition. This will be discussed by assuming that renormalized on-site Coulomb repulsion and attractive electron-electron interaction completely cancel each other. We suggest that the resulting superconducting state with polaronic Cooper pairs is mediated by the exchange of phonon clouds during the hopping processes of the electrons.

## II Theoretical approach

### II.1 The Lang-Firsov transformation of the Hubbard-Holstein model

The initial Hamiltonian of correlated electrons coupled to optical phonons with bare frequency $`\omega _0`$ is given by

$``$ $`=`$ $`_e+_{ph}^0+_{eph},`$ (1)
$`_e`$ $`=`$ $`{\displaystyle \underset{ij\sigma }{}}\left\{t(ji)ϵ_0\delta _{ij}\right\}a_{j\sigma }^{}a_{i\sigma }+U{\displaystyle \underset{i}{}}n_i_{}n_i_{},`$ (2)
$`_{ph}^0`$ $`=`$ $`{\displaystyle \underset{i}{}}\mathrm{}\omega _0\left(b_i^{}b_i+\frac{1}{2}\right),_{eph}=g{\displaystyle \underset{i}{}}n_iq_i,`$ (3)
$`n_i`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}n_{i\sigma },n_{i\sigma }=a_{i\sigma }^{}a_{i\sigma },q_i={\displaystyle \frac{1}{\sqrt{2}}}\left(b_i+b_i^{}\right).`$ (4)
Here $`a_{i\sigma }^{}`$ ($`a_{i\sigma }`$) and $`b_i^{}`$ ($`b_i`$) are creation (annihilation) operators of electrons and phonons, respectively; $`i`$ refers to the lattice site and $`\sigma `$ to the spin; $`q_i`$ is the phonon coordinate and $`g`$ the electron-phonon interaction constant; $`ϵ_0=\overline{ϵ}_0\mu `$ with local energy $`\overline{ϵ}_0`$ and chemical potential $`\mu `$; $`U`$ the on-site Coulomb repulsion; $`t(ji)`$ is the two-center transfer integral. The Fourier representation of $`t(ji)`$ is connected to the tight-binding dispersion $`\epsilon (𝐤)`$ of the bare electrons,

$$t(ji)=\frac{1}{N}\underset{𝐤}{}\epsilon (𝐤)\mathrm{exp}\{i𝐤(𝐑_j𝐑_i)\},$$

with band width $`W`$. The energy scale of this model is fixed by the parameters $`W`$, $`U`$, $`g`$ and $`\mathrm{}\omega _0`$. As an additional parameter we have the band filling.

After applying the Lang-Firsov transformation Lang\_63 (24)

$$_p=e^Se^S,c_{i\sigma }=e^Sa_{i\sigma }e^S,c_{i\sigma }^{}=e^Sa_{i\sigma }e^S$$
(5)
with

$$S=i\overline{g}\underset{i}{}n_ip_i,\overline{g}=\frac{g}{\mathrm{}\omega _0},p_i=\frac{i}{\sqrt{2}}\left(b_i^{}b_i\right),$$
(6)
where $`p_i`$ is the phonon momentum and $`\overline{g}`$ the dimensionless interaction constant, the polaron Hamiltonian is obtained as

$`_p`$ $`=`$ $`_p^0+_{ph}^0+_{int},`$ (7)
$`_p^0`$ $`=`$ $`{\displaystyle \underset{i}{}}_{ip}^0,_{ip}^0=ϵ{\displaystyle \underset{\sigma }{}}n_{i\sigma }+\overline{U}n_in_i,`$ (8)
$`_{int}`$ $`=`$ $`{\displaystyle \underset{i,j,\sigma }{}}t(ji)c_{j\sigma }^{}c_{i\sigma },`$ (9)
$`c_{i\sigma }^{}`$ $`=`$ $`a_{i\sigma }^{}e^{i\overline{g}p_i},c_{i\sigma }=a_{i\sigma }e^{i\overline{g}p_i},`$ (10)
$`ϵ`$ $`=`$ $`\overline{ϵ}_0\overline{\mu },\overline{\mu }=\mu +\alpha \mathrm{}\omega _0,\overline{U}=U2\alpha \mathrm{}\omega _0,\alpha =\frac{1}{2}\overline{g}^2.`$ (11)
In order to derive the polaron Hamiltonian it was necessary to include the shift of the phonon coordinate $`q__i`$ of the form $`e^Sq_ie^S=q_i\overline{g}n_i`$ which is responsible for the elimination of the linear electron-phonon interaction. The polaron Hamiltonian is by nature a polaron-phonon operator, i.e., the creation operator $`c_{i\sigma }^{}`$ and destruction operator $`c_{i\sigma }`$ in $`_p`$ should be interpreted as creation and destruction operators of polarons (electrons dressed with the displacements of the ions) which couple dynamically to the momentum of the optical phonon. In zero-order approximation (omitting $`_{int}`$) polarons and phonons are localized with strongly renormalized chemical potential $`\overline{\mu }`$ and on-site Coulomb interaction $`\overline{U}`$. The operator $`_{int}`$ describes the tunneling of polarons between the lattice sites, i.e., tunneling of electrons surrounded by clouds of phonons.

### II.2 Expansion about the atomic limit

The problem is now to deal properly with the impact of electronic correlations on the polaron problem. This can be done best by using Green’s functions provided one finds a key to deal with the spin and charge degrees of freedom. For the general case when $`\overline{U}`$ is different from zero, the Coulomb interaction has to be included in the zero-order Hamiltonian. As a consequence conventional perturbation theory of quantum statistical mechanics is not an adequate tool because it relies on the expansion of the partition function about the noninteracting state (use of traditional Wick’s theorem and conventional Feynman diagrams). Similar is the situation for composite particles like polarons, $`c_{i\sigma }=a_{i\sigma }\mathrm{exp}(i\overline{g}p_i)`$, involving operators for the electron and phonon subsystem.

It was Hubbard Hubbard\_67 (27) who proposed a graphical expansion for correlated electrons about the atomic limit in powers of the hopping integrals. This diagrammatic approach was reformulated in a systematic way for the single-band Hubbard model by Slobodyan and Stasyuk, Slobodyan\_74 (28) and, independently, by ZaitsevZaitsev\_76 (29) and further developed by Izyumov. Izyumov\_90 (30) In these approaches the complicated algebraical structure of projection or Hubbard operators was used. Therefore, it appeared to be more appropriate to develop a diagrammatic technique which uses more simple creation and annihilation operators for electrons at all intermediate stages of the theory (for details see Ref. Vladimir\_90, 10, 11). In the latter approach the averages of chronological products of interactions are reduced to $`n`$-particle Matsubara Green’s functions of the atomic system. These functions can be factorized into independent local averages by using a generalization of Wick’s theorem (GWT), which takes into account strong local correlations (details are given in Ref. Vladimir\_90, 10, 11, 25). The application of GWT yields new irreducible on-site many-particle Green’s functions or Kubo cumulants. These new functions contain all local spin and charge fluctuations. The analogical linked-cluster expansion for the Hubbard model around the atomic limit was recently reformulated by Metzner. Metzner\_91 (31)

### II.3 Averages of phonon operators

We define the temperature Green’s function for the polarons in (7) in the interaction representation by

$$𝒢(𝐱,\sigma ,\tau |𝐱^{},\sigma ^{},\tau ^{})=Tc_{𝐱\sigma }(\tau )\overline{c}_{𝐱^{}\sigma ^{}}(\tau ^{})U(\beta )_0^c,$$
(12)
with

$$c_{𝐱\sigma }(\tau )=e^{^0\tau }c_{𝐱\sigma }e^{^0\tau },\overline{c}_{𝐱\sigma }(\tau )=e^{^0\tau }c_{𝐱\sigma }^{}e^{^0\tau },$$

with $`^0=_p^0+_{ph}^0`$ and evolution operator given by

$$U(\beta )=T\mathrm{exp}\left(_0^\beta 𝑑\tau H_{int}(\tau )\right),$$
(13)
and where $`𝐱,𝐱^{}`$ are the site indices and $`\tau ,\tau ^{}`$ stand for the imaginary time with $`0<\tau <\beta `$; $`T`$ is the time ordering operator and $`\beta `$ the inverse temperature. The statistical average $`\mathrm{}_0^c`$ is evaluated with respect to the zero-order density matrix of the grand-canonical ensemble of the localized polarons and phonons

$$\frac{e^{\beta ^0}}{\mathrm{Tr}e^{\beta ^0}}=\underset{i}{}\frac{e^{\beta _{ip}^0}}{\mathrm{Tr}e^{\beta _{ip}^0}}\frac{e^{\beta _{iph}^0}}{\mathrm{Tr}e^{\beta _{iph}^0}}$$
(14)
The upper index $`\mathrm{`}c^{}`$ in (13) means that only connected diagrams must be taken into account. The density matrix (14) is factorized with respect to the lattice sites. The phonon part is easily diagonalized by using the free phonon operators, $`b_i`$ and $`b_i^{}`$, while the on-site polaron Hamiltonian contains the polaron-polaron interaction which is proportional to the renormalized parameter $`\overline{U}`$, which only can be diagonalized by using Hubbard operators. Hubbard\_63 (18) At this stage no special assumption is made about the quantity $`\overline{U}`$ and its sign and we will set up the equations of motion for the dynamical quantities for this general case. However, a detailed inspection of the equations will only be undertaken for the special case $`\overline{U}=0`$.

Wick’s theorem of weak-coupling quantum field theory can be used when evaluating statistical averages of phonon operators like, for example, the propagator for the phonon cloud,

$`\mathrm{\Phi }(\tau _1|\tau _2)=\mathrm{\Phi }(\tau _1\tau _2)`$ $``$ $`T\mathrm{exp}\{i\overline{g}[p(\tau _1)p(\tau _2)]\}_0`$ (15)
$`=`$ $`\mathrm{exp}\left(\frac{1}{2}\overline{g}^2T\left[p(\tau _1)p(\tau _2)\right]^2_0\right)`$
$`=`$ $`\mathrm{exp}\left(\sigma (\beta )+\sigma (|\tau _1\tau _2|)\right),`$
$`\mathrm{\Phi }(\tau _1,\tau _2|\tau _3,\tau _4)`$ $`T`$ $`\mathrm{exp}\{i\overline{g}[p(\tau _1)+p(\tau _2)p(\tau _3)p(\tau _4)]\}_0=`$ (16)
$`=`$ $`\mathrm{exp}\left(\frac{1}{2}\overline{g}^2T[p(\tau _1)+p(\tau _2)p(\tau _3)p(\tau _4)]^2_0\right)`$
$`=`$ $`\mathrm{exp}\{\sigma (|\tau _1\tau _3|)+\sigma (|\tau _1\tau _4|)+\sigma (|\tau _2\tau _3|)`$
$`+`$ $`\sigma (|\tau _2\tau _4|)\sigma (|\tau _1\tau _2|)\sigma (|\tau _3\tau _4|)`$
$``$ $`2\sigma (\beta )\},`$
where

$`\sigma \left(|\tau _1\tau _2|\right)`$ $`=`$ $`\overline{g}^2Tp(\tau _1)p(\tau _2)_0`$ (17)
$`=`$ $`\alpha {\displaystyle \frac{\mathrm{cosh}\left(\mathrm{}\omega _0\left\{{\displaystyle \frac{\beta }{2}}|\tau _1\tau _2|\right\}\right)}{\mathrm{sinh}\left({\displaystyle \frac{\beta \mathrm{}\omega _0}{2}}\right)}}.`$
We now have to discuss the problem of how to calculate chronological averages of combinations of polaron operators. Here we will make use of the above mentioned new diagram technique and the GWT.Vladimir\_90 (10, 11) This approach has many-particle on-site irreducible Green’s functions as main element of the diagrams.

## III Polaron and phonon Green’s functions

In zero-order approximation the one-polaron Green’s function is of the form

$`𝒢_p^0(x|x^{})`$ $`=`$ $`Tc_{𝐱\sigma }(\tau )\overline{c}_{𝐱^{}\sigma ^{}}(\tau ^{})_0`$ (18)
$`=`$ $`Ta_{𝐱\sigma }(\tau )\overline{a}_{𝐱^{}\sigma ^{}}(\tau ^{})_0\mathrm{\Phi }(\tau |\tau ^{})`$
$`=`$ $`𝒢^{(0)}(x|x^{})\mathrm{\Phi }(\tau |\tau ^{}),`$
with $`x=(𝐱,\sigma ,\tau )`$. The simplest new element of the diagram technique is the two-particle irreducible Green’s function or Kubo cumulant, which is equal to

$`𝒢_2^{(0)ir}(x_1,x_2|x_3,x_4)=`$ $`\delta _{𝐱_1,𝐱_2}\delta _{𝐱_1,𝐱_3}\delta _{𝐱_1,𝐱_4}`$
$`\times 𝒢_2^{(0)ir}(\sigma _1,\tau _1;\sigma _2,\tau _2|\sigma _3,\tau _3;\sigma _4,\tau _4),`$ (19)
where

$`𝒢_2^{(0)ir}`$ $`(\sigma _1,\tau _1;\sigma _2,\tau _2|\sigma _3,\tau _3;\sigma _4,\tau _4)`$
$`=`$ $`Tc_{\sigma _1}(\tau _1)c_{\sigma _2}(\tau _2)\overline{c}_{\sigma _3}(\tau _3)\overline{c}_{\sigma _4}(\tau _4)_0`$
$`Tc_{\sigma _1}(\tau _1)\overline{c}_{\sigma _4}(\tau _4)_0Tc_{\sigma _2}(\tau _2)\overline{c}_{\sigma _3}(\tau _3)_0`$
$`+Tc_{\sigma _1}(\tau _1)\overline{c}_{\sigma _3}(\tau _3)_0Tc_{\sigma _2}(\tau _2)\overline{c}_{\sigma _4}(\tau _4)_0.`$ (20)
The first term on the right-hand side of Eq. (21) is of the form

$``$ $`Tc_{\sigma _1}(\tau _1)c_{\sigma _2}(\tau _2)\overline{c}_{\sigma _3}(\tau _3)\overline{c}_{\sigma _4}(\tau _4)_0`$
$`=Ta_{\sigma _1}(\tau _1)a_{\sigma _2}(\tau _2)\overline{a}_{\sigma _3}(\tau _3)\overline{a}_{\sigma _4}(\tau _4)_0\mathrm{\Phi }(\tau _1\tau _2|\tau _3\tau _4)`$ (21)
When the number of polaron operators increases more complicated irreducible Green’s Functions like $`𝒢_n^{(0)ir}(x_1\mathrm{}x_n|x_1^{}\mathrm{}x_n^{})`$ with $`n3`$ and all possible terms of their products will appear. The sum of all strongly connected diagrams (i.e., those which cannot be divided into two parts by cutting a single hopping line) containing all kinds of irreducible Green’s functions in the perturbation expansion of the evolution operator, defines the special function $`Z(x|x^{})`$ (for details see Ref. Vladimir\_90, 10, 11). This function contains all contributions from charge and spin fluctuations. It allows us, together with the mass operator which is in our case the hopping matrix element, to formulate a Dyson-type of equation for the one-polaron Green’s function, Vladimir\_90 (10, 11, 12, 13, 14)

$$𝒢(x|x^{})=\mathrm{\Lambda }(x|x^{})+\underset{1,2}{}\mathrm{\Lambda }(x|1)t(12)𝒢(2|x^{}),$$
(22)
where

$`\mathrm{\Lambda }(x|x^{})`$ $`=`$ $`𝒢_p^{(0)}(x|x^{})+Z(x|x^{}),`$ (23)
$`t(xx^{})`$ $`=`$ $`\delta _{\sigma ,\sigma ^{}}\delta (\tau \tau ^{})t(𝐱𝐱^{}).`$ (24)
Here $`x`$ stands again for $`𝐱,\sigma ,\tau `$ and the sum is over the discrete indices including the integration over $`\tau `$. Using the Fourier representation for these quantities,

$`𝒢(𝐱|\tau )`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\omega _n}{}}e^{i\mathrm{𝐤𝐱}i\omega _n\tau }𝒢_\sigma (𝐤|i\omega _n),`$
$`\mathrm{\Lambda }(𝐱|\tau )`$ $`=`$ $`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\omega _n}{}}e^{i\mathrm{𝐤𝐱}i\omega _n\tau }\mathrm{\Lambda }_\sigma (𝐤|i\omega _n),`$
$`𝒢(𝐱|i\omega _n)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _\beta ^{+\beta }}𝑑\tau e^{i\omega _n\tau }𝒢_\sigma (𝐱|\tau ),`$
$`\mathrm{\Lambda }(𝐱|i\omega _n)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _\beta ^{+\beta }}𝑑\tau e^{i\omega _n\tau }\mathrm{\Lambda }_\sigma (𝐱|\tau ),`$ (25)
gives us the following form of the Dyson equation for the renormalized one-polaron Green’s function,

$$𝒢_\sigma (𝐤|i\omega _n)=\frac{\mathrm{\Lambda }_\sigma (𝐤)}{1\epsilon (𝐤)\mathrm{\Lambda }(𝐤|i\omega _n)}$$
(26)
Here $`\omega _n`$ stands for the odd Matsubara frequency, $`\omega _n=(2n+1)\pi /\beta `$.

For further discussion of $`𝒢_\sigma (k|i\omega )`$ we need the Fourier representation of the zero-order one-polaron Green’s function $`𝒢_p^{(0)}`$ defined in (18). In order to facilitate the investigation we have evaluated the propagator of the phonon cloud (16) in the strong-coupling limit $`\alpha 1`$, Moskalenko\_97a (15, 16, 26)

$`\mathrm{\Phi }(\tau )`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\mathrm{\Omega }_n}{}}e^{i\mathrm{\Omega }_n(\tau )}\overline{\mathrm{\Phi }}(i\mathrm{\Omega }_n),`$ (27)
$`\overline{\mathrm{\Phi }}(i\mathrm{\Omega }_n)`$ $`=`$ $`{\displaystyle \frac{e^{\sigma (\beta )}}{2}}{\displaystyle _\beta ^\beta }𝑑\tau e^{i\mathrm{\Omega }_n\tau +\sigma (|\tau |)}`$ (28)
with $`\mathrm{\Omega }_n=2n\pi /\beta `$. In order to find $`\overline{\mathrm{\Phi }}(i\mathrm{\Omega }_n)`$ we use the Laplace approximation Fedoryuk\_87 (32) for the integral (28) which contains an exponential function with the parameter $`\alpha `$. In the strong-coupling limit $`\alpha 1`$ we have obtained

$`\overline{\mathrm{\Phi }}(i\mathrm{\Omega }n){\displaystyle \frac{2\omega _c}{\mathrm{\Omega }_n^2+\omega _c^2}},\omega _c=\mathrm{}\alpha \omega _0=\mathrm{}{\displaystyle \frac{g^2}{2\mathrm{}\omega _0}}.`$ (29)
This term is the harmonic propagator of the collective mode of phonons belonging to the polaron clouds. There are further terms describing anharmonic deviations. For $`\alpha 1`$ these terms will be omitted since they are smaller compared to the harmonic contribution. We obtain then for the Fourier representation of the phonon correlation function by making use of the Laplace approximation Fedoryuk\_87 (32) and

$`\mathrm{\Phi }(\tau _1,\tau _2|\tau _3,\tau _4)`$ $`={\displaystyle \frac{1}{\beta ^4}}{\displaystyle \underset{\mathrm{\Omega }_1\mathrm{}\mathrm{\Omega }_4}{}}\overline{\mathrm{\Phi }}(i\mathrm{\Omega }_1,i\mathrm{\Omega }_2|i\mathrm{\Omega }_3,i\mathrm{\Omega }_4)e^{i\mathrm{\Omega }_1\tau _1i\mathrm{\Omega }_2\tau _2+i\mathrm{\Omega }_3\tau _3+i\mathrm{\Omega }_4\tau _4},`$ (30)
$`\overline{\mathrm{\Phi }}(i\mathrm{\Omega }_1,i\mathrm{\Omega }_2|i\mathrm{\Omega }_3,i\mathrm{\Omega }_4)`$ $`={\displaystyle _0^\beta }\mathrm{}{\displaystyle _0^\beta }𝑑\tau _1\mathrm{}𝑑\tau _4e^{i\mathrm{\Omega }_1\tau _1+i\mathrm{\Omega }_2\tau _2i\mathrm{\Omega }_3\tau _3i\mathrm{\Omega }_4\tau _4},`$ (31)
the result
$`\overline{\mathrm{\Phi }}(i\mathrm{\Omega }_1,i\mathrm{\Omega }_2|i\mathrm{\Omega }_3,i\mathrm{\Omega }_4)`$ $`\left[\delta _{\mathrm{\Omega }_1,\mathrm{\Omega }_3}\delta _{\mathrm{\Omega }_2,\mathrm{\Omega }_4}+\delta _{\mathrm{\Omega }_1,\mathrm{\Omega }_4}\delta _{\mathrm{\Omega }_2,\mathrm{\Omega }_3}\right]\overline{\mathrm{\Phi }}(i\mathrm{\Omega }_1)\overline{\mathrm{\Phi }}(i\mathrm{\Omega }_2),`$ (32)
which corresponds to
$`\mathrm{\Phi }(\tau _1,\tau _2|\tau _3,\tau _4)`$ $`\mathrm{\Phi }(\tau _1|\tau _3)\mathrm{\Phi }(\tau _2|\tau _4)+\mathrm{\Phi }(\tau _1|\tau _4)\mathrm{\Phi }(\tau _2|\tau _3).`$ (33)
This means that we keep in the following only the free collective oscillations of the phonon clouds (29) which surround the polarons and use the Hartree-Fock approximation (32) and (33) for their two-particle correlation functions. In particular we will investigate the influence of absorption and emission of this collective mode (by the polarons) on the superconducting phase transition.

By using the harmonic mode (29) the Fourier representation of the local polaron Green’s function

$$\overline{𝒢}_{p\sigma }^{(0)}(i\omega _n)=\frac{1}{2}_\beta ^\beta 𝑑\tau e^{i\omega _n\tau }\overline{𝒢}_{p\sigma }^{(0)}(\tau )$$
(34)
is equal to

$`\overline{𝒢}_{p\sigma }^{(0)}(i\omega _n){\displaystyle \frac{1}{Z_0}}`$ $`({\displaystyle \frac{e^{\beta E_0}+\overline{N}(\omega _c)(e^{\beta E_0}+e^{\beta E_\sigma })}{i\omega _n+E_0E_\sigma \omega _c}}`$
$`+{\displaystyle \frac{e^{\beta E_\sigma }+\overline{N}(\omega _c)(e^{\beta E_0}+e^{\beta E_\sigma })}{i\omega _n+E_0E_\sigma +\omega _c}}`$
$`+{\displaystyle \frac{e^{\beta E_\sigma }+\overline{N}(\omega _c)(e^{\beta E_\sigma }+e^{\beta E_2})}{i\omega _n+E_\sigma E_2\omega _c}}`$
$`+{\displaystyle \frac{e^{\beta E_2}+\overline{N}(\omega _c)(e^{\beta E_\sigma }+e^{\beta E_2})}{i\omega _n+E_\sigma E_2+\omega _c}}),`$ (35)
where

$`Z_0`$ $`=1+e^{\beta E_\sigma }+e^{\beta E_\sigma }+e^{\beta E_2},`$ (36a)
$`E_0`$ $`=0,E_{\pm \sigma }=ϵ,E_2=\overline{U}+2ϵ,`$ (36b)
$$\overline{n}(ϵ)=\left(e^{\beta ϵ}+1\right)^1,\overline{N}(\omega _c)=\left(e^{\beta \omega _c}1\right)^1.$$
(37)
Equation (35) shows that the on-site transition energies of the polarons are changed by the collective-mode energy $`\pm \omega _c`$ of the phonon clouds. The delocalization of the polarons due to their hopping between the lattice sites causes the broadening of the polaron energy levels. Equation (35) can be further simplified for the case of small on-site interaction energy $`\overline{U}`$ of polarons. For $`\overline{U}=0`$ we obtain

$`\overline{𝒢}_{p\sigma }^{(0)}(i\omega _n|ϵ)`$ $`={\displaystyle \frac{\overline{N}(\omega _c)+1\overline{n}(ϵ)}{i\omega _nϵ\omega _c}}+{\displaystyle \frac{\overline{N}(\omega _c)+\overline{n}(ϵ)}{i\omega _nϵ+\omega _c}}`$
$`={\displaystyle \frac{(i\omega _nϵ)\mathrm{coth}(\beta \omega _c/2)+\omega _c\mathrm{tanh}(\beta ϵ/2)}{(i\omega _nϵ)^2\omega _c^2}}.`$ (38)
This function has the antisymmetric property

$$\overline{𝒢}_{p\sigma }^{(0)}(i\omega _n|ϵ)=\overline{𝒢}_{p\sigma }^{(0)}(i\omega _n|ϵ)$$
(39)
which holds also for the renormalized polaron quantities,

$`\mathrm{\Lambda }_\sigma (𝐤,i\omega _n|ϵ)`$ $`=\mathrm{\Lambda }_\sigma (𝐤,i\omega _n|ϵ),`$
$`𝒢_\sigma (𝐤,i\omega _n|ϵ)`$ $`=𝒢_\sigma (𝐤,i\omega _n|ϵ).`$ (40)
When taking $`\overline{U}0`$ we assume that the strong on-site Coulomb repulsion of polarons can be canceled by the attraction induced by the strong electron-phonon interaction. We consider this as a model case which allows to discuss in a transparent manner the polarons exchanging phonon clouds during hopping between the lattice sites.

## IV Two-particle irreducible correlation functions

In the following we discuss the influence of strong electron-phonon interaction on the two-particle irreducible Green’s function. For $`\overline{U}=0`$ the electronic correlation function in (22) is equal to

$`T`$ $`a_{\sigma _1}(\tau _1)a_{\sigma _2}(\tau _2)\overline{a}_{\sigma _3}(\tau _3)\overline{a}_{\sigma _4}(\tau _4)_0`$
$`=`$ $`Ta_{\sigma _1}(\tau _1)\overline{a}_{\sigma _4}(\tau _4)_0Ta_{\sigma _2}(\tau _2)\overline{a}_{\sigma _3}(\tau _3)_0`$
$`Ta_{\sigma _1}(\tau _1)\overline{a}_{\sigma _3}(\tau _3)_0Ta_{\sigma _2}(\tau _2)\overline{a}_{\sigma _4}(\tau _4)_0`$ (41)
because usual Wick’s theorem works now. Making use of (33) we obtain for the two-particle irreducible Green’s function (21) the relation,

$`𝒢_2^{(0)ir}`$ $`(\sigma _1,\tau _1;\sigma _2,\tau _2|\sigma _3,\tau _3;\sigma _4\tau _4)`$
$`=`$ $`\delta _{\sigma _1,\sigma _4}\delta _{\sigma _2,\sigma _3}𝒢_{\sigma _1}^{(0)}(\tau _1\tau _4)𝒢_{\sigma _2}^{(0)}(\tau _2\tau _3)\mathrm{\Phi }(\tau _1\tau _3)\mathrm{\Phi }(\tau _2\tau _4)`$
$`\delta _{\sigma _1,\sigma _3}\delta _{\sigma _2,\sigma _4}𝒢_{\sigma _1}^{(0)}(\tau _1\tau _3)𝒢_{\sigma _2}^{(0)}(\tau _2\tau _4)\mathrm{\Phi }(\tau _1\tau _4)\mathrm{\Phi }(\tau _2\tau _3).`$ (42)
In the absence of exchange of phonon clouds by the polarons this quantity has to be zero. Indeed, if during the time of propagation of two polarons the electrons keep their initial phonon clouds, then the irreducible two-polaron Green’s function (21) will vanish for the case $`\overline{U}=0`$. But because two electrons can be exchanged (independently of the exchange of phonon clouds) we obtain new contributions corresponding to two polarons with exchanged phonon clouds. Alternatively we may say that, in the case $`\overline{U}=0`$, Wick’s theorem works separately for free electrons and free phonons; however, it does not apply for polarons as composite particles. Hence their cumulants do not vanish.

The Fourier representation of (42),

$`𝒢_2^{(0)ir}`$ $`(\sigma _1,i\omega _1;\sigma _2,i\omega _2|\sigma _3,i\omega _3;\sigma _4,i\omega _4)`$
$`=`$ $`{\displaystyle _0^\beta }\mathrm{}{\displaystyle _0^\beta }𝑑\tau _1\mathrm{}𝑑\tau _4𝒢_2^{(0)ir}(\sigma _1,\tau _1;\sigma _2,\tau _2|\sigma _3,\tau _3;\sigma _4,\tau _4)e^{i\omega _1\tau _1+i\omega _2\tau _2i\omega _3\tau _3i\omega _4\tau _4},`$ (43)
has the form

$`𝒢_2^{(0)ir}`$ $`(\sigma _1,i\omega _1;\sigma _2,i\omega _2|\sigma _3,i\omega _3;\sigma _4,i\omega _4)`$
$`=`$ $`\beta \delta _{\omega _1+\omega _2,\omega _3+\omega _4}\overline{𝒢}_2^{(0)ir}(\sigma _1,i\omega _1;\sigma _2,i\omega _2|\sigma _3,i\omega _3;\sigma _4,i\omega _4)`$
$`=`$ $`\beta \delta _{\omega _1+\omega _2,\omega _3+\omega _4}\{\delta _{\sigma _1,\sigma _4}\delta _{\sigma _2,\sigma _3}A_{\sigma _1\sigma _2}(\sigma _1,i\omega _1;\sigma _2,i\omega _2|\sigma _2,i\omega _3;\sigma _1,i\omega _4)`$
$`\delta _{\sigma _1,\sigma _3}\delta _{\sigma _2,\sigma _4}A_{\sigma _1\sigma _2}(\sigma _1,i\omega _1;\sigma _2,i\omega _2|\sigma _1,i\omega _4;\sigma _2,i\omega _3)\},`$ (44)
where

$`A_{\sigma _1\sigma _2}`$ $`(\sigma _1,i\omega _1;\sigma _2,i\omega _2|\sigma _2,i\omega _3;\sigma _1,i\omega _4)`$
$`={\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\mathrm{\Omega }}{}}{\displaystyle \frac{(2\omega _c)^2}{[i\omega _1\mathrm{\Omega }ϵ][i\omega _3\mathrm{\Omega }ϵ][\mathrm{\Omega }^2+\omega _c^2][(\mathrm{\Omega }+\mathrm{\Omega }_1)^2+\omega _c^2]}}`$ (45)
with $`\mathrm{\Omega }_1=\omega _2\omega _3`$. The summation leads to

$`A_{\sigma _1\sigma _2}`$ $`(\sigma _1,i\omega _1;\sigma _2,i\omega _2|\sigma _2,i\omega _3;\sigma _1,i\omega _4)=2(\omega _c)^2`$
$`\times \{{\displaystyle \frac{\mathrm{tanh}(\beta ϵ/2)[i\omega _1+i\omega _22ϵ][2\omega _c^2(i\omega _1ϵ)(i\omega _4ϵ)(i\omega _2ϵ)(i\omega _3ϵ)]}{[(i\omega _1ϵ)^2\omega _c^2][(i\omega _2ϵ)^2\omega _c^2][(i\omega _3ϵ)^2\omega _c^2][(i\omega _4ϵ)^2\omega _c^2]}}`$
$`{\displaystyle \frac{2\omega _c[i\omega _1+i\omega _22ϵ]^2\mathrm{coth}(\beta \omega _c/2)[2\omega _c^2(i\omega _1ϵ)(i\omega _3ϵ)(i\omega _2ϵ)(i\omega _4ϵ)]}{[(i\mathrm{\Omega }_1)^2(2\omega _c)^2][(i\omega _1ϵ)^2\omega _c^2][(i\omega _2ϵ)^2\omega _c^2][(i\omega _3ϵ)^2\omega _c^2][(i\omega _4ϵ)^2\omega _c^2]}}`$
$`{\displaystyle \frac{\mathrm{coth}(\beta \omega _c/2)}{\omega _c[(i\mathrm{\Omega }_1)^2(2\omega _c)^2]}}[{\displaystyle \frac{(i\omega _1ϵ)(i\omega _3ϵ)+3\omega _c^2}{[(i\omega _1ϵ)^2\omega _c^2][(i\omega _3ϵ)^2\omega _c^2]}}`$
$`+{\displaystyle \frac{(i\omega _2ϵ)(i\omega _4ϵ)+3\omega _c^2}{[(i\omega _2ϵ)^2\omega _c^2][(i\omega _4ϵ)^2\omega _c^2]}}]\}`$ (46)
The function $`A_{\sigma _1\sigma _2}`$ contains the contributions to the two-particles on-site Green’s function from the different spin channels. The spin structure in Eq. (44) is due to the conversation law for the spins of the polarons.

## V Superconducting phase transition

In the following we check whether the polaronic system may have a superconducting instability in the absence of a direct attractive interaction for the polarons, i.e., for $`\overline{U}=0`$. In this case the attraction is only brought about dynamically by polarons exchanging phonon clouds. With respect to superconductivity we need in addition to the normal state Green’s function (13) the anomalous propagators. Abrikosov\_62 (33) For simplicity we limit the discussion to $`s`$-wave superconductivity as in previous investigations of superconducting instabilities in the Hubbard model Bogoliubov\_91 (13, 14) and Hubbard-Holstein model in the strong-coupling limit, $`\alpha 1`$. Moskalenko\_99b (26)

For a description of the superconducting state we need the three irreducible functions $`\mathrm{\Lambda }_\sigma `$, $`Y_{\sigma ,\sigma }`$ and $`\overline{Y}_{\sigma ,\sigma }`$ which represent infinite sums of diagrams containing irreducible many-particle Green’s functions. In order to obtain a close set of equations we will restrict ourselves to a class of rather simple contributions which, however, contain the most important charge, spin and pairing correlations; for details see Ref. Moskalenko\_99b, 26. This class of diagrams is obtained by neglecting contributions for which the Fourier representation of the superconducting order parameters, $`Y_{\sigma ,\sigma }`$ and $`\overline{Y}_{\sigma ,\sigma }`$, depend on the polaron momentum $`𝐤`$. In this approximation $`Y_{\sigma ,\sigma }`$ has to be obtained from

$`Y_{\sigma ,\sigma }(i\omega )=`$ $`{\displaystyle \frac{1}{\beta N}}{\displaystyle \underset{𝐤,\omega _l}{}}{\displaystyle \frac{\epsilon (𝐤)\epsilon (𝐤)Y_{\sigma ,\sigma }(i\omega _l)}{D_\sigma (𝐤,i\omega _l)}}`$
$`\times \overline{𝒢}_2^{(0)ir}(\sigma ,i\omega ;\sigma ,i\omega |\sigma ,i\omega _l;\sigma ,i\omega _l).`$ (47)
In the same approximation $`\mathrm{\Lambda }_\sigma `$ has to be computed from

$`\mathrm{\Lambda }_\sigma (i\omega )=`$ $`𝒢_{p\sigma }^{(0)}(i\omega )`$
$`{\displaystyle \frac{1}{\beta N}}{\displaystyle \underset{𝐤,\omega _l}{}}{\displaystyle \frac{\epsilon ^2(𝐤)}{D_\sigma (𝐤,i\omega _l)}}\{\mathrm{\Lambda }_\sigma (i\omega _l)[1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega _l)]`$
$`\epsilon (𝐤)Y_{\sigma ,\sigma }(i\omega _l)\overline{Y}_{\sigma ,\sigma }(i\omega _l)\}\overline{𝒢}_2^{(0)ir}(\sigma ,i\omega ;\sigma ,i\omega _l|\sigma ,i\omega _l;\sigma ,i\omega )`$
$`{\displaystyle \frac{1}{\beta N}}{\displaystyle \underset{𝐤,\omega _l}{}}{\displaystyle \frac{\epsilon ^2(𝐤)}{D_\sigma (𝐤,i\omega _l)}}\{\mathrm{\Lambda }_\sigma (i\omega _l)[1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega _l)]`$
$`\epsilon (𝐤)Y_{\sigma ,\sigma }(i\omega _l)\overline{Y}_{\sigma ,\sigma }(i\omega _l)\}\overline{𝒢}_2^{(0)ir}(\sigma ,i\omega ;\sigma ,i\omega _l|\sigma ,i\omega _l;\sigma ,i\omega )`$ (48)
with

$$D_\sigma (𝐤,i\omega )=[1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega )][1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega )]+\epsilon (𝐤)\epsilon (𝐤)Y_{\sigma ,\sigma }(i\omega )\overline{Y}_{\sigma ,\sigma }(i\omega ).$$
(49)
The corresponding equation for $`\overline{Y}_{\sigma ,\sigma }(i\omega )`$ can be obtained from the expression for $`Y_{\sigma ,\sigma }(i\omega )`$.

The last equations together with the equations for the one- and two-particle Green’s functions determine completely the properties of the superconducting phase, provided it exists. In order to gain further insight into the physics contained in (47) and (48) we will linearize the equations in terms of the order parameter $`Y_{\sigma ,\sigma }(i\omega )`$ which will determine the critical temperature $`T_c`$. The resulting equation for the order parameter is of the form

$`Y_{\sigma ,\sigma }(i\omega )=`$ $`{\displaystyle \frac{1}{\beta N}}{\displaystyle \underset{𝐤,\omega _l}{}}{\displaystyle \frac{\epsilon (𝐤)\epsilon (𝐤)Y_{\sigma ,\sigma }(i\omega _l)}{[1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega _l)][1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega _l)]}}`$
$`\times \overline{𝒢}_2^{(0)ir}(\sigma ,i\omega ;\sigma ,i\omega |\sigma ,i\omega _l;\sigma ,i\omega _l).`$ (50)
This equation must be solved together with the equation for $`\mathrm{\Lambda }_\sigma (i\omega )`$ which may be approximated by setting the order parameters to zero giving,

$`\mathrm{\Lambda }_\sigma (i\omega )=𝒢_{p\sigma }^{(0)}(i\omega )`$ $`{\displaystyle \frac{1}{\beta N}}{\displaystyle \underset{𝐤,\omega _l}{}}{\displaystyle \frac{\epsilon ^2(𝐤)\mathrm{\Lambda }_\sigma (i\omega _l)}{1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega _l)}}\overline{𝒢}_2^{(0)ir}(\sigma ,i\omega ;\sigma ,i\omega _l|\sigma ,i\omega _l;\sigma ,i\omega )`$
$`{\displaystyle \frac{1}{\beta N}}{\displaystyle \underset{𝐤,\omega _l}{}}{\displaystyle \frac{\epsilon ^2(𝐤)\mathrm{\Lambda }_\sigma (i\omega _l)}{1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega _l)}}\overline{𝒢}_2^{(0)ir}(\sigma ,i\omega ;\sigma ,i\omega _l|\sigma ,i\omega _l;\sigma ,i\omega ).`$ (51)
In order to determine $`T_c`$ we must solve (51) for $`\mathrm{\Lambda }_\sigma `$ and insert the result in (50). The irreducible functions in (50) and (51) can be written as

$`\overline{𝒢}_2^{(0)ir}`$ $`(\sigma ,i\omega ;\sigma ,i\omega _l|\sigma ,i\omega _l;\sigma ,i\omega )={\displaystyle \frac{(\omega \omega _l)^2}{\mathrm{\Delta }^2\mathrm{\Delta }_l^2}}\{2\omega _c^2(x+x_l)\mathrm{tanh}(\beta ϵ/2)`$
$`{\displaystyle \frac{\mathrm{coth}(\beta \omega _c/2)}{\omega _c[(i\omega i\omega _l)^24\omega _c^2]}}[(xx_l+\omega _c^2)(\mathrm{\Delta }\mathrm{\Delta }_l+8\omega _c^4)2\omega _c^2(\mathrm{\Delta }+\mathrm{\Delta }_l)(xx_l\omega _c^2)]\},`$ (52)
$`\overline{𝒢}_2^{(0)ir}`$ $`(\sigma ,i\omega ;\sigma ,i\omega _l|\sigma ,i\omega _l;\sigma ,i\omega )={\displaystyle \frac{2\omega _c}{\mathrm{\Delta }^2\mathrm{\Delta }_l^2}}\{\omega _c(x+x_l)(\mathrm{\Delta }+\mathrm{\Delta }_l)\mathrm{tanh}(\beta ϵ/2)`$
$`+\mathrm{coth}(\beta \omega _c/2)(x+x_l)^2(xx_l\omega _c^2)\}+{\displaystyle \frac{\mathrm{coth}(\beta \omega _c/2)((xx_l+3\omega _c^2)}{\omega _c\mathrm{\Delta }\mathrm{\Delta }_l}},`$ (53)
$`\overline{𝒢}_2^{(0)ir}`$ $`(\sigma ,i\omega ;\sigma ,i\omega |\sigma ,i\omega _l;\sigma ,i\omega _l)=`$
$`{\displaystyle \frac{2ϵ(2\omega _c)^2\mathrm{tanh}(\beta ϵ/2)[i\omega i\omega _l+ϵ^2\omega _c^2]}{[\omega ^2+(ϵ+\omega _c)^2][\omega ^2+(ϵ\omega _c)^2][\omega _l^2+(ϵ+\omega _c)^2][\omega _l^2+(ϵ\omega _c)^2]}}`$
$`+{\displaystyle \frac{2\omega _c\mathrm{coth}(\beta \omega _c/2)[i\omega i\omega _l+2\omega _c(ϵ\omega _c)(ϵ\omega _c)^2]}{[\omega ^2+(ϵ\omega _c)^2][\omega _l^2+(ϵ\omega _c)^2][(\omega \omega _l)^2+(ϵ+(2\omega _c)^2]}}`$
$`+{\displaystyle \frac{2\omega _c\mathrm{coth}(\beta \omega _c/2)[i\omega i\omega _l2\omega _c(ϵ+\omega _c)(ϵ+\omega _c)^2]}{[\omega ^2+(ϵ+\omega _c)^2][\omega _l^2+(ϵ+\omega _c)^2][(\omega \omega _l)^2+(ϵ+(2\omega _c)^2]}},`$ (54)
where

$`x`$ $`=i\omega ϵ,`$ $`\mathrm{\Delta }`$ $`=(i\omega ϵ)^2\omega _c^2`$ (55a)
$`x_l`$ $`=i\omega _lϵ,`$ $`\mathrm{\Delta }_l`$ $`=(i\omega _lϵ)^2\omega _c^2.`$ (55b)
For further discussion of (50) and (51) we introduce the following shorthand notation:

$`\varphi _\sigma (i\omega )`$ $`={\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\epsilon ^2(𝐤)\mathrm{\Lambda }_\sigma (i\omega )}{1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega )}}`$
$`={\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\epsilon (𝐤)}{1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega )}},`$ (56)
$`g_\sigma (i\omega )`$ $`=𝒢(𝐱=𝐱^{}|i\omega )`$
$`={\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\mathrm{\Lambda }_\sigma (i\omega )}{1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega )}},`$ (57)
$`\varphi _\sigma ^{sc}(i\omega )`$ $`={\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{\epsilon (𝐤)\epsilon (𝐤)}{[1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega )][1\epsilon (𝐤)\mathrm{\Lambda }_\sigma (i\omega )]}}`$
$`={\displaystyle \frac{\varphi _\sigma (i\omega )\varphi _\sigma (i\omega )}{\mathrm{\Lambda }_\sigma (i\omega )\mathrm{\Lambda }_\sigma (i\omega )}}.`$ (58)
Furthermore we assume that $`\epsilon (𝐤)=\epsilon (𝐤)`$ holds with $`_𝐤\epsilon (𝐤)=_𝐤\epsilon ^3(𝐤)=0`$. Sums will be replaced by integrals,

$`{\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐤}{}}`$ $`={\displaystyle 𝑑\epsilon \rho _0(\epsilon )},`$ (59)
$`\rho _0(\epsilon )`$ $`={\displaystyle \frac{4}{\pi W}}\sqrt{1\left({\displaystyle \frac{2\epsilon }{W}}\right)^2}\times \{\begin{array}{cc}1\hfill & |\epsilon |<\frac{W}{2},\hfill \\ 0\hfill & |\epsilon |<\frac{W}{2},\hfill \end{array}`$ (60)
where $`W`$ is the band width and $`\rho _0`$ a model density of states of semielliptic form. Since we do not consider magnetic states here, the spin index in the paramagnetic phase can be omitted:

$`\mathrm{\Lambda }_\sigma (i\omega )`$ $`=\mathrm{\Lambda }_\sigma (i\omega )=\mathrm{\Lambda }(i\omega ),`$ (61a)
$`\varphi _\sigma (i\omega )`$ $`=\varphi _\sigma (i\omega )=\varphi (i\omega ),`$ (61b)
$`\varphi _\sigma ^{sc}(i\omega )`$ $`=\varphi _\sigma ^{sc}(i\omega )=\varphi ^{sc}(i\omega ).`$ (61c)
However, the spin index is essential for the superconducting order parameter $`Y_{\sigma ,\sigma }(i\omega )`$,

$`Y_{\sigma ,\sigma }(i\omega )`$ $`=g_{\sigma ,\sigma }Y(i\omega ),`$ (62a)
$`g_{\sigma ,\sigma }`$ $`=\delta _{\sigma ,}\delta _{\sigma ,},`$ (62b)
where $`Y(i\omega )`$ is an even function of the frequency,

$$Y(i\omega )=Y(i\omega ).$$
(63)
Finally we have to add the equation which determines the chemical potential:

$$\frac{1}{\beta }\underset{\omega _n}{}\underset{\sigma }{}𝒢_\sigma (𝐱=𝐱^{}|i\omega )e^{i\omega _n0^+}=\frac{2}{\beta }\underset{\omega _n}{}g_\sigma (i\omega )e^{i\omega _n0^+}=\frac{N_p}{N}.$$
(64)
Here $`N_p`$ is the number of polarons and $`N`$ the number of lattice sites. With (59) and (60) the functions (56) and (57) can be written as

$`\varphi (i\omega )`$ $`={\displaystyle \frac{W}{2}}{\displaystyle \frac{\left(1\sqrt{1\lambda ^2(i\omega )}\right)^2}{\lambda ^3(i\omega )}}={\displaystyle \frac{W}{2}}{\displaystyle \frac{\lambda (i\omega )}{\left(1+\sqrt{1\lambda ^2(i\omega )}\right)^2}}`$ (65)
$`g(i\omega )`$ $`={\displaystyle \frac{4}{W}}{\displaystyle \frac{\left(1\sqrt{1\lambda ^2(i\omega )}\right)^2}{\lambda (i\omega )}}={\displaystyle \frac{4}{W}}{\displaystyle \frac{\lambda (i\omega )}{\left(1+\sqrt{1\lambda ^2(i\omega )}\right)^2}},`$ (66)
with $`\lambda (i\omega )=(W/2)\mathrm{\Lambda }(i\omega )`$. In order to check whether the state to be determined from Eq. (51) is metallic or dielectric one has to analyze the renormalized density of states given by

$$\rho (E)=\frac{1}{\pi }\mathrm{Im}g(E+i0^+)=\frac{1}{\pi }\mathrm{Im}\left(\frac{1\sqrt{1\lambda ^2(E+i0^+)}}{\lambda (E+i0^+)}\right),$$
(67)
where $`\lambda (E+i0^+)`$ is the analytical continuation of $`\lambda (i\omega )`$.

## VI Analytical solutions

The expressions for $`\mathrm{\Lambda }(i\omega )`$ and $`Y(i\omega )`$ simplify by using the shorthand notations (56) and (57) and the symmetry property (62):

$`Y(i\omega )`$ $`={\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\omega _l}{}}\varphi ^{sc}(i\omega _l)\overline{𝒢}^{(0)ir}(\sigma ,i\omega ;\sigma ,i\omega |\sigma ,i\omega _l;\sigma ,i\omega _l)Y(i\omega _l),`$ (68)
$`\mathrm{\Lambda }(i\omega )`$ $`=𝒢_p^{(0)}(i\omega ){\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\omega _l}{}}\varphi (i\omega _l)[\overline{𝒢}^{(0)ir}(\sigma ,i\omega ;\sigma ,i\omega _l|\sigma ,i\omega _l;\sigma ,i\omega )`$
$`+\overline{𝒢}^{(0)ir}(\sigma ,i\omega ;\sigma ,i\omega _l|\sigma ,i\omega _l;\sigma ,i\omega )].`$ (69)
In order to find a solution of Eq. (68) we insert (54), replace $`Y(i\omega _n)`$ by

$`Y(i\omega _n)`$ $`=\varphi ^{sc}(z_0)\chi (i\omega _n)Y(z_0)\text{with }z0,`$ (70)
$`\chi (i\omega _n)`$ $`={\displaystyle \frac{1}{\beta }}{\displaystyle \underset{\omega _l}{}}\overline{𝒢}^{(0)ir}(\sigma ,i\omega _n;\sigma ,i\omega _n|\sigma ,i\omega _l;\sigma ,i\omega _l)`$
$`=`$ $`{\displaystyle \frac{2\omega _c[\omega _cϵ\mathrm{tanh}(\beta ϵ/2)\mathrm{cosh}(\beta \omega _c/2)]+\mathrm{cosh}(\beta \omega _c)(\omega _c^2+ϵ^2+\omega _n^2)(\mathrm{cosh}(\beta \omega _c)1)^1}{[\omega _n^2+(\omega _c+ϵ)^2][\omega _n^2+(\omega _c+ϵ)^2]}}.`$ (71)
and make use of Poisson’s sum rule,

$$\frac{1}{\beta }\underset{\omega _n}{}f(i\omega _n)=\frac{1}{2\pi i}_C𝑑z\frac{f(z)}{e^{\beta z}+1},$$
(72)
where $`C`$ denotes the usual counterclockwise contour of the imaginary axis. We obtain then from Eq. (68) with help of the analytically continued function $`\chi (z)`$ for $`Z=Z_0=0`$ an equation for the critical temperature $`T_c`$,

$$\chi (0|ϵ)\varphi ^{sc}(0|ϵ)=1,$$
(73)
$$\chi (0|ϵ)=\left(2\omega _c[\omega _cϵ\mathrm{tanh}(\beta _cϵ/2)\mathrm{coth}(\beta _c\omega _c/2)]+\frac{(ϵ^2\omega _c^2)\mathrm{cosh}(\beta _c\omega _c)}{\mathrm{cosh}(\beta _c\omega _c)1}\right)\left(\omega _c^2ϵ^2\right)^2.$$
(74)
This quantity is even in $`ϵ`$ and therefore only the absolute value of $`ϵ=\overline{ϵ}_0\overline{\mu }`$ determines $`k_BT_c=\beta _c^1`$. From (58) and (65) one can make a rough guess for the quantity $`\varphi ^{sc}(0)`$:

$$\varphi ^{sc}(0)\left(\frac{W}{4}\right)^2\frac{1}{\gamma ^2},\gamma =\frac{1}{2}\left(1+\sqrt{1\lambda ^2(0+i\delta )}\right),$$
(75)
where $`\gamma `$ has to satisfy: $`\gamma (ϵ)=\gamma (ϵ)`$. This quantity can be obtained self-consistently from Eq. (64) for the chemical potential. For simplicity we replace here $`[1+\sqrt{1\lambda ^2(0+i\delta )}]`$ by $`2\gamma `$. Then (64) can be written as

$$\frac{2}{\gamma }\frac{1}{\beta }\underset{\omega _n}{}\mathrm{\Lambda }(i\omega _n)e^{i\omega _n0^+}=\frac{N_p}{N}.$$
(76)
Using (69) together with (52) and (53) we can express $`\mathrm{\Lambda }(i\omega _n)`$ as

$`\mathrm{\Lambda }(i\omega _n)=`$ $`{\displaystyle \frac{(i\omega _nϵ)A_1(ϵ)+\omega _cB_1(ϵ)}{(i\omega _nϵ)^2\omega _c^2}}`$
$`+{\displaystyle \frac{\omega _c^2[(i\omega _nϵ)A_2(ϵ)+\omega _cB_2(ϵ)]}{[(i\omega _nϵ)^2\omega _c^2]^2}}`$ (77)
with unknown coefficients $`A_i`$ and $`B_i`$. They can be found from Eq. (69) or more easily from the asymptotic behavior of (77) for $`|\omega _n|\mathrm{}`$,

$`\underset{|\omega _n|\mathrm{}}{\mathrm{\Lambda }(i\omega _n)}=`$ $`{\displaystyle \frac{A_1}{i\omega _n}}+{\displaystyle \frac{A_1ϵ+\omega _cB_1}{(i\omega _n)^2}}+{\displaystyle \frac{A_1(\omega _c^2+ϵ^2)+2ϵ\omega _cB_1+\omega _c^2A_2}{(i\omega _n)^3}}`$
$`+{\displaystyle \frac{A_1(ϵ^3+3ϵ\omega _c^2)+B_1(\omega _c^3+3ϵ^2\omega _c)+\omega _c^2(3ϵA_2+\omega _cB_2)}{(i\omega _n)^4}}+\mathrm{}`$ (78)
If we compare this with the asymptotic behavior of the full one-polaron Green’s function (Appendix A) by invoking the methods of moments together with the asymptotic behavior of $`g(i\omega _n)`$ in (66), we obtain

$`A_1(ϵ)`$ $`=1,`$ (79a)
$`B_1(ϵ)`$ $`={\displaystyle \frac{1}{\omega _c}}[M_1+ϵ],`$ (79b)
$`A_2(ϵ)`$ $`={\displaystyle \frac{1}{\omega _c^2}}\left[M_2+2ϵM_1+ϵ^2\omega _c^2(W/4)^2\right],`$ (79c)
$`B_2(ϵ)`$ $`={\displaystyle \frac{1}{\omega _c^3}}\left[M_33ϵM_2+M_1\left(\omega _c^23ϵ^2+3(W/4)^2\right)+ϵ\omega _c^2ϵ^3+3ϵ(W/4)^2\right],`$ (79d)
where $`M_i`$ is the $`i`$-th moment of the one-polaron Green’s function. The results in (A.5) for the moments in lowest order allow to evaluate $`A_i`$ and $`B_i`$, see (A.7). $`A_1=1`$ describes the asymptotic freedom of the polarons. $`B_1=\mathrm{tanh}(\beta ϵ/2)`$ is identical with its value in the zero-order polaron Green’s function of (38). The two new quantities $`A_2`$ and $`B_2`$ are small quantities being proportional to $`\omega _0/\omega _c=1/\alpha `$

Inserting now (77) into the left-hand part of Eq. (76) and performing the summation leads to

$`{\displaystyle \frac{1}{\beta }}`$ $`{\displaystyle \underset{\omega _n}{}}\mathrm{\Lambda }(i\omega _n)e^{i\omega _n0^+}=\overline{n}(ϵ)`$
$`+{\displaystyle \frac{\mathrm{tanh}(\beta ϵ/2)[\mathrm{tanh}(\beta \omega _c/2)1][1\mathrm{tanh}^2(\beta ϵ/2)]}{2[1\mathrm{tanh}^2(\beta \omega _c/2)\mathrm{tanh}^2(\beta ϵ/2)]}}`$
$`+{\displaystyle \frac{B_2(ϵ)}{4}}\mathrm{tanh}(\beta \omega _c/2){\displaystyle \frac{1\mathrm{tanh}^2(\beta ϵ/2)}{1\mathrm{tanh}^2(\beta ϵ/2)\mathrm{tanh}^2(\beta \omega _c/2)}}`$
$`{\displaystyle \frac{\beta \omega _c}{16}}\left[{\displaystyle \frac{A_2(ϵ)+B_2(ϵ)}{\mathrm{cosh}^2[\beta (\omega _c+ϵ)/2]}}{\displaystyle \frac{A_2(ϵ)B_2(ϵ)}{\mathrm{cosh}^2[\beta (\omega _cϵ)/2]}}\right],`$ (80)
which, according to Eq. (76), is equal to $`(\gamma /2)(N_p/N)`$. Because of the large collective frequency, $`\beta \omega _c1`$, we may omit exponentially small quantities like $`\mathrm{exp}(\beta \omega _c)`$. Since we are interested in results for electron numbers which are close to half filling ($`ϵ=0`$), also $`|ϵ|\omega _c`$ holds. Furthermore we will neglect contributions of the order $`1/\alpha `$. Then the equation for the chemical potential (74) is simply

$$\overline{n}(ϵ)=\gamma n_p/2,n_p=N_p/N.$$
(81)
If furthermore we use $`\gamma =1`$ (free polarons) we obtain with (75) for $`\varphi ^{sc}(0)`$ the value

$$\varphi ^{sc}(0=(W/4)^2$$
(82)
which allows to write the equation for the critical temperature $`T_c`$ in the following form

$$ϵ^2+\omega _c^22ϵ\omega _c\mathrm{tanh}(\beta ϵ/2)=(\omega _c^2ϵ^2)^2(4/W)^2.$$
(83)
In this approximation $`T_c`$ depends only on the local parameters; but we expect that close to half filling this should give an impression which of the local quantities is most important for superconductivity in the strong-coupling limit of the Hubbard-Holstein model. For strictly half filling Eq. (83) can only be fulfilled when $`\omega _c=W/4`$. This may perhaps be an unphysical large value for the renormalized quantity. It also shows that the specific limit $`\overline{U}=0`$ is probably the critical value for the occurence of superconductivity in the frame of the Hubbard-Holstein model. It is clear that for $`\overline{U}<0`$ superconductivity is possible; but in this case it would have to compete in energy with the energies of charge-ordered states.

For the special case $`\omega _c=W/4`$ we obtain

$$\left(ϵ/\omega _c\right)^2\left[3\left(ϵ/\omega _c\right)^2\right]=\left(2|ϵ|/\omega _c\right)\mathrm{tanh}(\beta _c|ϵ|/2).$$
(84)
Since $`(ϵ/\omega _c<3`$ holds (not discussed in detail) we may seek for solutions for the case $`|ϵ|\omega _c`$ leading to

$`k_BT_c`$ $`={\displaystyle \frac{\omega _c}{3}}\left[1{\displaystyle \frac{5}{12}}\left({\displaystyle \frac{ϵ}{\omega _c}}\right)^2+\mathrm{}\right]`$
$`={\displaystyle \frac{W}{12}}\left[1{\displaystyle \frac{20}{3}}\left({\displaystyle \frac{ϵ}{W}}\right)^2+\mathrm{}\right].`$ (85)
In spite of the many approximations (which however are all reasonable) used the result for $`T_c`$ is in so far remarkable as it shows that the critical temperature depends on the band width (corresponding to the largest cut-off energy of the model) and not on the effective mass of the ions. for small deviations from half filling $`T_c`$ decreases and does not depend on the sign of $`ϵ`$.

For different values of $`\omega _c`$,

$$\omega _c=W/4y,$$
(86)
with $`y0`$ there are only solutions for not half filling. In this case Eq. (84) may in this case be written as

$`\beta _c|ϵ|`$ $`=\mathrm{ln}{\displaystyle \frac{1+\kappa }{1\kappa }},`$ (87)
$`\kappa `$ $`={\displaystyle \frac{ϵ^2+\omega _c^2(4/W)^2(\omega _c^2ϵ^2)^2}{2|ϵ|\omega _c}},0<\kappa <1.`$ (88)
The condition $`\kappa <1`$ is equivalent to

$$|ϵ|+\omega _c<W/4$$
(89)
while $`\kappa >0`$ is related to the parameter $`y`$,

$`\omega _c<W/4,y>0`$ $`:(W/4\omega _c)^2<ϵ^2<ϵ_{\mathrm{max}}^2,`$ (90)
$`\omega _c>W/4,y<0`$ $`:ϵ_{min}^2<ϵ^2<ϵ_{max}^2,`$ (91)
$`ϵ_{max,min}^2=\omega _c^2`$ $`+(1/2)(W/4)^2\pm (W/8)\sqrt{(W/4)^2+8\omega _c^2}.`$ (92)
For small $`y`$ we can simplify (87) and (88),

$`\kappa `$ $`\mathrm{\hspace{0.33em}\hspace{0.17em}2}/(W|ϵ|)\{ϵ^2[3ϵ^2(4/W)^2]+y[W/24ϵ^2/Wϵ^4(4(W)^3]`$
$`+y^2[3+ϵ^2(4/W)^2ϵ^4(4/W)^4],\}`$ (93)
with the following restrictions for $`ϵ`$:

$`y>0`$ $`:y^2<ϵ^2<3(W/4)^2(5/6)Wy+(29/27)y^2,`$ (94)
$`y<0`$ $`:(W/6)|y|+(25/27)y^2<ϵ^2<3(W/4)^2+(5/6)W|y|+(29/27)y^2.`$ (95)
Large values for $`T_c`$ can be achieved for $`\kappa ⩽̸1`$ and in the vicinity of half filling ($`ϵ0`$):

$$k_BT_c\frac{W\delta }{12(\delta 1)},\delta =\frac{6ϵ^2}{W|y|}<1,y=\frac{W}{4}\omega _c,$$
(96)
but only for $`y<0`$ and hence $`\omega _c>W/4`$.

## VII summary

We have discussed the occurence of superconductivity for the Hubbard-Holstein model in the strong-coupling limit ($`\overline{g}1`$). Strong coupling leads to a renormalization of the one-polaron Green’s function already in the local approximation. For $`\overline{g}1`$ we find a collective mode for the phonon clouds estimated by evaluating integrals in the Laplace approximation. Due to the absorption and emission of this mode by the polarons, the on-site energies of polarons are renormalized. Similarly the irreducible two-particle Green’s are renormalized. Allowing for the exchange of polarons including their phonon clouds leads to a new irreducible Green’s function which has been used to study spin-singlet pairing of polarons. Analytical results for the superconducting phase have been obtained for the limiting case, for which the local repulsion of polarons is exactly canceled by their attractive interaction. The resulting equation for the critical temperature has been obtained by assuming a large collective-mode frequency and a nearly half-filled band case. The parameters which determine $`T_c`$ are $`\omega _c(W/4)`$, $`ϵ`$ ($`ϵ=0`$ corresponds to half filling) and the band width $`W`$ with $`T_c`$. In the strong-coupling limit we obtain a critical temperature which is of the order of $`W/12`$.

###### Acknowledgements.

This work was supported by the Heisenberg-Landau Program and by the Program of the High Council for Science and Technical Development of Moldova. It is a pleasure to acknowledge discussions with Prof. Plakida and Prof. Crisan. V. A. M. would like to thank the Universities of Duisburg and Salerno for financial support. P. E. would like to thank the Bogoliubov Laboratory of Theoretical Physics, JINR, for the hospitality he received during his stay in Dubna.

*

## Appendix A Method of moments

By using the Heisenberg representation of the one-polaron Green’ function

$`𝒢_\sigma (𝐱𝐱^{}|\tau \tau ^{})`$ $`=T\widehat{c}_{𝐱\sigma }(\tau )\widehat{\overline{c}}_{𝐱^{}\sigma ^{}}(\tau ^{})_H,`$ (97)
with

$`\widehat{c}_{𝐱\sigma }(\tau )`$ $`=e^\tau c_{𝐱\sigma }e^\tau ,`$ (98a)
$`\widehat{\overline{c}}_{𝐱\sigma }(\tau )`$ $`=e^\tau c_{𝐱\sigma }^{}e^\tau ,`$ (98b)
we can write the asimptotic expansion of the Fourier representation (25) for $`|\omega _n|\mathrm{}`$ as

$`𝒢_\sigma (𝐱=0|i\omega _n)`$ $`=g_\sigma (i\omega _n)={\displaystyle \frac{1}{i\omega _n}}{\displaystyle \frac{M_1}{(i\omega _n)^2}}+{\displaystyle \frac{M_2}{(i\omega _n)^3}}{\displaystyle \frac{M_3}{(i\omega _n)^4}}+\mathrm{},`$ (99)
$`M_n`$ $`=\{c_{𝐱\sigma }^{},\underset{n}{\underset{}{[[\mathrm{}[,c_{𝐱\sigma }]}}\mathrm{}\mathrm{}]\}_H.`$ (100)
Here the statistical average $`\mathrm{}_H`$ is defined with respect to the full density matrix of the grand canonical ensemble. In the simplest approximation we obtain for the first three moments of the Green’s functions the relations,

$`M_1`$ $`=(ϵ+\omega _c\mathrm{tanh}(\beta ϵ/2),`$ (101a)
$`M_2`$ $`=ϵ^2+\omega _c^2+(W/4)^2+2ϵ\omega _c\mathrm{tanh}(\beta ϵ/2)+\omega _0\omega _c\mathrm{coth}(\beta \omega _c/2),`$ (101b)
$`M_3`$ $`=\{ϵ^3+3ϵ[\omega _c^2+\omega _c\omega _0\mathrm{coth}(\beta \omega _c/2)+(W/4)^2]`$
$`+\omega _c\mathrm{tanh}(\beta ϵ/2)[\mathrm{\hspace{0.17em}3}ϵ^2+3(W/4)^2+\omega _c^2+\omega _0^2+3\omega _0\omega _c\mathrm{coth}(\beta \omega _0/2)]\}.`$ (101c)
The expressions for the moments can be used to determine the unknown coefficients, $`A_n(ϵ)A_n(ϵ)`$, in the relation for $`\mathrm{\Lambda }_\sigma (i\omega )`$, see Eq. (78), by considering also the asymptotic behavior of $`g_\sigma (i\omega )`$ in (66) for small values of $`\lambda _\sigma (i\omega )`$,

$$g_\sigma (i\omega )=(2/W)\lambda _\sigma (i\omega )[\mathrm{\hspace{0.17em}1}+(\lambda ^2/4)+2(\lambda ^2/4)^2+\mathrm{}].$$
(102)
We then insert the asymptotic form for $`\lambda _\sigma (i\omega )`$ from (78). Comparing the corresponding equations fixes the coefficients $`A_n(ϵ)`$ and $`A_n(ϵ)`$,

$`A_1(ϵ)`$ $`=1,`$ (103a)
$`B_1(ϵ)`$ $`={\displaystyle \frac{1}{\omega _c}}[M_1+ϵ],\mathrm{tanh}(\beta ϵ/2),`$ (103b)
$`A_2(ϵ)`$ $`={\displaystyle \frac{1}{\omega _c^2}}\left[M_2+2ϵM_1+ϵ^2\omega _c^2(W/4)^2\right]{\displaystyle \frac{\omega _0}{\omega _c}}\mathrm{coth}(\beta \omega /2),`$ (103c)
$`B_2(ϵ)`$ $`={\displaystyle \frac{1}{\omega _c^3}}\left[M_33ϵM_2+M_1\left(\omega _c^23ϵ^2+3(W/4)^2\right)+ϵ\omega _c^2ϵ^3+3ϵ(W/4)^2\right]`$
$`{\displaystyle \frac{\omega _0}{\omega _c}}\mathrm{tanh}(\beta ϵ/2)\left[{\displaystyle \frac{\omega _0}{\omega _c}}+3\mathrm{coth}(\beta \omega _0/2)\right].`$ (103d)