text:
# Superconductivity in an Exactly Solvable Model of the Pseudogap State: Absence of Self-Averaging.

## I Introduction

Among the anomalies of electronic properties of high – temperature superconducting cuprates (HTSC) especially interesting are the properties of the pseudogap state, observed in a wide region of the phase diagram . From our point of view, the preferential scenario for the formation of the pseudogap state is based on the existence (mainly in the underdoped region of the phase diagram) of strong scattering of current carriers on well developed fluctuations of “dielectric” short range order (of antiferromagnetic (AFM) or charge density wave (CDW) type) . This scattering leads to essentially non – Fermi liquid like renormalization of electronic spectrum in certain parts of momentum space close to the Fermi surface and near the so called “hot” spots or “hot” (flat) patches . The preferential nature of “dielectric” and not of “superconducting” pseudogap formation  follows from the number experiments, the appropriate discussion can be found in Ref. .

The major part of existing theoretical papers is dealing with pseudogap formation and its influence on the properties of the system in normal (non superconducting) state, while only few are consider the anomalies of superconductivity in the pseudogap state . In particular, in Ref.  we have analyzed superconductivity in a simple exactly solvable model of the pseudogap state, based on the model Fermi surface (in two dimensions) with “hot” patches . For the description of the pseudogap state we have used an exactly solvable model first developed for the one dimensional case in Ref.  and for the limit of very large correlation lengths of “dielectric” short range order. It was shown that the superconducting energy gap averaged over the fluctuations of short range order is, in general, non zero also in the temperature region above the “mean-field” temperature of superconducting transition $`T_c`$, corresponding (according to Ref. ) to the appearance of homogeneous superconducting state in a sample as a whole. This lead us to the conclusion , that for temperatures $`T>T_c`$ there appear superconducting “drops”, which exist up to the temperature $`T_{c0}`$ of superconducting transition in the absence of the pseudogap of “dielectric” nature. This effect was attributed in Ref.  to the absence of self averaging of superconducting order parameter (gap) in situation when the correlation length of short range order fluctuations is larger than the coherence length of superconductivity (the size of Cooper pairs).

Under the assumption of self averaging superconducting energy gap the effects of finiteness of the correlation length of short range order fluctuations were analyzed in Ref. , where we have considered the pseudogap influence on $`T_c`$, calculated the temperature dependence of the gap for $`T<T_c`$, and derived the Ginzburg – Landau expansion for $`TT_c`$. In this case we have used our “nearly exact” solution for the pseudogap state induced by Gaussian fluctuations of short range order, proposed first in Refs.  for one dimensional case and generalized for two dimensions in Refs. . Within this approach it seems rather difficult to find the solution without the assumption of the self averaging superconducting energy gap. Note that the problem of the presence or absence of self averaging of superconducting gap was not studied well enough in most of the previous work on disordered superconductors. In most studies the self averaging property was just assumed on “physical grounds”, with the reference on significantly different length scales on which the superconducting gap changes (coherence length $`\xi _0`$) and characteristic scales for the electronic subsystem (interatomic distance or the inverse Fermi momentum), e.g. in the impurity problem , or the correlation length $`\xi _{corr}`$ of short range order in our problem ), where it may be expected that the complete self averaging of the superconducting energy gap appears for $`\xi _{corr}\xi _0`$ . We are unaware of any work, where the problem of self averageness was analyzed within some exactly solvable model of disorder.

One of the main aims of the present work is to perform precisely this type of analysis within very simple (though probably not realistic enough) one dimensional model of the pseudogap state, induced by fluctuations of “dielectric” short range order with finite correlation length, proposed in a recent paper of Bartosch and Kopietz . An exact solution found in this work is very similar to that of our previous studies  and allows to perform a complete analysis of the self averaging properties of superconducting energy gap for two dimensional model of “hot” patches, proposed in Refs. . Besides we shall study the full temperature dependences of superconducting energy gap for a superconductor with the pseudogap of “dielectric” nature.

## II Simplified Model of the Pseudogap State.

Let us consider an exactly solvable model of the pseudogap state proposed in Ref.  using slightly different approach. Consider an electron in one dimension in a periodic potential of the following form:

$$V(x)=2D\mathrm{cos}(Qx+\varphi )$$
(1)
with $`Q=2p_Fk`$, where $`p_F`$ – is Fermi momentum, and $`kp_F`$ – is some small deviation from the scattering vector $`2p_F`$ ***This choice of the vector of AFM or CDW superstructure corresponds, in general, to the case of incommensurate ordering and fluctuations.. Electronic spectrum close to the Fermi level is taken in the usual linearized form:

$`\xi _1\xi _p=v_F(|p|p_F)\xi _{p2p_F}=\xi _p\text{(“ ”)}`$ (2)
$`\xi _2\xi _{pQ}=\xi _pv_Fk\xi _p\eta `$ (3)
where we introduced the notation $`\eta =v_Fk`$ ($`v_F`$ – Fermi velocity), which will be widely used in the following. The potential (1) can be rewritten as:

$$V(x)=De^{i2p_Fxikx}+D^{}e^{i2p_Fx+ikx}$$
(4)
where we have introduced the complex amplitude $`DDe^{i\varphi }`$.

The solution of this problem is elementary. In the “two-wave” approximation of the usual band theory the one-electron (normal) Green’s function, corresponding to the diagonal transition $`pp`$, takes the following form (in Matsubara’s representation):

$`g_{11}(i\epsilon _npp)={\displaystyle \frac{1}{i\epsilon _n\xi _1}}+{\displaystyle \frac{1}{i\epsilon _n\xi _1}}D^{}{\displaystyle \frac{1}{i\epsilon _n\xi _2}}D{\displaystyle \frac{1}{i\epsilon _n\xi _1}}+\mathrm{}=`$ (5)
$`={\displaystyle \frac{i\epsilon _n\xi _2}{(i\epsilon _n\xi _1)(i\epsilon _n\xi _2)|D|^2}}={\displaystyle \frac{i\epsilon +\xi +\eta }{(i\epsilon \xi )(i\epsilon +\xi +\eta )|D|^2}}`$ (6)
where in the last equality we have introduced the following short notations: $`\xi _p=\xi `$, $`\epsilon _n=\epsilon `$. We can also introduce the non diagonal (anomalous) Green’s function, corresponding to the Umklapp process $`ppQ`$:

$`g_{12}(i\epsilon _nppQ)={\displaystyle \frac{1}{i\epsilon _n\xi _1}}D^{}{\displaystyle \frac{1}{i\epsilon _n\xi _2}}+\mathrm{}=`$ (7)
$`={\displaystyle \frac{D^{}}{(i\epsilon _n\xi _1)(i\epsilon _n\xi _2)|D|^2}}={\displaystyle \frac{D^{}}{(i\epsilon \xi )(i\epsilon +\xi +\eta )|D|^2}}`$ (8)
Let now (1) be random field. Following Ref.  we shall consider the very special model of disorder, when we shall assume the randomness of the scattering vector deviation $`k`$, with Lorentzian distribution In fact this corresponds to a specific model of phase fluctuations of the potential (1).:

$$𝒫_k(k)=\frac{1}{\pi }\frac{\kappa }{k^2+\kappa ^2}$$
(9)
where $`\kappa \xi _{corr}^1`$ and $`\xi _{corr}`$ – is the correlation length of short range order. Phase $`\varphi `$ in (1) is also considered as random with uniform distribution between $`0`$ and $`2\pi `$:

$$𝒫_\varphi (\varphi )=\{\begin{array}{c}\frac{1}{2\pi }\text{for}0\varphi 2\pi \hfill \\ 0\text{for other values}\hfill \end{array}$$
(10)
Correlation function of random fields $`V(x)`$ at different points can be calculated directly and is given by:

$$<V(x)V(x^{})>=2D^2\mathrm{cos}[2p_F(xx^{})]\mathrm{exp}[\kappa |xx^{}|]$$
(11)
where the angular brackets denote averaging over (9) and (10). The random field with precisely this form of pair correlator was treated first in Ref. , as well as in Refs. , though in these works the Gaussian statistics of the random field was also assumed For the Gaussian random field all higher correlators of the random field $`V(x)`$ are factorized “a’la Wick” into products of pair correlators (11).. In the model under consideration the random field $`V(x)`$ is, in general, non Gaussian . The Fourier transformation of (11) takes the form of characteristic Lorentzian, which defines the effective interaction of an electron with fluctuations of short range order :

$$V_{eff}(q)=2D^2\left\{\frac{\kappa }{(q2p_F)^2+\kappa ^2}+\frac{\kappa }{(q+2p_F)^2+\kappa ^2}\right\}$$
(12)
This type of effective interaction was introduced in all papers on the pseudogap of “dielectric” nature cited above.

Green’s functions averaged over the ensemble of random fields (1) with distributions (9) and (10) are calculated by elementary integration. The average value of the anomalous Green’s function (8) is simply zero (after the averaging over (10)), which corresponds to the absence of “dielectric” long range order. The average of the Green’s function (6) is obtained by term by term integration of the perturbation series (6) with (9), so that:

$`G(i\epsilon _np)={\displaystyle \frac{1}{i\epsilon _n\xi _p}}+{\displaystyle \frac{1}{i\epsilon _n\xi _p}}D^{}{\displaystyle \frac{1}{i\epsilon _n+\xi _p+iv_F\kappa }}D{\displaystyle \frac{1}{i\epsilon _n\xi _p}}+`$ (13)
$`+{\displaystyle \frac{1}{i\epsilon _n\xi _p}}D^{}{\displaystyle \frac{1}{i\epsilon _n+\xi _p+iv_F\kappa }}D{\displaystyle \frac{1}{i\epsilon _n\xi _p}}D^{}{\displaystyle \frac{1}{i\epsilon _n+\xi _p+iv_F\kappa }}D{\displaystyle \frac{1}{i\epsilon _n\xi _p}}+\mathrm{}=`$ (14)
$`={\displaystyle \frac{i\epsilon _n+\xi _p+iv_F\kappa }{(i\epsilon _n\xi _p)(i\epsilon _n+\xi _p+iv_F\kappa )|D|^2}}`$ (15)
This is precisely the exact solution for the Green’s function given in Ref. .

In the following we shall also consider the model with fluctuating amplitude $`D`$ of the field (1), so that the appropriate Green’s function can be obtained by the averaging of (15) with some amplitude distribution $`𝒫_D(D)`$. In particular, we can take the amplitude distribution in the Rayleigh form :

$$𝒫_D(D)=\frac{2D}{W^2}\mathrm{exp}\left(\frac{D^2}{W^2}\right)$$
(16)
In this case the additional averaging of correlators (11) and (12) leads to the simple replacement $`DW`$. The average Green’s function of an electron in this case becomes:

$`G(i\epsilon _np)={\displaystyle _0^{\mathrm{}}}𝑑D𝒫_D(D){\displaystyle \frac{i\epsilon _n+\xi _p+iv_F\kappa }{(i\epsilon _n\xi _p)(i\epsilon _n+\xi _p+iv_F\kappa )|D|^2}}=`$ (17)
$`={\displaystyle _0^{\mathrm{}}}𝑑\zeta e^\zeta {\displaystyle \frac{i\epsilon _n+\xi _p+iv_F\kappa }{(i\epsilon _n\xi _p)(i\epsilon _n+\xi _p+iv_F\kappa )\zeta W^2}}`$ (18)
where $`W`$ determines now the effective width of the pseudogap. In the limit of the large correlation length of fluctuations of (1), i.e. for $`\xi _{corr}\mathrm{}`$ ($`\kappa 0`$), the solution (18) coincides with that found earlier in Refs.  for the case of the Gaussian random field. For finite $`\kappa `$ this solution coincides with that proposed in Ref.  during the formal analysis of approximations, used in Refs.  in the analysis of the general problem of an electron in the Gaussian random field with pair correlator given by (11). In Refs.  it was shown that the density of states, corresponding to the Green’s function (18), possesses a characteristic smooth pseudogap in the vicinity of the Fermi level, and the values of the density of states are quantitatively very close  (practically for all energies in the incommensurate case) to the values obtained in Ref. , as well as to the results of exact numerical calculation for the Gaussian random field, performed in Refs.  §§§Using the approach of Ref.  in the present model it is easy to find also the two particle Green’s function and accordingly the frequency dependence of conductivity . Unfortunately, the specific form of “disorder” (random field) leads to unphysical behavior at zero frequency, corresponding to an “ideal” conductor..

If the random field (1) is created by fluctuations of some kind of “dielectric” order parameter (e.g. CDW or AFM), distribution (16) may describe Gaussian fluctuations, existing at high enough temperatures . For lower temperatures, even before the appropriate long range order appears, the amplitude fluctuations of the order parameter are “frozen out” (cf. ) and we may assume the amplitude to be non random and put $`D=W`$, while phase fluctuations remain important up to much lower temperatures. Accordingly, the solution of the form of (15), leading to sharp enough pseudogap for large correlation lengths $`\xi _{corr}`$ , can be used to describe the low temperature region of fluctuations of short range order. As we do not analyze the microscopics of “dielectric” fluctuations, all the parameters, characterizing these fluctuations, such as correlation length $`\xi _{corr}=\kappa ^1`$ and amplitudes $`D`$ and $`W`$ (the energy width of the pseudogap) are considered as phenomenological. The “low temperature” and “high temperature” regime of fluctuations of short range order can, in principle, be realized at different temperatures in comparison to the temperature of superconducting transition.

Generalization to the case of two dimensional system of electrons, typical for HTSC – cuprates, can be done in the spirit of “hot patches” model for the Fermi surface, proposed in Refs. . In this case we shall assume the existence of two independent types of fluctuations of the type of (1) Note the crude analogy of this picture to the concept of phase separation in HTSC – cuprates (stripes) , if we associate the correlation length $`\xi _{corr}`$ with characteristic size (period) of stripes ., oriented along the orthogonal axes $`x`$ and $`y`$, strongly interacting only with electrons from flat parts of the two dimensional Fermi surface, orthogonal to these axes. Accordingly, we assume the factorized form of two dimensional (random) potential for electrons: $`V(x,y)=V(x)V(y)`$ . The size of these flat (“hot”) parts of the Fermi surface is determined by an additional parameter $`\alpha `$, so that $`2\alpha `$ gives the angular size of the flat part, as seen from the center of the Brillouin zone . In particular, the value of $`\alpha =\pi /4`$ corresponds to the square like Fermi surface (complete nesting), when all Fermi surface becomes “hot”. For $`\alpha <\pi /4`$ the “cold” parts of the Fermi surface appear, where we neglect the scattering on fluctuations of the “dielectric” order parameter and electrons are considered as free. In this model all the physical characteristics, determined by the integrals over the Fermi surface, consist of additive contributions from “hot” and “cold” parts. Pseudogap renormalization of electronic spectrum takes place only on “hot” parts (and close to these parts), while on “cold” part the usual Fermi liquid (gas) behavior remains .

This picture is in qualitative agreement with a numerous ARPES experiments on underdoped cuprates , which show that pseudogap anomalies appear in the vicinity of $`(0,\pi )`$ point of the Brillouin zone, and vanish as we move to the diagonal direction. The presence of flat parts on the Fermi surface of HTSC – cuprates was also reliably observed in the ARPES experiments by several independent groups .

## III Gor’kov’s equations in the Pseudogap State.

To study superconductivity for the system with pseudogap due to fluctuations of “dielectric” short range order we shall assume the simplest BCS form of pairing interaction, characterized by the attraction constant $`V`$, which is non zero in some energy interval $`2\omega _c`$ in the vicinity of the Fermi level ($`\omega _c`$ \- is the characteristic frequency of the quanta, responsible for electron attraction). We have already used the same approach in Refs. . In the present work we shall limit ourselves only to the analysis of $`s`$-wave pairing. There are no serious obstacles within our approach for the analysis of $`d`$-wave pairing, typical for HTSC – cuprates, though in this case the presence of the angular dependence (anisotropy) of superconducting gap leads  to the additional integration over the angle and to a significant growth of time necessary for numerical computations. At the same time, it was shown in Refs.  , that the pseudogap influence on superconductivity is essentially similar both for $`s`$ and $`d`$-wave cases, differing only by the scale of parameters leading to the same changes of basic characteristics of superconducting state ($`d`$-wave pairing is less stable to the dielectrization of electronic spectrum in comparison to $`s`$-wave pairing).

On “cold” parts of the Fermi surface superconductivity is described by standard BCS equations. In the following we shall concentrate on the derivation of Gor’kov’s equations for the one dimensional model, which is equivalent to the analysis of the situation at “hot” parts of the Fermi surface for two dimensional case . Green’s functions (6), (8) for one dimensional system in periodic field (1), can be written as matrix elements:

$`g_{11}={\displaystyle \frac{i\epsilon _n\xi _2}{(i\epsilon _n\xi _1)(i\epsilon _n\xi _2)|D|^2}}g_{12}={\displaystyle \frac{D^{}}{(i\epsilon _n\xi _1)(i\epsilon _n\xi _2)|D|^2}}`$ (19)
$`g_{21}={\displaystyle \frac{D}{(i\epsilon _n\xi _1)(i\epsilon _n\xi _1)|D|^2}}g_{22}={\displaystyle \frac{i\epsilon _n\xi _1}{(i\epsilon _n\xi _1)(i\epsilon _n\xi _2)|D|^2}}`$ (20)
In the presence of Cooper pairing the Gor’kov’s equations, constructed on “free” Green’s functions of the type of (20), are shown in diagrammatic form in Fig. 1. In analytic form we have:

$`G_{11}=g_{11}g_{11}\mathrm{\Delta }F_{11}^+g_{12}\mathrm{\Delta }F_{21}^+`$ (21)
$`F_{11}^+=g_{11}^{}\mathrm{\Delta }^{}G_{11}+g_{12}^{}\mathrm{\Delta }^{}G_{12}`$ (22)
$`G_{21}=g_{21}g_{21}\mathrm{\Delta }F_{11}^+g_{22}\mathrm{\Delta }F_{21}^+`$ (23)
$`F_{21}^+=g_{21}^{}\mathrm{\Delta }^{}G_{11}+g_{22}^{}\mathrm{\Delta }^{}G_{21}`$ (24)
where superconducting energy gap is determined, as usual, from:

$$\mathrm{\Delta }^{}=VT\underset{np}{}F_{11}^+(\epsilon _np)=\lambda T\underset{n}{}_{\mathrm{}}^{\mathrm{}}𝑑\xi _pF_{11}^+(\epsilon _n\xi _p)\lambda T\underset{n}{}\overline{F_{11}^+(\epsilon _n)}$$
(25)
where we have introduced the dimensionless pairing constant $`\lambda =N_0(0)V`$ ($`N_0(0)`$ – is the free electron density of states at the Fermi level.).

The solution of Eqs. (24) gives:

$`G_{11}={\displaystyle \frac{1}{Det}}[(i\epsilon +\xi _1)(\epsilon ^2+\xi _2^2+D^2+\mathrm{\Delta }^2)D^2(\xi _1+\xi _2)]=`$ (26)
$`={\displaystyle \frac{1}{Det}}\{(i\epsilon +\xi )[\epsilon ^2+(\xi +\eta )^2+D^2+\mathrm{\Delta }^2]+D^2\eta \}`$ (27)
$`F_{11}^+={\displaystyle \frac{1}{Det}}\mathrm{\Delta }^{}(\epsilon ^2+\xi _2^2+D^2+\mathrm{\Delta }^2)=`$ (28)
$`={\displaystyle \frac{1}{Det}}\mathrm{\Delta }^{}[\epsilon ^2+(\xi +\eta )^2+D^2+\mathrm{\Delta }^2]`$ (29)
where

$`Det=(\epsilon ^2+\xi _1^2+D^2+\mathrm{\Delta }^2)(\epsilon ^2+\xi _2^2+D^2+\mathrm{\Delta }^2)(\xi _1+\xi _2)^2D^2=`$ (30)
$`=(\epsilon ^2+\xi ^2+D^2+\mathrm{\Delta }^2)(\epsilon ^2+(\xi +\eta )^2+D^2+\mathrm{\Delta }^2)\eta ^2D^2`$ (31)
where $`D`$ denotes the real amplitude of fluctuation field (1). In accordance with Eq. (25) Gor’kov’s Green’s function $`F_{11}^+`$ determines the superconducting energy gap. Taking into account the random nature of “dielectric” fluctuations Eq. (25) should be averaged over the fluctuations of both “phase” $`\eta =v_Fk`$ and amplitude $`D`$, using distributions (9) and (for high temperature fluctuations) (16).

Rather long, though direct, calculation of the integral in (25) via residues gives:

$`\overline{F_{11}^+(\epsilon )}={\displaystyle \frac{\pi \mathrm{\Delta }^{}}{\sqrt{2}}}{\displaystyle \frac{1}{\sqrt{\sqrt{(\stackrel{~}{\epsilon }^2+D^2+\frac{\eta ^2}{4})^2\eta ^2D^2}+\stackrel{~}{\epsilon }^2+D^2\frac{\eta ^2}{4}}}}\left\{1+{\displaystyle \frac{\stackrel{~}{\epsilon }^2+D^2+\frac{\eta ^2}{4}}{\sqrt{(\stackrel{~}{\epsilon }^2+D^2+\frac{\eta ^2}{4})^2\eta ^2D^2}}}\right\}`$ (32)
$`\pi \mathrm{\Delta }^{}(\epsilon ,\mathrm{\Delta },\eta ,D)`$ (33)
where we have introduced:

$$\stackrel{~}{\epsilon }=\sqrt{\epsilon ^2+\mathrm{\Delta }^2}$$
(34)
Then from Eq. (25) we can immediately obtain the equation for superconducting energy gap in two dimensional “hot patches” model :

$$1=2\pi \lambda T\underset{n=0}{\overset{\left[\frac{\omega _c}{2\pi T}\right]}{}}\left\{\stackrel{~}{\alpha }(\epsilon ,\mathrm{\Delta },\eta ,D)+\frac{1\stackrel{~}{\alpha }}{\stackrel{~}{\epsilon }}\right\}$$
(35)
where we have introduced the relative fraction of “hot” parts on the Fermi surface $`\stackrel{~}{\alpha }=\frac{4}{\pi }\alpha `$. The second term in (35) gives the standard BCS-like contribution from “cold” parts, occupying $`(1\stackrel{~}{\alpha })`$ part of the Fermi surface. Summation over $`n`$ in (35) is performed up to some maximal value determined by the integer part of the ratio $`\frac{\omega _c}{2\pi T}`$.

Numerical solution of (35) allows to find the value of the gap $`\mathrm{\Delta }(\eta ,D)`$ for fixed $`\eta `$ and $`D`$ (i.e. for the given value of the random field of fluctuations (1)) for any given temperature. After that we can perform averaging over (9) and (16) and find in this way temperature dependence of the average energy gap. In particular, for the “low temperature” regime of “dielectric” fluctuations it is sufficient to average only over the “phase” $`\eta `$, so that the superconducting gap is given by:

$$<\mathrm{\Delta }>=\frac{1}{\pi }_{\mathrm{}}^{\mathrm{}}𝑑\eta \frac{v_F\kappa }{\eta ^2+v_F^2\kappa ^2}\mathrm{\Delta }(\eta ,D)$$
(36)
In “high temperature” regime we have to add also the averaging over the amplitude $`D`$ with distribution function (16):

$$<\mathrm{\Delta }>=\frac{2}{W^2}_0^{\mathrm{}}𝑑DD\mathrm{exp}\left(\frac{D^2}{W^2}\right)\frac{1}{\pi }_{\mathrm{}}^{\mathrm{}}𝑑\eta \frac{v_F\kappa }{\eta ^2+v_F^2\kappa ^2}\mathrm{\Delta }(\eta ,D)$$
(37)
As a result we shall find the temperature dependences of the average superconducting gap $`<\mathrm{\Delta }>`$ without any statistical assumptions like the self averaging nature of the order parameter. Analogously we can calculate the temperature dependence of dispersion $`<\mathrm{\Delta }^2><\mathrm{\Delta }>^2`$, allowing to estimate the randomness of $`\mathrm{\Delta }`$, i.e. the presence or absence of self averaging. Results of these calculations will be presented in the next section.

As we have already noted in the Introduction, most papers on superconductivity in disordered systems actually assume the self averaging nature of the superconducting gap $`\mathrm{\Delta }`$. In this case $`\mathrm{\Delta }`$ is treated, in fact, as non random and independent of parameters of the random field in which the electrons forming the Cooper pairs actually propagate. In our case these parameters are the amplitude $`D`$ and “phase” $`\eta `$ of (1), accordingly the self averaging over these characteristics may be studied separately.

Let $`\mathrm{\Delta }`$ be self averaging over the fluctuations of $`\eta `$. In this case we can treat $`\mathrm{\Delta }`$ in Eq. (29) is independent of $`\eta `$. Then the anomalous Gor’kov’s function averaged over the fluctuations of $`\eta `$ can be written as:

$$<F_{11}^+>=\frac{\mathrm{\Delta }^{}}{\pi }_{\mathrm{}}^{\mathrm{}}𝑑\eta \frac{v_F\kappa }{\eta ^2+v_F^2\kappa ^2}\frac{\epsilon ^2+(\xi +\eta )^2+D^2+\mathrm{\Delta }^2}{(\epsilon ^2+\xi ^2+D^2+\mathrm{\Delta }^2)[\epsilon ^2+(\xi +\eta )^2+D^2+\mathrm{\Delta }^2]\eta ^2D^2}$$
(38)
This integral can be directly calculated, so that after long calculations we get:

$$<F_{11}^+>=\mathrm{\Delta }^{}\frac{\stackrel{~}{\epsilon }^2\left(1+\frac{v_F\kappa }{\stackrel{~}{\epsilon }}\right)^2+D^2\left(1+\frac{v_F\kappa }{\stackrel{~}{\epsilon }}\right)+\xi ^2}{\left[\left(1+\frac{v_F\kappa }{\stackrel{~}{\epsilon }}\right)\stackrel{~}{\epsilon }^2+\xi ^2+D^2\right]^2+v_F^2\kappa ^2\xi ^2}$$
(39)
where $`\stackrel{~}{\epsilon }`$ was introduced in (34). Accordingly we can calculate the integral of (39), entering the gap equation:

$$\overline{<F_{11}^+>}_{\mathrm{}}^{\mathrm{}}𝑑\xi <F_{11}^+>=\frac{\pi \mathrm{\Delta }^{}\left(1+\frac{v_F\kappa }{2\stackrel{~}{\epsilon }}\right)}{\sqrt{D^2+\stackrel{~}{\epsilon }^2\left(1+\frac{v_F\kappa }{2\stackrel{~}{\epsilon }}\right)^2}}$$
(40)
So, despite rather complicated form of the anomalous Green’s function (39), in the gap equation the account of the interaction with fluctuations on “hot” (flat) parts of the Fermi surface reduces to more or less “standard” renormalization:

$`\epsilon \epsilon \left(1+{\displaystyle \frac{v_F\kappa }{2\stackrel{~}{\epsilon }}}\right)=\epsilon \left(1+{\displaystyle \frac{v_F\kappa }{2\sqrt{\epsilon ^2+\mathrm{\Delta }^2}}}\right)`$ (41)
$`\mathrm{\Delta }\mathrm{\Delta }\left(1+{\displaystyle \frac{v_F\kappa }{2\stackrel{~}{\epsilon }}}\right)=\mathrm{\Delta }\left(1+{\displaystyle \frac{v_F\kappa }{2\sqrt{\epsilon ^2+\mathrm{\Delta }^2}}}\right)`$ (42)
analogous to that appearing in the problem of impurity scattering in superconductors . Similar renormalization was already noted for the variant of the present problem in Ref. . The analogy with impurity problem here is rather natural as our parameter $`v_F\kappa =v_F\xi _{corr}^1`$ actually determines the characteristic inverse time of flight of an electron through the region of short range order of the size of $`\xi _{corr}`$. Of course, the additional pseudogap influence is also connected with the appearance in Eqs. (39), (40) of the square of dielectric gap $`D^2`$.

Finally, the gap equation determining superconductivity in the “hot patches” model with the assumption of self averaging over “phase” fluctuations takes the following form:

$$1=2\pi \lambda T\underset{n=0}{\overset{\left[\frac{\omega _c}{2\pi T}\right]}{}}\left\{\stackrel{~}{\alpha }\frac{1+\frac{v_F\kappa }{2\stackrel{~}{\epsilon }}}{\sqrt{D^2+\stackrel{~}{\epsilon }^2\left(1+\frac{v_F\kappa }{2\stackrel{~}{\epsilon }}\right)^2}}+\frac{1\stackrel{~}{\alpha }}{\stackrel{~}{\epsilon }}\right\}$$
(43)
The solution of this equations, naturally, is simpler than that of Eq. (35) with later averaging over (36). In the absence of fluctuations of the amplitude of “dielectric” field $`D`$, which is valid for the “low temperature” region of fluctuations of short range order, it is Eq. (43) that determines “mean field” (in terms of Ref. ) behavior of $`\mathrm{\Delta }(T)`$ with respect to fluctuations of the random field (1).

In “high temperature” region of fluctuations of short range order, assuming for $`D`$ the distribution function (16) and self averaging over the fluctuations of $`D`$, we obtain the following equation for the average superconducting gap:

$$1=2\pi \lambda T\underset{n=0}{\overset{\left[\frac{\omega _c}{2\pi T}\right]}{}}\left\{\frac{2\stackrel{~}{\alpha }}{W^2}_0^{\mathrm{}}𝑑DD\mathrm{exp}\left(\frac{D^2}{W^2}\right)\frac{1+\frac{v_F\kappa }{2\stackrel{~}{\epsilon }}}{\sqrt{D^2+\stackrel{~}{\epsilon }^2\left(1+\frac{v_F\kappa }{2\stackrel{~}{\epsilon }}\right)^2}}+\frac{1\stackrel{~}{\alpha }}{\stackrel{~}{\epsilon }}\right\}$$
(44)
describing the situation analogous to that studied in our previous work , where we have analyzed the influence of Gaussian fluctuations of “dielectric” short range order using the methods of Refs. . Here the fluctuations of (1) are treated exactly, but $`D`$ is assumed to be self averaging. We shall see below that the results following from the solution of Eq. (44) are very close to that obtained in Ref. . For $`\kappa 0(\xi _{corr}\mathrm{})`$ Eq. (44) reduces to the similar “mean field” equation of Ref. . The temperature of superconducting transition, determined by Eq. (43) or Eq. (44), can apparently be identified as the temperature of the appearance of infinitesimally small superconducting gap, homogeneous over the whole sample .

In the next section we shall present the results of numerical solution of Eqs. (43), (44) in comparison with the results of an exact analysis, based upon the approach using Eqs. (35), (36), (37).

## IV Main Results and Discussion.

Let us discuss the results of numerical analysis of gap equations derived in the previous sectionDuring our calculations we have assumed for the relative part of flat patches on the Fermi the value of $`\stackrel{~}{\alpha }=2/3`$.

In Figs. 2 and 3 we show the dependences of the critical temperature $`T_c`$ of superconducting transition for the “low temperature” regime of “dielectric” fluctuations (at this temperature the “mean field” gap, determined by Eq. (43) becomes zero) on the width of the pseudogap $`W`$ (which in this case coincide with the amplitude of dielectric gap $`D`$) and correlation length of short range order. These results are in qualitative agreement with similar dependences for the “high temperature” regime of “dielectric” fluctuations (when $`T_c`$ is determined by Eq. (44)), as well as with dependences, obtained by us earlier for the different model of “dielectric” fluctuations with finite correlation length in Ref. . With the growth of the pseudogap width $`W`$ the “mean field” value of $`T_c`$ is suppressed. Diminishing correlation length leads to the “filling” of the pseudogap  and $`T_c`$ suppression becomes less effective.

In Fig. 4 the curves show the temperature dependences of superconducting gap $`<\mathrm{\Delta }>`$, averaged over both amplitude $`D`$ and “phase” $`\eta `$ (“high temperature” region of fluctuations of short range order, where $`<\mathrm{\Delta }>`$ is given by (37)), for different values of $`v_F\kappa `$. Dashed curves give the appropriate “mean field” temperature dependences of superconducting gap, obtained assuming the self averaging of superconducting order parameter over both the fluctuations of amplitude and phase, as determined by Eq. (44).

Superconducting gap, averaged over the fluctuations of short range order, is non zero in the temperature region above $`T_c`$, corresponding to the zero of the “mean field” gap (i.e. the gap homogeneous over the sample). More so, it is seen that the average gap is non zero also in a narrow region of temperatures larger than the transition temperature in the absence of short range order fluctuations (pseudogap) $`T_{c0}`$. This effect is due to the existence of fluctuations of the “phase” $`\eta `$, when the Fermi level is in the energy interval, corresponding to the peaks of the density of states due to formation of dielectric gap. To understand this, note that for the given realization of “phase” $`\eta `$ and of amplitude $`D`$, the density of states has the form:

$$\frac{N(E)}{N_0(0)}=\frac{1}{\pi N_0(0)}Im\underset{𝐩}{}g_{11}^R(Epp)=\{\begin{array}{c}\frac{|E+\frac{\eta }{2}|}{\sqrt{(E+\frac{\eta }{2})^2D^2}}\text{for}|E+\frac{\eta }{2}|>D\hfill \\ 0\text{for other values}\hfill \end{array}$$
(45)
where $`g_{11}^R(Epp)`$ – is the retarded Green’s function, obtained from (6) via standard analytical continuation $`i\epsilon _nE+i0`$, and $`N_0(0)`$ – is the density of states at the Fermi level in the absence of short range order fluctuations. Thus for $`\frac{\eta }{2}D`$ the Fermi level position approximately coincides with the peak of the density of states, leading to larger values of superconducting gap $`\mathrm{\Delta }(\eta ,D)`$. The growth of the dielectric gap $`D`$ leads also to the larger width of peaks in the density of states (45), so that for $`\frac{\eta }{2}D`$ superconducting gap $`\mathrm{\Delta }(\eta ,D)`$ grows with $`D`$. Then, for any temperature above $`T_{c0}`$ and for large enough amplitudes of dielectric gap $`D>D^{}(T)`$, on the phase diagram in the variables $`\eta `$ and $`D`$ there is always a narrow region close to the line $`\frac{\eta }{2}=D`$, where superconducting gap $`\mathrm{\Delta }(\eta ,D)`$ is different from zero (cf. Fig.5). This leads to the appearance in the temperature dependence of the averaged (over random field fluctuations) superconducting gap $`<\mathrm{\Delta }>`$ of an infinite, though exponentially small “tail” in the temperature region above $`T_{c0}`$.

At the insert in Fig. 4 we show the temperature dependence of the relative mean square fluctuation of the superconducting gap $`\delta \mathrm{\Delta }/\mathrm{\Delta }=\sqrt{<\mathrm{\Delta }^2><\mathrm{\Delta }>^2}/<\mathrm{\Delta }>`$ for the “high temperature” regime of “dielectric” fluctuations. For large correlation lengths of short range order ($`\xi _0/\xi _{corr}1`$) these fluctuations of superconducting order parameter are very strong for all temperatures, signifying the obvious absence of self averaging. Surprisingly, these fluctuations of superconducting energy gap are strong enough also for small enough correlation lengths, at least in the region of $`T>T_c`$. In particular, the “tail” in the temperature dependence of $`<\mathrm{\Delta }>`$ for $`T>T_c`$ is observed even for $`v_F\kappa /T_{c0}=100`$, when $`\xi _0/\xi _{corr}301`$.

The full curves in Fig. 6 show the temperature dependence of superconducting gap $`<\mathrm{\Delta }>`$, averaged over the “phase” $`\eta `$ (cf. (36)), for the “low temperature regime” of “dielectric” fluctuations, when the amplitude fluctuations of dielectric gap are frozen out and $`D=W`$. Dashed curves show the appropriate temperature dependences of the “mean field” superconducting gap, obtained with the assumption of self averaging superconducting order parameter over fluctuations of the “phase” $`\eta `$, and defined by Eq. (43). For large enough correlation lengths of the short range order the average gap for $`T<T_c`$ is very close to its “mean field” values and its “tail” in the region of $`T>T_c`$ is relatively small. This behavior for the “low temperature” region of “dielectric” fluctuations is due to the fact, that for $`\xi _{corr}\mathrm{}`$ there is no randomness in this model at all ($`\eta =0`$, $`D=W`$). Accordingly, the mean square fluctuation of the gap, shown at the insert in Fig. 6, is rather small for large correlation lengths and $`T<T_c`$, but grows sharply for $`T>T_c`$. As correlation length becomes smaller, fluctuations of superconducting gap $`\delta \mathrm{\Delta }`$ for $`T<T_c`$ at first grow , mainly due to the growth of randomness (parameter $`v_F\kappa `$ determines the width of the distribution of fluctuations of $`\eta `$), but afterwards diminish in the region of $`\xi _0/\xi _{corr}1`$. In the “tail” region of the averaged superconducting gap ($`T>T_c`$) fluctuations of superconducting gap are quite large. Though these fluctuations diminish for smaller correlation lengths of the short range order $`\xi _{corr}`$, they still remain significant even for small enough correlation lengths, i.e. for $`\xi _0/\xi _{corr}1`$.

As well as in the “high temperature” region of “dielectric” fluctuations the “tail” in the temperature dependence of the average gap is observed here also in the region of $`T>T_{c0}`$. It is explained by the same reasons as discussed above. However, for the “low temperature” region the amplitude of the dielectric gap is non random and fixed at $`D=W`$. Thus, for $`T_{c0}<T<T_c^{}`$, where $`T_c^{}`$ is defined by $`D^{}(T_c^{})=W`$, there exists a narrow region of “phases” close to $`\eta =2W`$, where the superconducting gap $`\mathrm{\Delta }(\eta ,W)`$ is different from zero, while for $`T>T_c^{}`$ such region is absent (cf. Fig.5). The value of $`T_c^{}`$ determines the temperature up to which exists the “tail” of the average gap, i.e. the critical temperature for the average gap $`<\mathrm{\Delta }>`$. From the definition of $`T_c^{}`$ it is obvious that this temperature is independent of correlation lengths and depends only on $`W`$. As the width of the peaks in the density of states (45), as well as $`\mathrm{\Delta }(\eta ,D)`$, grows with the growth of $`D`$ if condition $`\frac{\eta }{2}D`$ is fulfilled, the value of $`T_c^{}`$ grows with $`W`$. The dependence of $`T_c^{}`$ on $`W`$ is shown at the insert in Fig. 6.

## V Conclusion

In this work we have studied superconductivity within very simple model of the pseudogap in two dimensional electronic system which allows an exact solution. The central result is the explicit demonstration of the absence\[21 e of complete self averaging of superconducting order parameter (energy gap) over the random field of “dielectric” fluctuations, leading to the formation of the pseudogap. This is rather surprising from the point of view of the standard theory of superconductivity in disordered systems . The absence of self averaging is reflected in the appearance of fluctuations of the gap, especially strong for the temperatures larger than the “mean field” temperature of superconducting transition $`T_c`$, following from the standard equations, assuming the the self averaging nature of the order parameter. We identify this temperature with the critical temperature for the appearance of homogeneous superconducting state in the whole sample, while for the real disordered system the superconducting state is inhomogeneous and for $`T>T_c`$ superconductivity exists in separate regions (“drops”), appearing due to random fluctuations of the local density of electronic states. The difference with our previous work , where this picture was discussed in the limit of very large correlation lengths of short range order $`\xi _{corr}\mathrm{}`$, the use of the model of Ref.  allowed to find the complete solution for arbitrary values of $`\xi _{corr}`$. This solution demonstrates the absence of complete self averaging of superconducting gap even for $`\xi _{corr}<\xi _0`$, which contradicts the naive expectations of the standard approach . As we noted above, we are unaware of works, where the problem of self averaging of $`\mathrm{\Delta }`$ was analyzed within any exactly solvable model of disorder. Here we presented precisely this type of analysis. It is unclear at present whether the results obtained will be conserved in more realistic models of disorder.

It will be very interesting to analyze also the behavior of the spectral density and tunneling density of states, analogous to that done in our previous work  in the limit of $`\xi _{corr}\mathrm{}`$. It will be especially important to consider the self averaging properties of the density of states, which is the standard assumption in the theory of disordered systems.

As to comparison with experiments on high – temperature superconductors, we note Refs. , where the scanning tunneling microscopy (measuring the local density of states) on films of $`Bi_2Sr_2CaCu_2O_{8+\delta }`$ clearly demonstrated the existence in this system of microscopic superconducting regions along with dominant semiconducting regions with typical pseudogap in the electronic spectrum. These observations are in qualitative agreement with the main conclusions of the present study.

This work was partly supported by the grants 99-02-16285 of the Russian Foundation of Basic Research and CRDF No. REC-005, as well as by the program of fundamental research of the Presidium of the Russian Academy of Sciences “Quantum macrophysics”, and by the State Contracts of the Ministry of Industry and Science 108-11(00) (program “Statistical Physics”) and 107-1(00) (program of HTSC research).


format: text

text:
# Ab initio Calculations of Multilayer Relaxations of Stepped Cu Surfaces

## I INTRODUCTION

Structural properties of regularly stepped metal surfaces have been the focus of a broad range of theoretical and experimental studies because of the eminent role they play in technologically important phenomena such as thin film growth, epitaxial layer formation, nanostructuring of material, and catalysis . According to crystallographic notation, these surfaces are denoted by high Miller indices and are called vicinals of their low Miller index counterparts (flat surfaces). The presence of arrays of atomic steps separated by flat terraces creates regions of differing local coordination and makes the microscopic structure of a vicinal surface distinct from that of a flat surface. According to Smoluchowski’s idea of charge smoothing , for example, electronic charge densities are expected to rearrange in the vicinity of the steps, thereby causing the ion cores to relax to new configurations. The modified electronic structure may also be expected to impact the reactivity and the nature of the force fields in the region around the steps. Knowledge of atomic relaxations in the equilibrium positions near the step and kink sites is thus a step towards understanding the novel vibrational and electronic properties of vicinal surfaces.

Fortunately with advances in atomic scale experimental techniques, there has been a surge in investigations of the structure of vicinal surfaces in recent years. The majority of the experimental data have undoubtedly come  from the low-energy electron diffraction (LEED) technique which is now capable of detecting changes even for interlayer spacings smaller than 1$`\AA `$. For some surfaces the X-ray scattering technique has provided much needed complementary structural data . An impressive number of theoretical calculations of multilayer relaxations  have also helped in bringing several issues related to the characteristics of vicinal surfaces to the forefront. Of particular interest here are experimental and theoretical studies of a set of vicinals of Cu(100) and Cu(111) which have addressed the question of the impact of local coordination on the structural and dynamical properties of the surface. In an earlier paper , a comparative study of the local structural and vibrational properties of Cu(211), Cu(511), and Cu(331) was performed using empirical potentials from the embedded atom method (EAM) . This study found that the first two surfaces displayed similar local characteristics, while the third surface was somewhat different. An explanation provided for this behavior was the similarity in the local environment of the (211) and (511) surfaces of fcc metals (a combination of (100) and (111), terrace geometry, and step face), and its consequent difference from that of the (331) surface( a (111) terrace geometry and a (111)-microfacetted step face). The issue of the impact of the local geometry was further raised in a joint theoretical and experimental study of the vibrational dynamics of Cu(211) and Cu(511) together with those of the kinked surface Cu(532) . Experimental data from Electron Energy Loss Spectroscopy (EELS) found modes above the bulk band on Cu(211) but not on Cu(511) (or on Cu(17,1,1) which has the same step geometry as Cu(511)), but theoretical calculations based on EAM potentials predicted modes (slightly) above the bulk phonon spectrum for each of these surfaces. While the similarity between the calculated structural relaxation patterns of Cu(211) and Cu(511) argues in favor of a similarity in the local vibrational dynamics of these two surfaces, the disagreement between the experimental and the theoretical results for the high frequency modes on Cu(511) (and Cu(17,1,1)) remains unreconciled. For Cu(211) agreement of the EAM based results with available structural data from LEED  and with ab initio calculations for both the structure and the dynamics  provides considerable confidence in its predicted properties. The case of Cu(511) is not as simple because of lack of calculations based on potentials more accurate than EAM, and because of conflicting conclusions from the analysis of experimental data from LEED  and X-ray scattering measurements . The most striking difference in these two sets of data is the relaxation pattern for the second layer which is inwards in LEED and outwards in the X-ray data. The oscillatory pattern found in the X-ray data is also in disagreement with the conclusion from a series of previous experimental and theoretical findings on stepped surfaces. Based on these studies , there is a definite symmetry in the relaxation patterns of stepped surfaces. All terrace atoms, save for the corner one, display inward relaxations. The EAM based calculations  further predict this oscillatory relaxation pattern to continue into the bulk with a damping in the amplitude . Thus the expected relaxation pattern for the (211), (511), (331) surfaces, each with 3-atom wide terraces, would be (- - + - - +…), although questions have been raised whether Cu(331) follows this rule . Similarly, the patterns for (711) and (911) with, respectively, 4 and 5 atoms on the terrace, would be predicted to be (- - - + - - - +…) and (- - - - + - - - - +…). The LEED data on the first three surfaces follow these predicted trend in relaxations, atleast for the top 3 layers. The very recent LEED data  for Cu(711) also displays the pattern (- - - +) for the top layers, in good agreement with EAM based predictions. However, a small discrepancy in the sign of the relaxation is found for both Cu(711) and Cu(511), for a particular layer separation (d56 for Cu(511) and d78 for Cu(711)) , beyond that expected from the error bars. Arguably the actual numbers involved in these comparisons are small, but the systematic nature of the discrepancies and the fact that it negates the prediction of a periodicity in the oscillatory relaxation pattern , raise interesting questions about the complexities of the atomic displacements in these systems. Given the above uncertainties arising from experimental observations, it is opportune to carry out more accurate calculations of these relaxation patterns using techniques which are capable of revealing the accompanying changes in the surface electronic structure. It is with this goal in mind that we have carried out ab initio electronic structure calculations of the surface geometry and interlayer spacing for a set of vicinals of Cu(100) and Cu(111). In addition to Cu(211), Cu(331) and Cu(511) which are included to address the question of the influence of the local geometry on the structure, we have extended the investigation to Cu(711) and Cu(911) to examine the influence of increasing terrace width of the relaxation pattern. Of course, for all surfaces comparison with available experimental data is of prime concern.

The rest of this report is organized as follows. In Section II, the system geometries are presented together with some computational details. Section III contains the results and their discussion. Concluding remarks are presented in Section IV.

## II SURFACE GEOMETRIES

Vicinal surfaces can easily be constructed by cutting the crystal at an angle slightly away from the lower-index crystal planes (i.e. (100), (111), (110)). For reasons discussed above, we are interested here in the vicinals of the (100) and (111) surfaces of fcc metals of which the most tightly packed steps are along the $`<110>`$ direction. In the case of the (111) surface, however, the $`<110>`$ direction is not parallel to any plane of symmetry and there are two different ways of generating monoatomic stepped surfaces. In one type of such vicinals, the step edge has a (100)-microfacet, while the other has the (111)-microfacet (these are the so called A and B types, respectively). In the standard nomenclature, the vicinals of fcc(111) surface with monoatomic steps and (100) step edges are denoted by $`(n,n,n+2)`$, while those with (111) step edges are labeled as $`(n,n,n2)`$, where $`n`$ is the number of atoms on the terrace. The B-type vicinal Cu(331) considered here is named accordingly, while the A-type vicinal Cu(211) seems to be a misnomer. Similarly, the vicinals of fcc(100) surface consisting of monoatomic step edges with (111) microfacet are labeled $`(2n1,1,1)`$. The Cu(111) vicinals considered here are created by cutting the crystal at an angle of 19.5o and 22o away from the (111) plane towards the \[$`2\overline{1}\overline{1}`$\] and \[$`\overline{2}11`$\] direction to produce the (211) and (331) surfaces, respectively, whereas the three vicinals of (100), (511), (711), and (911), are constructed by slicing the crystal at angles of 15.8o, 11.4o and 8.9o, respectively, from the (100) plane towards the  direction. To facilitate the discussion we have also labeled the atoms that play the dominant role in our calculation . For the three surfaces consisting of three chains of atoms on the terrace we label them as corner-chain (CC), terrace-chain (TC1), and step-chain (SC). The chain just underneath the step-chain is called a bulk nearest neighbor chain (BNN). The other two surfaces, Cu(711) and Cu(911), contain, respectively, one and two extra chains of terrace atoms, labeled accordingly as (TC2) and (TC3). We have taken the $`x`$ and $`y`$ axes to lie in the surface plane, the $`x`$-axis being perpendicular to the step and the $`y`$-axis along the step, and the $`z`$-axis is along the surface normal. In Fig. 1, we display a side view of the (511) surface with the appropriate labeling of the atoms and interlayer spacing.

## III SOME DETAILS OF THEORETICAL CALCULATIONS

The ab initio electronic structure calculations are performed within a pseudopotential approach to density-functional theory in the local density approximation , numerical implementation of the technique is based on a computer code developed by B. Meyer et al.  The local-density approximation is applied using the Hedin-Lundqvist form of the exchange-correlation functional . A norm conserving pseudopotential for Cu constructed according to a scheme proposed by Hamann-Schlueter-Chiang  has been used which has already been successfully employed for calculations of the structure and the phonons of low index surfaces of Cu . A mixed basis set is applied to represent the valence states consisting of five d-type local functions at each Cu site, smoothly cut off at a radius of 2.3 a.u., and of plane waves with kinetic energy of 11 Ry. The Brillouin-zone(BZ) integration was carried out using the special point sampling technique  together with a Gaussian broadening of the energy levels of 0.2 eV. For simulating surfaces we used the supercell approach with cells containing 21 to 35 atoms (1 atom per layer), depending on the surface orientation. The z-dimension of all cells was 47.7155 a.u. The distance between the top and the bottom layer of the slabs were thus 31.2761 a.u. for Cu(331), 27.8277 a.u. for Cu(211), 26.2340 a.u. for Cu(511), 25.7715 a.u. for Cu(711), and 25.4390 a.u. for Cu(911). Structure optimization was carried out until forces on all atoms were smaller than $`10^3`$ Ry/a.u., which is two orders of magnitude smaller than the forces present on the unrelaxed surfaces. With increasing terrace width, the calculations become increasingly tedious since the reduction in interlayer spacing makes it more difficult to achieve geometries converged to 1% of the interlayer spacing. We also find that results for terraces with (100) geometry are more sensitive to the number of k-points sampled, as compared to those with (111) geometry. For the latter case 30 points in the BZ are sufficient for the determination of the equilibrium structure, while for surfaces with (100) terraces at least twice as many points are needed to get converged results.

## IV RESULTS AND DISCUSSIONS

Our results for the multilayer relaxations of Cu(211), Cu(331), Cu(511), Cu(711) and Cu(911) are summarized in Table I. As in previous theoretical studies  of relaxations on stepped Cu surfaces, we find changes in interlayer separations, from bulk terminated configurations, to persist on all surfaces for a large number of layers. Of course, the concept of layers is different here from that on flat surfaces (see Fig. 1). The first n  layers, for example, are all exposed to the vacuum, where n  is the number of atoms on the terrace. Correspondingly, the interlayer separations are small and even large percentage changes in interlayer separations correspond to small numbers in distances. It is important to bear this point in mind when comparing the results for a particular surface either with those for flat surfaces or with the results from other methods.

A common feature of all surfaces examined in Table I is that all terrace atoms except for the ones in the corner chain (CC) undergo significantly large inwards relaxations. The corner atoms are always found to relax outwards. Additionally, the atoms in the terrace adjacent to CC exhibit comparatively large inward relaxation whose magnitude maybe larger than that of the step atoms (SC). For example, for Cu(511), Cu(711), and Cu(911) changes, respectively, in $`d_{23}`$ (involving TC1), $`d_{34}`$ (involving TC2), and $`d_{45}`$ (involving TC3), are larger than that of $`d_{12}`$. Thus, in keeping with Schmolkowski’s  ideas of charge smooting, the maximum relative change in interlayer separation is focussed around the corner atoms. This is particularly true for the vicinals of Cu(100). The situation with the more closepacked surface Cu(331) is somewhat different, As seen in Table I the outward relaxation of the corner atom and the inward relaxation of the preceeding atom on the terrace on Cu(331) are less than half of that for similar atoms on the other surfaces considered here. Incidently, this conclusion is in good agreement with results from previous studies which were based on semi-empirical potentials .

There is, however, a disconcerting difference in the results obtained here from ab initio electronic structure calculations and those from semi-empirical potentials. An intriguing result for multilayer relaxations of the vicinals of Cu(100) obtained with EAM potentials in Ref.  was that the pattern of inward and outward relaxations continued well into the bulk with an expontially decreasing amplitude. Thus for Cu(711), the relaxation pattern predicted by EAM was (- - - +,- - - +,- - - +,…) with an eventual dampimg of the relaxations. The pattern for Cu(711) from Table I is instead (- - - +,- - + +,- - + +). That is our present calculations do not predict a periodically oscillatory relaxation pattern with a decaying amplitude as we move into the bulk. As we shall see, this particular feature is more in agreement with experimental data and help remove the slight discrepancy between experiment and theory presented by the EAM result pointed by Walter et al. . Again, it should be recalled that the numbers involved are very small and within the limits of accuracy of ab initio calculations. In particular, the small numbers for the relaxations of the inner layers of Cu(911) have to be taken with caution as our convergence criteria for this surface were not as good as that for the others because of the demands on computational resources imposed by a system as large as this one.

Unlike flat surfaces, vicinal surfaces relax in both $`x`$ and $`z`$ directions, since the existence of steps at the surface leads to broken symmetry in both of these directions. While relaxations along the $`z`$-direction yield characteristic interlayer separations that we have discussed above, those along the $`x`$-direction provide new registries of atoms, as compared to those in the bulk. Our calculated percentage intralayer registries for five surfaces are summarized in Table II. As in the observations from EAM calculations, the changes in the registries of the atoms are small. It is thus not useful to make a one-to-one comparison with results from semiempirical calculations. However, the changes in registries of the atoms are not inconsequential since they affect the changes in the bond lengths between the atoms in these regions of low coordinations. In Table III, we tabulate our results for the total changes in the distances between the step atoms and their nearest neighbors. For comparison we have included in parenthesis the results obtained earlier for the same quantities with EAM potentials . The largest changes in the bond lengths (from unrelaxed configurations) are for those between the step atoms and their bulk nearest neighbour (BNN) which lies right below them. The bonds between CC and BNN show small enlargement, while all other bonds in Table III are found to undergo shortening. In Figs. 2 and 3, we have drawn the actual displacements of the atoms on the five surfaces obtained from our calculations. While the size of the arrows are exagerated, it is their relative length and direction that is of consequence. As already noted by Durukanoglu et al. , all atoms in low coordinated sites move to enhance their local coordination. The complex displacement pattern that emerges is thus the net outcome of the competition between the different directions in which the various atoms would like to relax to enhance thein own coordination. For readers who are interested in the exact positions of the atoms on the relaxed surfaces, we have summarized them in Table IV. The unusual behavior of Cu(331) terrace atom is more apparent from this Table than the earlier one on changes in the bond lengths. The TC1 atom of Cu(331) undergoes the least displacement among its counterparts. Its displacement is also smaller than that of TC2 (1.4, 0.0, -2.6) on Cu(711), and of both TC2 (0.8, 0.0, -2.0) and TC3 (0.1, 0.0, -2.3) on Cu(911). The coordinates of the displacements above in paranthesis are in the same units as those in Table IV.

We now turn to comparisons of the results obtained here for individual surfaces with those available from experimental measurements. In Table V we show that for Cu(211) the salient features in the trends in the relaxation patterns predicted by our calculations are observed in the experimental data. Apart from the large inward relaxation of the step atoms, the major change occurs at the corner atom and its adjacent terrace atom. Our results are in good agreement with previous DFT/LDA calculations , based on the pseudopotential approximation and with results from EAM based method. Theoretical calculations using the full potential linearized augmented plane wave (FLAPW) method , however, predict much larger relaxation(-28.4%) of the step atom than any of the previous theoretical or experimental studies. This brings us to the discussion of ionic relaxation on Cu(331) in Table VI for which also we do not get the large relaxation reported in Ref\[\]. The present results for Cu(331) are, however, in good agreement with the LEED data . The discrepancy with the results from LEED for $`d_{23}`$ should not be taken too seriously, given an error bar of 4% in the analysis of the LEED data. With respect to EAM based results , we find a noticeable difference for d34, for which the present results agree better with the LEED data and also preserve the predicted relaxation pattern (- - +) for the terrace atoms. This trend is inkeeping with what was reported in calculations on Al(331) . In trying to reconcile our results with those of Geng et al. , we note that the latter predict an outward displacement of the TC1 atomic chain for both Cu(211) and Cu(331), while we find this not to be the case. As already mentioned, while the changes in the bond lengths of the terrace atoms of Cu(331) are no different from those of the other surfaces, the displacement of TC1 is strikingly smaller than that of the TC’s on other surfaces.

The case of multilayer relaxations for Cu(511) is interesting because of the differences in the published data from LEED  and X-ray measurements . These are displayed in Table VII. Except for the displacement of the step atoms, for which all results point to a large inward relaxation, the results from X-ray scattering measurments are in disagreement with present results and with those from LEED, as well as, from EAM calculations. We do not understand the reasons for this disagreement but for the notion that X-ray measurements may be very sensitive to the quality of the crystal surface. It should be noted that the differences with the X-ray results are both qualitative and quantitative, beyond the established error bars in the experiments and calculations. Because of the controversy in the experimentally determined multilayer relaxations of Cu(511), we have carried out an extensive analysis of the dependence of the theoretical results on the approximations necessary to produce computational feasibility: choice of pseudopotentials, maximum kinetic energy of the plane-waves (Ecut), the number of layers in the supercell, and the number of points used to sample the surface Brillouin zone. As for the dependence of the results on the type of pseudopotentials and Ecut, we have carried out calculations with three different pseudopotentials (on three sets of codes) to find that the results which lie within 3% of each other. We also find our choice of supercell size to be adequate. There is, however, a strong dependence of the results on the number of BZ points sampled. For the case of Cu(511) this dependence is illustrated in Table VIII. Calculations performed with few points could give erroneous relaxations as signified by the case of d23 in Table VIII. An inward relaxation of 1.8% is found with 4 points, while the converged result is 10.7%. Convergence in the calculated relaxation is reached once the number of points is increased to 24 and beyond. Thus, when comparing results from ab initio calculations, one has to keep these technical points in mind. Unless sufficient checks are made for convergence in the reported values, quantities like equilibrium positions of surface atoms may differ in different calculations and lead to disagreement in the calculated relaxations. It would be worthwhile to clarify whether the differences between our results and those from the FLAPW method for Cu(211) and Cu(331) could be attributed to k-points sampling.

Finally, we come to the comparison of our results for Cu(711) with those from experiments ( we are not aware of any data on Cu(911), so far). The LEED data  for this surface has been very carefully analyzed and compared to existing calculations. Table IX shows that the ab initio results obtained here are in excellent agreement with the data, and that the small differences with the EAM results that the authors had noted, is removed by the present calculations. As in the case of Cu(511), the largest percentage change in the interlayer spacing is not for d12. In this case it is for d34 which separates CC from TC2. As before, there is outward relaxation of the spacing between CC and BNN. The fact that relaxations near CC persist on being strong even as the terrace width increases, is interesting in itself. This particular argument has not been made in any previous theoretical result. Our calculated values for Cu(911) further support this argument as the largest percentage change is found for d45, the interlayer spacing between CC and TC3 (in this case). While these results are intriguing the main outcome of relaxations that ensue when a surface is created is in the actual displacements of the atoms from their bulk terminated positions to the new equilibrium positions. As already stated, these values are summarized in Table V and the related patterns presented in Figs. 2 and 3. Obviously, for stepped surfaces there is a complex rearrangement of most terrace atoms. Our calculations show that despite this complexity, all terrace atoms except for CC move inwards.

## V CONCLUSION

In summary, we have performed a comparative study of multilayer and atomic relaxations of five stepped Cu surfaces which are vicinals of Cu(100) and Cu(111) using ab initio electronic structure calculations based on density functional theory and non-local, norm-conserving pseudopotentials. The set of three of these surfaces: Cu(211), Cu(331), and Cu(511), provides a comparison of structural changes from bulk termination, for vicinals of similar terrace widths but differing local geometry. The other set consisting of Cu(511), Cu(711) and Cu(911) offers a comparative study of relaxation patterns with changing terrace width. In each case we find the relaxation of the step atoms to be pronounced inwards and that of the corner atom to be outwards. The other terrace atoms and their nearest neighbors also undergo relaxations following a complex displacement pattern. Subsequently, the bond lengths between all the surface atoms and their nearest neighbors change from the bulk terminated values while the bond length between CC and BNN atoms experiences an elongation (about 1%) all other surface length shrink anywhere from 1% to 4%. Most of our findings are in agreement with previous calculations which were based on semiempirical model potentials except that we do not find the pattern of inward relaxations of SC, TC1, TC2 etc followed by outward relaxation of CC atoms to continue into inner layers. We also find that the percentage contraction of the spacing between the TC and CC atoms is generally larger than that between SC and the TC atoms. While the actual magnitudes of the changes in the spacing considered here are small, there is a systematic trend in the relaxation pattern which points to significant rearrangements in the electronic charge densities near the SC and CC atoms. By and large our results are in good agreement with available structural data on these surfaces, except for the case of Cu(511) for which we favor the LEED results over those from X-ray scattering meaurements. We believe our results will help settle the issues that have emerged on this particular surface. Our systematic examination of five surfaces, also helps address the question wether the relaxations on Cu(331) are anomalous. The only striking difference between this surface and the others is in the relaxation of TC1 which is very small. Otherwise the relaxation pattern and the changes in bond lengths are similar to those on the other surfaces.

The main message from these observations is that the important quantity to examine is the displacement pattern of the surface atoms as they relax to their equilibrium positions from their bulk terminated configurations. The deeper question, of course, is the nature of the accompanying changes in the surface electronic structure. It will be interesting to examine the characteristics of the local electronic densities of states in the different regions of low symmetry that are present naturally on the stepped surfaces considered here. We leave this as an exercise for the future.

## VI ACKNOWLEDGMENTS

The work of TSR and AK was supported in part by the US National Science Foundation, Grant CHE-9812397 and by the Basic Energy Research Division, Department of Energy, Grant DE-FG03-97ER45650. TSR also acknowledges the Alexander von Humboldt Foundation for the award of a Forschungspreis and thanks her colleagues at the Fritz Haber Institut, Berlin and at the Forschungszentrum, Karlsruhe for their warm hospitality.


text:
# The Nuts and Bolts of Diffraction

## I Hard diffraction: the question

The signature of a diffractive event in $`\overline{p}p`$ collisions is a leading proton or antiproton and/or a rapidity gap defined as a region of pseudorapidity, $`\eta \mathrm{ln}\mathrm{tan}\frac{\theta }{2}`$, devoid of particles. Hard diffraction is a term used to refer to a diffractive process containing a hard partonic scattering (Fig. 1). In deep inelastic scattering (DIS), diffraction is identified by a leading proton in the final state adjacent to a rapidity gap (Fig. 2). The rapidity gap is presumed to be due to the exchange of a Pomeron, whose generic QCD definition is a color-singlet combination of quarks and/or gluons carrying the quantum numbers of the vacuum.  xxx   SD DD DPE FIG. 1: Dijet production in $`\overline{p}p`$ single (a) and double (b) diffraction, and in double Pomeron exchange (c).   xxx   FIG. 2: Diagram for diffractive deep inelastic scattering at HERA.

In addition to its dependence on $`x`$-Bjorken and $`Q^2`$, the diffractive structure function (DSF) of the leading nucleon could also depend on the nucleon’s fractional momentum loss $`\xi `$ and its 4-momentum transfer squared $`t`$. The central question in diffraction is the validity of QCD factorization, i.e. whether hard diffraction processes can be described in terms of parton level cross sections convoluted with a universal DSF.

## II Data: the answer

The question about QCD factorization in diffraction was addressed “head on” by a comparison jj1800  between the DSF measured by CDF in dijet production at the Tevatron and the prediction based on parton densities extracted from diffractive DIS at HERA. The DSF at the Tevatron was found to be suppressed relative to the prediction from HERA by a factor of $`10`$. This result confirmed previous CDF results from diffractive $`W`$ W , dijet jj  and $`b`$-quark b  production at $`\sqrt{s}`$=1800 GeV, and was corroborated by more recent CDF results on diffractive $`J/\psi `$ jpsi  production at 1800 GeV and dijet production at 630 GeV jj630 .

Although factorization breaks down severely between HERA and the Tevatron, it nevertheless holds within the HERA data and within the single-diffractive data at the Tevatron at the same center of mass collision energy. This is demonstrated by the fact that the gluon parton distribution function (PDF) derived from DIS adequately describes diffractive dijet production at HERA jjhera , while at the Tevatron a consistent gluon PDF is obtained from the measured rates of diffractive $`W`$, dijet, $`b`$-quark and $`J/\psi `$ production jpsi . However, factorization was found to break down at the Tevatron between the structure functions measured in single-diffraction and in double-Pomeron exchange (DPE) at $`\sqrt{s}=`$1800 GeV DPE .

## III Soft diffraction: the explanation

The breakdown of QCD factorization observed in hard diffraction is related to the breakdown of Regge factorization responsible for the suppression of soft diffraction cross sections at high energies RR . This may seem paradoxical, but since the rapidity gap formation is a non-perturbative effect it should not come as a surprise. Thus, “the nuts and bolts of diffraction” are contained in soft diffraction processes.

Soft diffraction has been traditionally treated theoretically in the framework of Regge theory. For large rapidity gaps ($`\mathrm{\Delta }\eta 3`$), the cross sections for single and double (central) diffraction can be written as dd

$`{\displaystyle \frac{d^2\sigma _{SD}}{dtd\mathrm{\Delta }\eta }}=\left[{\displaystyle \frac{\beta ^2(t)}{16\pi }}e^{2[\alpha (t)1]\mathrm{\Delta }\eta }\right]\left[\kappa \beta ^2(0)\left({\displaystyle \frac{s^{}}{s_{}}}\right)^{\alpha (0)1}\right]`$ $`{\displaystyle \frac{d^3\sigma _{DD}}{dtd\mathrm{\Delta }\eta d\eta _c}}=\left[{\displaystyle \frac{\kappa \beta ^2(0)}{16\pi }}e^{2[\alpha (t)1]\mathrm{\Delta }\eta }\right]\left[\kappa \beta ^2(0)\left({\displaystyle \frac{s^{}}{s_{}}}\right)^{\alpha (0)1}\right]`$
where $`\alpha (t)=1+ϵ+\alpha ^{}t`$ is the Pomeron trajectory, $`\beta (t)`$ the coupling of the Pomeron to the (anti)proton, and $`\kappa g(t)/\beta (0)`$ the ratio of the triple-Pomeron to the Pomeron-proton couplings. The above two equations are remarkably similar. In each case, there are two terms:

* the first term, which is $`\{\mathrm{exp}[(ϵ+\alpha ^{}t)\mathrm{\Delta }\eta ]\}^2`$ and thus depends on the rapidity gap,
* and the second term, which is $``$exp$`\left[ϵ\mathrm{ln}\left(\frac{s^{}}{s_0}\right)\right]`$ and depends on the pseudorapidity interval $`\mathrm{\Delta }\eta ^{}=\mathrm{ln}\left(\frac{s^{}}{s_0}\right)`$ within which there is particle production.

In the parton model, the second term is interpreted as the sub-energy total cross section, while the first term is the square of the elastic scattering amplitude between the diffractively excited state and the nucleon in SD, or between the two diffractive states in DD. The factor $`\kappa `$ in the second term may then be interpreted as being due to the color matching required for a diffractive rapidity gap to occur. Since the sub-energy cross section is properly normalized, the first factor in the equations may be thought of as the rapidity gap probability and renormalized to unity. A model based on such a renormalization procedure RR ; KGgap  has yielded predictions in excellent agreement with measured SD and DD cross sections, as seen in Figs. 3 and 4.

The renormalized rapidity gap probability is by definition energy independent and thus represents a scaling behaviour. This procedure has the added advantage of preserving unitarity, which otherwise would be violated. Convoluting the gap probability with partonic level cross sections yields hard diffractive cross sections in general agreement with observations, explaining the factorization properties discussed in the previous section RR .

FIG. 3: The $`\overline{p}p`$ total SD cross section exhibits an $`s`$-dependence consistent with the renormalization procedure of Ref. RR , contrary to the $`s^{2ϵ}`$ behaviour expected from Regge theory (figure from Ref. RR ).

FIG. 4: The $`\overline{p}p`$ total DD (central gap) cross section agrees with the prediction of the renormalized rapidity gap model KGgap , contrary to the $`s^{2ϵ}`$ expectation from Regge theory (figure from Ref. dd ).

## IV Multiple rapidity gaps in diffraction

The renormalization method used to calculate the SD and DD cross sections can be extended to multi-gap diffractive events. Below, we outline the procedure for calculating the differential cross section for a 4-gap event (we use rapidity, $`y`$, and pseudorapidity, $`\eta `$, interchangeably).

FIG. 5: Topology of a 4-gap event in pseudorapidity space.

The calculation of the differential cross section is based on the parton-model scattering amplitude:

$$\mathrm{Im}\mathrm{f}(t,\mathrm{\Delta }y)e^{(ϵ+\alpha ^{}t)\mathrm{\Delta }y}$$

For the rapidity regions $`\mathrm{\Delta }y_i^{}`$, where there is particle production, the $`t=0`$ parton model amplitude is used and the sub-energy cross section is $`e^{ϵ\mathrm{\Delta }y^{}}`$. For rapidity gaps, $`\mathrm{\Delta }y`$, which can be considered as resulting from elastic scattering between diffractively excited states, the square of the full parton-model amplitude is used, $`e^{2(ϵ+\alpha ^{}t_i)\mathrm{\Delta }y_i}`$, and the form factor $`\beta ^2(t)`$ is included for a surviving nucleon. The gap probability (product of all rapidity gap terms) is then normalized to unity, and a color matching factor $`\kappa `$ is included for each gap.

For the 4-gap example of Fig. 5, which has 10 independent variables, $`V_i`$ (shown below the figure), we have:

* $`\frac{d^{10}\sigma }{\mathrm{\Pi }_{i=1}^{10}dV_i}=P_{gap}\times \sigma (\mathrm{sub}\mathrm{energy})`$
* $`\sigma (\mathrm{sub}\mathrm{energy})=\kappa ^4\left[\beta ^2(0)e^{ϵ\mathrm{\Delta }y^{}}\right]`$ ($`\mathrm{\Delta }y^{}=_{i=1}^3\mathrm{\Delta }y_i^{}`$)
* $`P_{gap}=N_{gap}\times \mathrm{\Pi }_{i=1}^4\left[e^{(ϵ+\alpha ^{}t_i)\mathrm{\Delta }y_i}\right]^2\times [\beta (t_1)\beta (t_4)]^2=N_{gap}e^{2ϵ\mathrm{\Delta }y}f(V_i)|_{i=1}^{10}`$ ($`\mathrm{\Delta }y=_{i=1}^4\mathrm{\Delta }y_i`$)
where $`N_{gap}`$ is the factor that normalizes $`P_{gap}`$ over all phase space to unity.

The last equation shows that the renormalization factor of the gap probability depends only on $`s`$ (since $`\mathrm{\Delta }y_{max}=\mathrm{ln}s`$) and not on the number of diffractive gaps. Thus, the ratio of DPE to SD cross sections is expected to be $`\kappa `$, with no additional energy dependent suppression. This can be tested at the Tevatron, where one may also study events with a central rapidity gap within single-diffractive clusters. The study of events with more than two gaps will have to await the commissioning of the LHC.