text:
# Optimal basis set for electronic structure calculations in periodic systems

## I Introduction

First-principles methods have become a widespread tool for studying the microscopic aspects of a wide class of condensed systems, including their crystalline phases. No empirical modelling of the interatomic interactions is required in first-principles simulations, since the electronic ground state is evaluated at each atomic configuration. In a large set of systems, mostly in metals, but also in some molecular crystals like hydrogen, a very fine sampling of the Brillouin Zone (BZ) is required to provide a sufficiently accurate description of the electronic ground state. Since the computational load for refining BZ sums grows linearly with the number of sampling points (k-points), special k-point grids that minimize sampling errors, as well as interpolation techniques (the so called tetrahedron method) have been introduced in the past. However, little attention has been paid to the fact that Bloch states at nearby k-points are not totally uncorrelated. In fact, once the Bloch states at a given k-point are known, those at a closeby k-point can be easily estimated by perturbation theory, within the so-called $`kp`$ approximation. This extremely helpful piece of information has been little exploited so far in standard first-principles approaches, where for each k-point a full solution of the Kohn-Sham equations is instead calculated independently, e.g. through diagonalization of the Kohn-Sham hamiltonian on a finite basis set (typically plane waves in the case of crystalline systems). An approximate technique based on the exact solution of the Kohn-Sham hamiltonian at a few k-points has been proposed, but the analysis of the errors associated with such method, as well as their correction, has proved rather difficult.

In this work we propose a simple and efficient method for solving the Kohn-Sham equations for a very large number of k-points in the BZ, with an effort that in most cases is comparable to that of solving the problem for just a few k-points. The method is based on a rigorous formulation of the ideas sketched above, and relies on the construction of a single, optimal basis set for the expansion of the Bloch orbitals, which is common to all k-points. The application of this method to the study of the compressed phases of molecular hydrogen has been recently reported . In Section II we formulate the method. In Section III we analyse the reliability of the method on a test case (16 Na atoms in the bcc lattice), and in Section IV we analyse the computational cost. Section V contains the conclusions.

## II The method

In first-principles simulations the electronic ground state is calculated within the Kohn-Sham, self-consistent one-electron formalism of density functional theory . Translational invariance is introduced by averaging the electronic density and the total energy of the system on a finite number of k-points in the BZ:

$$\rho (𝐫)=\underset{𝐤}{\overset{BZ}{}}\omega _𝐤\underset{i}{}f_i(𝐤)\psi _i^𝐤(𝐫)^2,$$
(1)
$$E[\rho ]=T+E_{\mathrm{ext}}[\rho ]+E_\mathrm{H}[\rho ]+E\mathrm{xc}[\rho ]$$
(2)
with

$$T=\underset{𝐤}{\overset{BZ}{}}\omega _𝐤\underset{i}{}f_i(𝐤)\psi _i^𝐤(𝐫)\left(\frac{^2}{2}\right)\psi _i^𝐤(𝐫)𝑑𝐫,$$
(3)
where $`i`$ labels bands, k a k-vector in the BZ, $`\omega _𝐤`$ its weight, $`f_i(𝐤)`$ are occupation numbers, and $`\psi _i^𝐤(𝐫)`$ are the Kohn-Sham orbitals. The last three terms in (2) are the external, Hartree, and exchange-correlation contributions to the energy, and depend on the BZ averaged electronic density (1). Atomic units are used throughout this work.

Variational equations for the periodic part of the Kohn-Sham orbitals $`u_i^𝐤(𝐫)=\mathrm{exp}(i𝐤𝐫)\psi _i^𝐤(𝐫)`$ are readily obtained from expression (2):

$$\left(\frac{^2}{2}i𝐤+\frac{k^2}{2}+V[\rho (𝐫)]\right)u_i^𝐤(𝐫)=ϵ_i^𝐤u_i^𝐤(𝐫),$$
(4)
where $`V[\rho ]`$ is the Kohn-Sham potential, containing external, Hartree and exchange-correlation contributions. Normally, equations (4) are solved independently for each k-point, a new density is constructed via (1), and the cycle is iterated to self-consistency. Equations (4) are often solved by expanding the orbitals corresponding to the various k-points in a basis set of plane waves, which remains unchanged when the atomic configuration is changed. The choice of a plane-wave basis set offers some practical advantages, including the convenience of calculating a number of quantities (matrix elements, form factors, etc) a single time throughout the calculation. However, working with a plane-wave basis set becomes particularly demanding from the computational point of view when a large supercell and/or a large number of $`k`$-points has to be considered, a situation which is often encountered when dealing with metals. Moreover, a plane-wave basis set cannot straightforwardly exploit the similarity of Bloch states at closeby k-points.

In the following we will show that relaxing the requirement of a “simple” but inmutable basis set in favor of a variational search for the “best” basis set at each atomic configuration may improve substantially the efficiency of BZ summations. To this aim, let us expand the Kohn-Sham orbitals on a generic basis set $`\varphi _i(i=1,\mathrm{},N)`$ of size $`N`$:

$$u_i^𝐤(𝐫)=\underset{j=1}{\overset{N}{}}a_{ij}(𝐤)\varphi _j(𝐫).$$
(5)
In the case of standard calculations $`\varphi _i`$ would be a plane wave – $`\varphi _i(𝐫)\mathrm{exp}(i𝐠𝐫)`$, being g a reciprocal lattice vector – and solving the Kohn-Sham equations would amount to minimize the energy functional (2) with respect to the expansion coefficients $`a_{ij}`$, without modifying the basis set. In the present approach we minimize the energy functional (2) with respect to both the expansion coefficients $`\{a_{ij}\}`$ and the basis set functions $`\{\varphi _i\}`$. In other words, we look simultaneously for the ground state orbitals and for basis set that better describes them. For the sake of clarity, we note that this would be a redundant request for calculations with a single k-point, where a single set of Kohn-Sham equations have to be solved. The request of an optimal basis set becomes non redundant when the basis set is used to describe Kohn-Sham equations at many k-points.

Explicit minimization of (2) with respect to, simultaneously, $`a_{ij}`$ and $`\varphi _i`$ leads to the following system of equations. For the expansion coefficients we have the usual eigenvalue equations:

$$\underset{j=1}{\overset{N}{}}\left(\lambda _{ij}^0+\frac{k^2}{2}\delta _{ij}+𝐤𝐩_{ij}\right)a_{jl}(𝐤)=ϵ_i^𝐤a_{il}(𝐤),$$
(6)
with

$$\lambda _{ij}^0=\varphi _i\left|\left(\frac{^2}{2}+V_{\mathrm{KS}}[\rho ]\right)\right|\varphi _j$$
(7)
and

$$𝐩_{ij}=i\varphi _i||\varphi _j.$$
(8)
For the basis-set functions we have

$`{\displaystyle \frac{\delta E[\varphi ,𝐚]}{\delta \varphi _i^{}(𝐫)}}`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}B_{ij}^{(0)}\left[{\displaystyle \frac{^2}{2}}+v_{\mathrm{KS}}[\rho ]\right]\varphi _j^{}(𝐫)+`$ (9)
$`+`$ $`{\displaystyle \underset{j=1}{\overset{N}{}}}\left\{i𝐁_{ij}^{(1)}+B_i^{(2)}\delta _{ij}\right\}\varphi _j^{}(𝐫)+`$ (10)
$`+`$ $`{\displaystyle \frac{\delta E_{NL}[\varphi ,𝐛]}{\delta \varphi _i^{}(𝐫)}}={\displaystyle \underset{l=1}{\overset{N}{}}}\mathrm{\Lambda }_{il}\varphi _l(𝐫),`$ (11)
with

$`B_{ij}^{(0)}`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}\omega _𝐤b_{ij}(𝐤)`$ (12)
$`𝐁_{ij}^{(1)}`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}\omega _𝐤𝐤b_{ij}(𝐤)`$ (13)
$`B_i^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{𝐤}{}}\omega _𝐤k^2b_{ii}(𝐤)`$ (14)
and

$$b_{jl}(𝐤)=\underset{i=1}{\overset{N}{}}f_i(𝐤)a_{ij}^{}(𝐤)a_{il}(𝐤).$$
(15)
The $`\mathrm{\Lambda }`$ coefficients in (11) are Lagrange multipliers introduced to enforce the orthonomalization of the basis functions . The electron density in (11) is constructed from the Kohn-Sham orbitals (5) using (15), as

$$\rho (𝐫)=\underset{j=1}{\overset{N}{}}\underset{l=1}{\overset{N}{}}\varphi _j^{}(𝐫)\beta _{jl}\varphi _l(𝐫),$$
(16)
where

$$\beta _{jl}=\underset{𝐤}{}\omega _𝐤b_{jl}(𝐤).$$
(17)
More detailed expressions for the above equations as well as for atomic forces and stress are given in the Appendix, for the case of pseudopotential plane wave (PPW) calculations.

The problem of simultaneously solving Kohn-Sham equations for $`N_k`$ k-points independently, as done in standard electronic structure calculations, is here replaced by the problem of solving a single set of Kohn-Sham-like equations, eqs. (11), supplemented with $`N_k`$ diagonalizations of the $`N\times N`$ matrix in (6). As we will show in the next Section, the number of basis functions $`N`$ required for the optimal basis set $`\{\varphi _i\}`$ can be chosen to be much smaller than that in the basis set used to describe each single $`\varphi _i`$, e. g. the number of plane wave components in PPW calculations. This implies that, if $`N_k`$ is large, the computational effort required by this procedure is much smaller than the effort required to solve the standard set of Kohn-Sham equations at all k-points. The main limitation of the method resides, then, in the size $`N`$ of the optimal basis set. If $`N`$ turns out to be comparable to the original number of basis functions, then there is no point on pursuing this strategy.

While in standard calculations the number of Kohn-Sham equations is roughly (exactly for insulators) half the number of electrons, in our approach $`N`$ has to be sufficiently large to guarantee the proper description of the Kohn-Sham orbitals (5) at all k-points. A full analysis of the convergence of the method as a function of $`N`$ will be presented later, for a specific case. There are however reasons to expect that $`N`$ need not be too large. In fact, a reasonable approximation to the exact optimal basis set is given by the Kohn-Sham orbitals at one selected k-point, say at zone center ($`\mathrm{\Gamma }`$). We call this basis set: $`\{\varphi _i^\mathrm{\Gamma }\}`$. Since the optimal basis set is defined as the best basis set in a variational fashion,$`\{\varphi _i^\mathrm{\Gamma }\}`$ will necessarily provide a worse description of the Kohn-Sham orbitals (except, of course, at $`\mathrm{\Gamma }`$) than $`\{\varphi _i\}`$. The convergence properties of $`\{\varphi _i\}`$ will thus be better than those of $`\{\varphi _i^\mathrm{\Gamma }\}`$. If the BZ of the system of interest is sufficiently small (which is equivalent to ask that the unit cell of the calculation is large), then $`kp`$ perturbation theory suggests that Kohn-Sham orbitals at any k-point other than $`\mathrm{\Gamma }`$ may be expressed in terms of a few number of energetically close $`\mathrm{\Gamma }`$-point orbitals $`\{\varphi _i^\mathrm{\Gamma }\}`$ . In other words, only a few excited $`\mathrm{\Gamma }`$-point orbitals are required for a reasonable description of the orbitals at all other k-points. This property has been already noted and used in the context of first-principles electronic structure calculations. Our approach differs from that of Robertson and Payne  in the fact that their basis set is given by the Kohn-Sham orbitals at $`\mathrm{\Gamma }`$ (or at a small set of k-points), while our method includes a variational search for the best basis set. A comparison of the convergence properties of the two approaches is given in the next Section.

We also observe that our method has only one variational parameter, namely the size $`N`$ of the basis set. This implies that systematic improvements of the accuracy can be achieved by simply increasing $`N`$.

Finally, we notice that eqs. (11), like ordinary Kohn-Sham equations, have to be solved self-consistently, since the Kohn-Sham potential $`v_{KS}`$ depends on the density, and thus on wavefunctions. However, unlike ordinary Kohn-Sham equations, eqs. (11) depend also on the expansion coefficients $`a_{ij}`$. The ordinary self-consistent iterative procedure has, consequently, been modified in the following way:

1. an initial guess of the density is provided, e.g. as a sum of atomic densities; the $`a_{ij}`$ coefficients are initially set to $`\delta _{ij}`$,
2. eqs. (11) are solved,
3. the orbitals (which are no longer Kohn-Sham orbitals) are used to construct $`\lambda _{ij}^0`$ and $`𝐩_{ij}`$,
4. eq. (6) is solved for each k-point, via direct diagonalization,
5. a new density is constructed,
6. back to 2.

Notice that this procedure does not involve a major implementation effort, but only a few modifications to the usual self-consistency methods.

## III Convergence of the optimal basis set

In order to analyse the convergence properties of the method we have chosen the case of Na metal, which is properly described with a soft, local pseudopotential. The Ceperley-Alder local density functional was adopted for the exchange-correlation term . Sixteen Na atoms in the bcc structure have been placed in a simple cubic supercell ($`2\times 2\times 2`$ bcc conventional cells), and the full Brillouin zone of the supercell has been sampled with 108 k-points. A plane wave energy cutoff of 15 Ry was chosen, amounting to 1596 plane waves in the expansion of the optimal basis set. The convergence of the method has been analysed by performing calculations with $`N=19`$, 32, 51, 81, and 1596 . Since our aim is at checking the convergence properties of the method, and not the overall accuracy as determined by the other approximations (LDA, k-point grid, etc.), we will consider as converged the values obtained with $`N=1596`$ and $`N_k=108`$. As a comparison, we have also determined the convergence properties of a calculation where the basis set, instead of being “optimal”, is chosen to be the set of Kohn-Sham orbitals resulting from a self-consistent calculation using the $`\mathrm{\Gamma }`$-point only. This correponds to the approach of Ref. .

In Figs. 1 to 3 we show the convergence of the different physical properties as a function of $`N`$, namely the total energy, stress, and electronic density of states. The total energy (Fig. 1) shows an extremely good convergence even for the smallest value of $`N`$. The accuracy is better than 0.5 % for $`N=19`$, and better than 0.1 % for $`N=32`$ states. The same behavior is observed in the force constant calculated for a distortion along an optical phonon: the force is already converged to better than 0.1 % using 32 states in the expansion. Thus, the errors introduced in energies and forces due to truncation of the expansion are much smaller than, e.g., the typical accuracy of the usual density functionals used, namely the LDA or the GGA. Stresses also show good convergence properties, as reported in Fig. 2. Variations of stress ($`\mathrm{\Delta }\sigma `$) with respect to $`N`$ should in this case be compared with the value of the bulk modulus of Na ($`B=64`$ kbar). In fact, an error $`\mathrm{\Delta }\sigma `$ in the stress yields an error in the equilibrium lattice spacing that can be estimated as $`\mathrm{\Delta }a/a\mathrm{\Delta }\sigma /3B`$. Thus, even with $`N=19`$ the equilibrium lattice spacing is converged within less than 1 % .

The electronic density of states (EDOS) is also sufficiently accurate for the relevant portion of the spectrum, namely the one below the Fermi energy ($`E_F`$), using $`N=32`$. In contrast, $`N=19`$ appears to give an inaccurate EDOS, particularly in the vicinity of $`E_F`$, as expected in a $`kp`$ picture. In fact, as discussed in Section II, a small basis set is more likely to affect the accuracy in the description of states with higher energy, i.e. closer to $`E_F`$, due to the insufficient number of high energy states in the optimal basis set.

We now compare the results obtained with our method against those obtained with the method of Ref. . Our energies are always lower, for a given basis set size $`N`$, than energies calculated using the method of Ref. . This is not unexpected, since the optimal basis set is variationally the best basis set. Moreover, we find that the method of Ref.  provides results with an accuracy similar to ours only when the size of the basis set is 30-50 % larger than ours. The most striking difference, however, is in the convergence of the stress tensor. While our method converges very smoothly to the exact value, the method of Ref.  shows an erratic behavior, that clearly prevents every attempt to determine the exact value on the basis of size scaling. We believe that such a marked difference in the convergence properties has to be traced to the variational nature of our approach. Therefore, our method not only improves on that of Ref.  in the rate of convergence (which, in the case of energy, results in a 30-50 % saving in computational time, as shown in the next Section), but shows a much smoother behavior as a function of $`N`$, resulting in a more accurate estimate of the error introduced by the finite size of the basis set.

## IV Estimate of the computational load

The computational effort in a standard plane wave calculation scales as $`N_kN_{bands}N_{pw}^2`$, where $`N_{bands}`$ is the number of occupied bands (equal or slightly larger than half the number electrons) and $`N_{pw}`$ is the size of the plane wave basis set. This estimate is based on the assumption that standard diagonalization methods are used, such as iterative diagonalization, where the effort of determining the first $`N_{bands}`$ eigenvalues of a $`N_{pw}\times N_{pw}`$ matrix scales like $`N_{bands}N_{pw}^2`$. Since every k-point requires a separate diagonalization, the total effort scales as $`N_kN_{bands}N_{pw}^2`$.

In the present approach, if we assume that the effort of diagonalizing $`N_k`$ times the $`N\times N`$ matrix (6) is negligible with respect to the effort required to solve eqs. (11), then the load scales as $`NN_{pw}^2`$. Thus, a gain of a factor $`N_kN_{bands}/N`$ is expected with our method. Since in typical applications $`N2÷4N_{bands}`$, and $`N_k10^2`$, a computational gain of more than one order of magnitude is obtained. The neglected term becomes increasingly important as the size of the system ($`N`$) increases, scaling as $`N_k\times N^3`$. Therefore, it becomes dominant for sufficiently large $`N`$, thus discouraging the use of the present methodology in the case of large unit cells and BZs sampled in a few k-points.

The above estimate is based on the assumption that solving eqs. (11) is as time consuming as solving the standard Kohn-Sham equations, except for the fact that the requested number of states is larger. Our experience on Na (this work) and molecular hydrogen  shows that the computational overload of solving eqs. (11) basically scales linearly with the number of states; the extra effort required for solving a single state is negligible when compared to that of solving one of the usual Kohn-Sham states. Nevertheless, we are currently not in a position to generalize this conclusion to a generic system.

## V Summary and conclusions

We have presented and analysed the convergence properties of a method for performing electronic structure calculations and first-principles molecular dynamics simulations for systems that need a very fine sampling of the electronic Brillouin zone. The method is based on the construction of a basis set which optimally describes the Kohn-Sham orbitals at all k-points. The rapid convergence of the method on the size $`N`$ of the basis set is connected with the properties of the $`kp`$ approximation and has been fully analysed in the case of a sixteen Na atoms supercell. We have shown that $`N`$ can be kept reasonably small, of the order of a few times the number of occupied bands (2 to 4 times, depending on the accuracy desired), resulting in a gain in computational load of about one order of magnitude, with respect to standard methods.

Even if our analysis was mainly done in connection with PPW calculations, the same philosophy can be adopted in other electronic structure schemes like the LAPW method , provided that the number of orbitals required in the optimal basis set remains significantly smaller than the number of basis functions in the original basis set.

In conclusion, this technique constitutes a useful tool for studying the electronic, structural and dynamical properties of metals and more generally for any system whose proper description requires a fine sampling of the BZ, like, for example, molecular hydrogen.

## VI Acknowledgements

SS acknowledges partial support from INFM through “Iniziativa Calcolo Parallelo”, and from MURST through COFIN 1999. JK ackowledges the hospitality of the Universidad Nacional de San Luis, Argentina, where part of this work was carried out.

## VII Appendix

In this appendix we give some useful expressions for the total energy, atomic forces and stress, in the case of pseudopotential plane wave calculations.

The different contributions to the total energy (2) can be divided into the class of terms that only depend on the density (the hartree, the exchange-correlation, and the local part of the pseudopotential), and the class of terms that explicitly depend on the wavefunctions (the kinetic and the nonlocal pseudopotential contributions). While the former can be calculated in the usual way once the density is known (see below for an efficient way to evaluate density), the latter can be expressed more conveniently as follows. For the kinetic energy we have

$`T`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}\omega _𝐤{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{j=1}{\overset{N}{}}}b_{ij}^𝐤\left\{t_{ij}^0[\varphi ]+𝐤𝐩_{ij}+{\displaystyle \frac{k^2}{2}}\delta _{ij}\right\}=`$ (18)
$`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}\left\{B_{ii}^{(0)}t_{ii}^0[\varphi ]+B_i^{(2)}\right\}+`$ (19)
$`+`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \underset{ji}{\overset{N}{}}}\left\{B_{ij}^{(0)}t_{ij}^0[\varphi ]i𝐁_{ij}^{(1)}𝐩_{ij}\right\},`$ (20)
with

$$t_{ij}^0[\varphi ]=<\varphi _i\frac{^2}{2}\varphi _j>,$$
(21)
and for the non-local contribution we have

$`E_{NL}^l`$ $`=`$ $`{\displaystyle \underset{s=1}{\overset{\sigma }{}}}{\displaystyle \underset{m=l}{\overset{l}{}}}\alpha _{lm}^s{\displaystyle \underset{ij}{}}{\displaystyle \underset{𝐤}{}}\omega _𝐤b_{ij}(𝐤)\times `$ (22)
$`\times `$ $`{\displaystyle \underset{I=1}{\overset{N_s}{}}}F_{Iilm}^{}(𝐤)F_{Ijlm}(𝐤),`$ (23)
where, in a plane waves basis set, $`F_{Iilm}(𝐤)`$ assumes the following form:

$$F_{Ijlm}(𝐤)=\underset{𝐠}{}w_{lm}^s(𝐠+𝐤)\mathrm{exp}\left(i𝐠𝐑_I\right)\varphi _i(𝐠),$$
(24)
being $`\{𝐠\}`$ the reciprocal lattice vectors and $`𝐑_I`$ the coordinates of atom $`I`$ belonging to species $`s`$. The coefficient $`w_{lm}^s`$ is the projection of the $`𝐠`$-th plane wave component expanded around the k-vector $`𝐤`$, onto the $`lm`$-th subspace according to the Kleinman-Bylander prescription:

$$w_{lm}^s(𝐠+𝐤)=<\delta v_l^s\phi _{lm}^sY_{lm}𝐠+𝐤>,$$
(25)
with $`\delta v_l^s`$ the $`l`$-th angular momentum pseudopotential component for species $`s`$, $`\phi _{lm}^s`$ the corresponding pseudoatomic orbital, and $`Y_{lm}`$ the spherical harmonic functions. Equation (23) also contains a normalization factor, which is given by the expression $`\alpha _{lm}^s=<\delta v_l^s\phi _{lm}^s\phi _{lm}^s>^1`$. $`N_s`$ is the number of atoms of species $`s`$, and $`\sigma `$ is the number of species described through angular-dependent pseudopotentials.

It is interesting to remark that, since now the factors $`w_{lm}^s`$ depend explicitly on the k-vector, it is no longer possible to perform the BZ summation on the coefficients $`b_{ij}(𝐤)`$ once and for all, as it is done for the kinetic and local potential terms. Now this sum has to be carried out for each PW component $`𝐠`$, and this is likely to result in a costly computational scheme. So far we have not studied this issue in detail, because this method has only been applied to hydrogen and sodium, where a local pseudopotential description is feasible, but we hope that alternatives can be found so that the non-local part does not become the bottleneck of the calculation.

The calculation of the density is a subtle issue, because expression (16) involves a double sum of matrix products over all states, occupied and empty. Even if it has been suggested  that this operation is computationally convenient, we have found it more efficient to first transform the orbitals at $`\mathrm{\Gamma }`$ by multiplying them with the rotation matrices $`a_{ij}(𝐤)`$, and then carry out a single summation over the occupied states.

The last term in (11) is easier to compute in reciprocal space, where it assumes the following expression:

$`{\displaystyle \frac{\delta E_{NL}[\varphi ,𝐛]}{\delta \varphi _i^{}(𝐠)}}`$ $`=`$ $`{\displaystyle \underset{s=1}{\overset{\sigma }{}}}{\displaystyle \underset{lm}{}}\alpha _{lm}^s{\displaystyle \underset{j}{}}{\displaystyle \underset{𝐤}{}}\omega _𝐤b_{ij}(𝐤)\times `$ (26)
$`\times `$ $`w_{lm}^s(𝐠+𝐤){\displaystyle \underset{I=1}{\overset{N_s}{}}}\mathrm{exp}\left(i𝐠𝐑_I\right)F_{Ijlm}(𝐤).`$ (27)
The forces on the nuclei are trivially unchanged from the original ones, as long as local pseudopotentials are used. This is because the electronic kinetic term does not depend explicitly on the nuclear positions. In the case of angular-dependent pseudopotentials the contribution from the last term in (23) to the force on atom $`I`$ belonging to species $`s`$, is modified from the standard result in the following way:

$`{\displaystyle \frac{E_{NL}}{𝐑_I}}`$ $`=`$ $`2{\displaystyle \underset{lm}{}}\alpha _{lm}^s{\displaystyle \underset{ij}{}}{\displaystyle \underset{𝐤}{}}\omega _𝐤\times `$ (28)
$`\times `$ $`\mathrm{Re}\left(b_{ij}(𝐤){\displaystyle \frac{F_{Iilm}^{}(𝐤)}{𝐑_I}}F_{Ijlm}(𝐤)\right).`$ (29)
The calculation of the stress matrix is also modified with respect to standard calculations only for the kinetic (kin) and nonlocal pseudopotential (NL) contributions. The resulting expression for the kinetic contribution is:

$`\sigma _{\alpha \beta }^{kin}`$ $`=`$ $`{\displaystyle \underset{ij}{}}B_{ij}^{(0)}t_{ij\alpha \beta }^0+{\displaystyle \underset{ij}{}}\left(B_{ij\beta }^{(1)}p_{ij\alpha }+B_{ij\alpha }^{(1)}p_{ij\beta }\right)+`$ (30)
$`+`$ $`{\displaystyle \underset{i}{}}B_{i\alpha \beta }^{(2)},`$ (31)
with

$$t_{ij\alpha \beta }^0=<\varphi _i\frac{^2}{x_\alpha x_\beta }\varphi _j>$$
(32)
and

$`B_{ij\alpha }^{(1)}`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}\omega _𝐤k_\alpha b_{ij}(𝐤)`$ (33)
$`B_{i\alpha \beta }^{(2)}`$ $`=`$ $`{\displaystyle \underset{𝐤}{}}\omega _𝐤k_\alpha k_\beta b_{ii}(𝐤).`$ (34)
The nonlocal contribution is instead:

$`\sigma _{\alpha \beta }^{NL}`$ $`=`$ $`2{\displaystyle \underset{s=1}{\overset{\sigma }{}}}{\displaystyle \underset{lm}{}}\alpha _{lm}^s{\displaystyle \underset{ij}{}}{\displaystyle \underset{𝐤}{}}\omega _𝐤b_{ij}(𝐤)\times `$ (35)
$`\times `$ $`{\displaystyle \underset{I=1}{\overset{N_s}{}}}\mathrm{Re}\left({\displaystyle \frac{F_{Iilm}^{}(𝐤)}{𝐡_{\alpha \beta }}}F_{Ijlm}(𝐤)\right),`$ (36)
where $`𝐡_{\alpha \beta }`$ is the matrix of the three primitive Bravais vectors , and the partial derivative in (36) is given in Ref. .


format: text

text:
# Width and Magnetic Field Dependence of Transition Temperature in Ultranarrow Superconducting Wires

## I Introduction

Recent experiments on disordered superconducting wires have led to renewed interest in the interplay of interaction and disorder in low-dimensional systems. These experiments have been made possible by the development of fabrication techniques using either stencils or carbon nanotubes, and involve wires of width as low as $`50\AA `$ and thickness as low as $`10\AA `$. An obvious question to ask in such systems is how the transition temperature is affected by wire width and thickness. Decreasing the thickness of the wire increases the normal state resistance per square which is a measure of the amount of disorder in the system. This increase in disorder reduces the transition temperature by causing the electrons to move diffusively rather than ballistically, as discussed below. Decreasing the width of the wire leads to a crossover from two- to one-dimensional behavior, where increased fluctuation effects further reduce the transition temperature. Ultimately the combined disorder and fluctuation effects completely suppress superconductivity, and a superconductor-insulator transition results. The effect of dimensionality on this quantum phase transition is obviously very interesting from a theoretical point of view. As a first step to understanding this transition, we will compare predictions of transition temperature to experimental results. Unfortunately there is very little data available for comparison because each data point requires the fabrication of a new wire which is very costly in terms of experimental effort. We therefore suggest the application of a perpendicular magnetic field as a way to generate more data points for each wire and make predictions for this effect.

The most basic model used consists of electrons moving in a disordered background interacting via the Coulomb repulsion and a featureless phonon-mediated attraction. Simple as this model is, significant difficulties arise even when calculating the suppression of transition temperature in two-dimensional films. These difficulties will be compounded when looking at quasi one-dimensional wires, which are in a crossover regime, so it is most important that they are fully understood in the simpler two-dimensional case. The basic physics involved is not in dispute, and can be understood from the BCS-MacMillan formula for the transition temperature, $`T_c`$,

$$T_c=1.13\omega _D\mathrm{exp}\left(\frac{1}{N(0)[\gamma \mu ^{}]}\right),$$
(1)
where $`\omega _D`$ is the Debye frequency (a typical phonon frequency), $`N(0)`$ is the electronic density of states per spin at the Fermi surface, $`\gamma `$ is the featureless attraction, and $`\mu ^{}`$ is the Coulomb pseudopotential, a measure of the Coulomb repulsion. The disorder causes the electrons to move diffusively rather than ballistically making them less efficient at screening the Coulomb repulsion. This increases $`\mu ^{}`$ and decreases $`N(0)`$, both effects leading to a decrease in $`T_c`$. We shall assume that the disorder does not affect phonon properties such as $`\omega _D`$ and $`\gamma `$, an assumption that is backed up by detailed calculation on at least one model of electron-phonon interaction. Other authors have questioned this assumption, and we intend to address this point in future work.

The first difficulty which occurs in two dimensions relates to the low-momentum singularity of the disorder-screened Coulomb potential, which is given by

$$V_C(q,\omega _m)=\frac{4\pi e^2}{q^2}\frac{Dq^2+|\omega _m|}{Dq_{TF}^2+Dq^2+|\omega _m|}\frac{1}{2N(0)}\frac{Dq^2+|\omega _m|}{Dq^2},$$
(2)
where $`q_{TF}`$ is the Thomas-Fermi screening wavenumber defined by $`q_{TF}^2=8\pi e^2N(0)`$, $`D=v_F^2\tau /2`$ is the diffusion constant, $`v_F`$ is the Fermi velocity, and $`\tau `$ is the elastic scattering time. This singularity leads to a strong suppression of the one-particle density of states which is seen experimentally in tunneling measurements. Naively inserting this density of states into Eqn.  predicts a large suppression of $`T_c`$ which is not seen experimentally. This prediction is, in fact, incorrect due to a cancellation between the various diagrammatic contributions in perturbation theory. The apparently mysterious cancellation turns out to be due to gauge invariance and is thus expected to be a very robust phenomenon. The end result is that the Coulomb interaction is effectively featureless with strength $`N(0)V_C=1/2`$. The fact that we can legitimately use a featureless interaction makes calculations much simpler. Our first task is therefore to check that this replacement is still valid in quasi-1D systems.

The other difficulty is that perturbation theory cannot take us to the superconductor-insulator transition as it is only valid for weak disorder, and ad hoc attempts to extend its range of validity lead to unphysical reentrance problems. Oreg and Finkel’stein overcame this problem by developing a resummation technique, which has recently been shown to be related to strong-coupling theory. We reexamine their predictions for wires, and find, rather surprisingly, that the problem is extremely sensitive to the boundary conditions used. This unexpected sensitivity leads us to check all details and approximations used in the resummation technique. Upon satisfactorily completing the necessary checks, we compare theory and experiment. Although good agreement is obtained when we focus on the suppression of $`T_c`$ due to increasing disorder, it is less good when we focus on the additional suppression due to decreasing wire width. This is largely due to the limited amount of experimental data (only about four data points for each of four wire widths). We propose the application of a perpendicular magnetic field as a means of generating more data with a minimum of effort. We extend the resummation technique to include a magnetic field, and make predictions for transition temperature as a function of magnetic field for different wire widths and thicknesses. The most interesting feature which emerges is a crossover from two- to one-dimensional behavior which can be tuned by either decreasing wire width or increasing the resistance per square. We hope that the appropriate experiments will be carried out in the near future.

The remainder of the paper is organised in the following manner. In section (II) we consider first-order perturbation theory in the absence of a magnetic field. This is the simplest approach possible, and we will use it as a testing ground for the validity of the approximations used. Specifically we will show that it is legitimate to use a featureless Coulomb interaction and ignore self-energy effects, but that the details of boundary conditions are important. In section (III) we make predictions for the transition temperature using the resummation technique, all the necessary approximations having been justified in section (II), and then compare theory and experiment. In section (IV) we consider the effect of an applied magnetic field on the transition temperature. The dimensions of the wires are such that we must consider both orbital and Zeeman (spin) effects. As yet there is no experimental data available for comparison to the theory.

## II First Order Perturbation Theory

We will now calculate the transition temperature, $`T_c`$, which we identify as the temperature at which the pair propagator, $`L(q,\mathrm{\Omega }=0)`$, diverges. The latter is defined in Fig. 1(a) in terms of the BCS interaction, $`\gamma `$, and the pair polarization function, $`P(q,0)`$, by

$$L(q,0)^1=\gamma ^1P(q,0).$$
(3)
Let us first consider mean field theory, which corresponds to zeroth order in disorder perturbation theory, for which the expansion parameter is $`1/k_F\mathrm{}`$, where $`k_F`$ is the Fermi wavenumber and $`\mathrm{}=v_F\tau `$ is the elastic mean free path. The mean field pair polarization function is shown in Fig. 1(b), and leads to the pair propagator

$$L_0(q,0)^1=N(0)\left[\mathrm{ln}\left(\frac{T}{T_{c0}}\right)+\psi \left(\frac{1}{2}+\frac{Dq^2}{4\pi T}\right)\psi \left(\frac{1}{2}\right)\right].$$
(4)
The first-order perturbative corrections to the pair polarization function, $`\delta P(q,0)`$, are shown in Fig. 1(c), and lead to a suppression in transition temperature

$$\mathrm{ln}\left(\frac{T_c}{T_{c0}}\right)=\frac{\delta P(q,0)}{N(0)}.$$
(5)
Note that we evaluate the pair propagator at non-zero external momentum, $`q`$, since we will eventually consider the effect of a magnetic field.

The five diagrams of Fig. 1(c) can be understood in physical terms: diagrams 1 and 2 are self-energy diagrams and give corrections to the electronic density of states, $`N(0)`$; diagrams 3 and 4 are pseudopotential diagrams and give corrections to the effective Coulomb repulsion, $`\mu ^{}`$; diagram 5 involves the effect of Coulomb interaction on superconducting fluctuations. We will first perform the calculations exactly, and demonstrate the cancellation of the low-momentum singularity in the Coulomb interaction which occurs between diagrams 1–4 and 5. We show that this is equivalent to considering only diagrams 1–4 with a featureless Coulomb interaction of strength $`N(0)V_C=1/2`$. In other words, the only purpose of diagram 5 is to remove the low-momentum singularity, as required by gauge invariance. Finally we show that the self-energy diagrams 1 and 2 can be ignored as they give a much smaller contribution than diagrams 3 and 4. The contributions from the five diagrams are as follows:

$`\delta P_1`$ $`=`$ $`4\pi N(0)T{\displaystyle \underset{ϵ_l}{}}T{\displaystyle \underset{\omega _m}{}}{\displaystyle \frac{1}{aL}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{Dq^2+Dq^2+3|ϵ_l|+|ϵ_l+\omega _m|}{(Dq^2+2|ϵ_l|)^2(Dq^2+|ϵ_l|+|ϵ_l+\omega _m|)^2}}V_C(q^{},|\omega _m|)\theta (ϵ_l(ϵ_l+\omega _m))`$ (6)
$`\delta P_2`$ $`=`$ $`4\pi N(0)T{\displaystyle \underset{ϵ_l}{}}T{\displaystyle \underset{\omega _m}{}}{\displaystyle \frac{1}{aL}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{1}{(Dq^2+2|ϵ_l|)^2(D(q^{}+q)^2+|ϵ_l|+|ϵ_l+\omega _m|)}}V_C(q^{},|\omega _m|)\theta (ϵ_l(ϵ_l+\omega _m))`$ (7)
$`\delta P_3`$ $`=`$ $`4\pi N(0)T{\displaystyle \underset{ϵ_l}{}}T{\displaystyle \underset{\omega _m}{}}{\displaystyle \frac{1}{aL}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{Dq^2+Dq^2+|ϵ_l|+|ϵ_l+\omega _m|}{(Dq^2+2|ϵ_l|)(Dq^2+2|ϵ_l+\omega _m|)(Dq^2+|ϵ_l|+|ϵ_l+\omega _m|)^2}}V_C(q^{},|\omega _m|)\theta (ϵ_l(ϵ_l+\omega _m))`$ (8)
$`\delta P_4`$ $`=`$ $`4\pi N(0)T{\displaystyle \underset{ϵ_l}{}}T{\displaystyle \underset{\omega _m}{}}{\displaystyle \frac{1}{aL}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{1}{(Dq^2+2|ϵ_l|)(Dq^2+2|ϵ_l+\omega _m|)(D(q^{}+q)^2+|ϵ_l|+|ϵ_l+\omega _m|)}}V_C(q^{},|\omega _m|)\theta (ϵ_l(ϵ_l+\omega _m))`$ (9)
$`\delta P_5`$ $`=`$ $`4\pi ^2N(0)^2T{\displaystyle \underset{\omega _m}{}}{\displaystyle \frac{1}{aL}}{\displaystyle \underset{q^{}}{}}[T{\displaystyle \underset{ϵ_l}{}}\{{\displaystyle \frac{|\omega _m|}{(Dq^2+2|ϵ_l|)(Dq^2+2|ϵ_l+\omega _m|)(D(q+q^{})^2+|ϵ_l|+|ϵ_l+\omega _m|)}}\theta (ϵ_l(ϵ_l+\omega _m))`$ (11)
$`+{\displaystyle \frac{Dq^2+|\omega _m|}{(Dq^2+2|ϵ_l|)(Dq^2+2|ϵ_l+\omega _m|)(D(q+q^{})^2+|ϵ_l|+|ϵ_l+\omega _m|)}}\theta (ϵ_l(ϵ_l+\omega _m))\}]^2V_C(q^{},|\omega _m|)L(q+q^{},|\omega _m|)`$
where $`ϵ_l=(2l+1)\pi T`$ is a Fermi Matsubara frequency, $`\omega _m=2m\pi T`$ is a Bose Matsubara frequency, and $`a`$, $`L`$ are the width and length of the wire respectively. The diffusive form of the denominators in Eqn. (6) is valid only for momentum and frequency values such that $`Dq^2,|ϵ_l|,|\omega _m|<1/\tau `$, so we introduce the upper cutoff, $`1/\tau `$, in all these sums.

For the rest of this section we will work at zero magnetic field, so we set $`q=0`$. The terms $`\delta P_{15}`$ have previously been summed for 2D films to give the expression

$`\mathrm{ln}\left({\displaystyle \frac{T_c}{T_{c0}}}\right)`$ $`=`$ $`T{\displaystyle \underset{\omega _m}{}}{\displaystyle \frac{1}{aL}}{\displaystyle \underset{q^{}}{}}\{{\displaystyle \frac{1}{2\pi T}}{\displaystyle \frac{Dq^2}{\omega _m^2(Dq^2)^2}}\psi ^{}({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\omega _m}{2\pi T}})+{\displaystyle \frac{2Dq^2[\omega _m^2+(Dq^2)^2]}{\omega _m[\omega _m^2(Dq^2)^2]^2}}[\psi ({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\omega _m}{2\pi T}})\psi \left({\displaystyle \frac{1}{2}}\right)]`$ (12)
$``$ $`{\displaystyle \frac{4(Dq^2)^2}{[\omega _m^2(Dq^2)^2]^2}}{\displaystyle \frac{\left[\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\omega _m}{2\pi T}}\right)\psi \left({\displaystyle \frac{1}{2}}\right)\right]^2}{\left[\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{Dq^2+\omega _m}{4\pi T}}\right)\psi \left({\displaystyle \frac{1}{2}}\right)\right]}}\}V_C(q^{},\omega _m)`$ (13)
The $`1/q^2`$ singularity in the screened Coulomb potential $`V_C(q^{},\omega _m)`$ is cancelled since there is an overall factor of $`Dq^2`$ in the above equation. Note that there is, in fact, no divergence at $`Dq^2=\omega _m`$, but this can only be seen by Taylor expanding to second order. The need to perform this expansion makes the calculation involved, especially at the level of numerics. To allow simplification, we define $`y=Dq^2/2\pi T`$ and $`m=\omega _m/2\pi T`$, which leads to the result

$`\mathrm{ln}\left({\displaystyle \frac{T_c}{T_{c0}}}\right)`$ $`=`$ $`{\displaystyle \frac{1}{4\pi ^2N(0)T}}{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{1}{aL}}{\displaystyle \underset{q^{}}{}}[\{{\displaystyle \frac{2}{m}}[\psi ({\displaystyle \frac{1}{2}}+m)\psi \left({\displaystyle \frac{1}{2}}\right)]+\psi ^{}({\displaystyle \frac{1}{2}}+m)\}{\displaystyle \frac{1}{y+m}}`$ (14)
$``$ $`{\displaystyle \frac{4y}{(m+y)(my)^2}}\{{\displaystyle \frac{\left[\psi \left({\displaystyle \frac{1}{2}}+m\right)\psi \left({\displaystyle \frac{1}{2}}\right)\right]^2}{\left[\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{m+y}{2}}\right)\psi \left({\displaystyle \frac{1}{2}}\right)\right]}}[\psi ({\displaystyle \frac{1}{2}}+m)\psi \left({\displaystyle \frac{1}{2}}\right)]+{\displaystyle \frac{1}{2}}(ym)\psi ^{}({\displaystyle \frac{1}{2}}+m)\}]`$ (15)
where we have split off the most divergent term at large momentum and frequency, and $`M=1/2\pi T\tau `$. To proceed further we must consider the $`q^{}`$-sum in various dimensions. In two dimensions the sum becomes an integral, and it is most sensible to perform the cut-off in a circularly symmetric manner

$$\frac{1}{aL}\underset{q^{}}{}_0^{1/\tau }\frac{2\pi q^{}dq^{}}{(2\pi )^2}=\frac{T}{2D}_0^M𝑑y.$$
(16)
In quasi-1D wires, we must replace the integral in the transverse direction by a sum,

$$\frac{1}{aL}\underset{q^{}}{}\frac{1}{a}\underset{q_a^{}}{}_0^{1/\tau }\frac{dq_L^{}}{2\pi },$$
(17)
where $`q^2=q_L^2+q_a^2`$. At this point we have to decide which boundary conditions to use in the transverse direction. To illustrate the strong effect this has we will compare periodic and zero-gradient boundary conditions, the details of which are listed in Table 1. The zero-gradient boundary conditions are the physically correct ones for a superconductor-insulator boundary as discussed in Ref. \[\]. To perform the sum over $`q^{}`$ note that we can write $`y=y_L+y_a`$ where $`y_L=Dq_L^2/2\pi T`$ and $`y_a=Dq_a^2/2\pi T`$, and if we define $`z=\sqrt{y_L}`$, the $`q_L^{}`$-integral becomes

$$_0^{1/\tau }\frac{dq_L^{}}{2\pi }=_0^\sqrt{M}\frac{dz}{L_T}.$$
(18)
where $`L_T=\sqrt{2\pi D/T}`$ is the thermal length. Note that $`y=z^2+4p(a/L_T)^2`$ for zero-gradient boundary conditions and $`y=z^2+p(a/L_T)^2`$ for periodic boundary conditions. The final result of perturbation theory is then given by Eqn. (14) with the prefactor and sums becoming

$$\frac{1}{4\pi ^2N(0)T}\underset{m=1}{\overset{M}{}}\frac{1}{aL}\underset{q^{}}{}\{\begin{array}{cc}t\underset{m=1}{\overset{M}{}}_0^M𝑑y\hfill & \text{ }\text{2D}\hfill \\ \frac{t}{\pi }\left(\frac{L_T}{a}\right)\underset{m=1}{\overset{M}{}}\underset{p}{}_0^\sqrt{M}𝑑z\hfill & \text{ }\text{quasi–1D}\hfill \end{array}$$
(19)
In the above, $`t`$ is a dimensionless measure of disorder given by

$$t=\frac{1}{8\pi ^2N(0)D}=\frac{R_{\text{ }\text{ }\text{ }\text{ }\text{ }}}{R_0}$$
(20)
where $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$ is the normal state resistance per square, and $`R_0=2\pi h/e^2162k\mathrm{\Omega }`$. Note that $`t`$ only differs from the previously discussed parameter $`1/k_F\mathrm{}`$ by a factor $`2\pi `$.

Having shown that the low-momentum singularity in the screened Coulomb potential is cancelled by diagram 5, we recalculate the suppression of transition temperature using a featureless Coulomb potential of magnitude $`N(0)V_C=1/2`$, and including only diagrams 1–4. We have previously shown that the above equations simplify to give

$`\delta P_{12}`$ $`=`$ $`T{\displaystyle \underset{\omega _m>0}{}}\left[\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\omega _m}{2\pi T}}\right)\psi \left({\displaystyle \frac{1}{2}}\right)\right]{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{1}{(Dq^2+\omega _m)^2}}`$ (21)
$`\delta P_{34}`$ $`=`$ $`2T{\displaystyle \underset{\omega _m>0}{}}\left[\psi \left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{\omega _m}{2\pi T}}\right)\psi \left({\displaystyle \frac{1}{2}}\right)\right]{\displaystyle \frac{1}{\omega _m}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{1}{(Dq^2+\omega _m)}}.`$ (22)
For convenience we define dimensionless sums $`I_1(m)`$ and $`I_2(m)`$ as

$`I_1(m)`$ $``$ $`{\displaystyle \frac{4\pi D}{aL}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{1}{Dq^2+|\omega _m|}}`$ (23)
$`I_2(m)`$ $``$ $`{\displaystyle \frac{4\pi D}{aL}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{2\pi T}{(Dq^2+|\omega _m|)^2}}`$ (24)
which are evaluated in two dimensions and quasi-one dimension in Appendix A. When we combine Eqns. (5) and (21) we obtain the formula for transition temperature suppression assuming a featureless Coulomb potential,

$$\mathrm{ln}\left(\frac{T_c}{T_{c0}}\right)=t\underset{m=1}{\overset{M}{}}\left[\psi \left(\frac{1}{2}+m\right)\psi \left(\frac{1}{2}\right)\right]\left[I_2(m)+\frac{2I_1(m)}{m}\right].$$
(25)
If we ignore the self-energy diagrams ($`\delta P_{12}`$), the $`I_2(m)`$ term in the above equation vanishes to leave

$$\mathrm{ln}\left(\frac{T_c}{T_{c0}}\right)=2t\underset{m=1}{\overset{M}{}}\left[\psi \left(\frac{1}{2}+m\right)\psi \left(\frac{1}{2}\right)\right]\frac{2I_1(m)}{m}.$$
(26)
Having derived expressions for the suppression of transition temperature using various approximations and boundary conditions, we perform a detailed numerical analysis, the results of which can be summarised as follows. The transition temperature is very sensitive to the boundary conditions used – periodic boundary conditions do not yield a dependence of transition temperature on wire width, whilst zero-gradient boundary conditions show a systematic suppression as wire width is decreased. The predictions are also sensitive to the value of the thermal length, $`L_T`$, since the width, $`a`$, enters only in the dimensionless combination $`a/L_T`$. Experimentally $`L_T`$ may be known only roughly and it may be necessary to fit its value. On the other hand, the predictions are not very sensitive to either the form of the Coulomb interaction used (Eqns. (14) and (25)), or the inclusion of self-energy corrections (Eqns. (25) and (26)). The slight difference between these approaches vanishes when we fit the initial slope of the $`T_c(R_{\text{ }\text{ }\text{ }\text{ }\text{ }})`$ curve to experiment, which sets the value of the BCS attraction, $`\gamma `$. This demonstrates that the simplest approach – featureless Coulomb potential and no self-energy diagrams – is valid. This is important since these simplifications are necessary for the resummation technique discussed in the next section. It is necessary to use this technique to cure the unphysical re-entrance problems which beset perturbation theory at low temperature.

## III The Resummation Technique

In this section we will review the resummation technique developed by Oreg and Finkel’stein (OF) as applied to the problem of $`T_c`$ suppression in quasi-1D wires. This method involves calculating the pair-amplitude, $`\mathrm{\Gamma }(ϵ_n,ϵ_l)`$, which is a matrix with Fermi Matsubara frequencies labelling its rows and columns. The transition temperature is identified as the temperature at which this matrix becomes singular. OF obtained an equation for $`\mathrm{\Gamma }`$ by including only pseudopotential diagrams, and solved it using matrix perturbation theory with zero-gradient boundary conditions. We have performed calculations to show that it is legitimate to include only these diagrams, and that it is essential to use zero-gradient boundary conditions, as described in section II for perturbation theory. In a previous paper we have shown that the equation for $`\mathrm{\Gamma }`$ can be solved exactly at no extra numerical cost, so the matrix perturbation theory is unnecessary, and indeed gives incorrect results at large disorder. Therefore we shall include only pseudopotential diagrams and solve the equation for $`\mathrm{\Gamma }`$ using an exact diagonalization method with zero-gradient boundary conditions.

Before we proceed further, let us briefly describe the details of the technique. The ladder summation involved is shown in Fig. 2. The block $`t\mathrm{\Lambda }`$, an effective Coulomb pseudopotential, contains diagrams 3 and 4 of Fig. 1(c); the block $`\mathrm{\Sigma }`$, the self-energy of the Cooperon impurity ladder, $`C`$, contains diagrams 1 and 2 of Fig. 1(c). We have included the self-energy diagrams for completeness, although their effects can be ignored in calculations, as discussed above. The equation for $`\mathrm{\Gamma }`$ is seen to be

$$\mathrm{\Gamma }(ϵ_n,ϵ_l)=|\gamma |+t\mathrm{\Lambda }(ϵ_n,ϵ_l)2\pi T\underset{m=0}{\overset{M}{}}\left[|\gamma |+t\mathrm{\Lambda }(ϵ_n,ϵ_m)\right]\frac{1}{|ϵ_m|}\mathrm{\Gamma }(ϵ_m,ϵ_l)$$
(27)
where the contributions to $`t\mathrm{\Lambda }`$ are

$`t\mathrm{\Lambda }_3(ϵ_n,ϵ_m)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi N(0)}}{\displaystyle \frac{1}{aL}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{Dq^2+Dq^2+|ϵ_n|+|ϵ_m|}{[Dq^2+|ϵ_n|+|ϵ_m|]^2}}\theta (ϵ_nϵ_m)`$ (28)
$`t\mathrm{\Lambda }_4(ϵ_n,ϵ_m)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi N(0)}}{\displaystyle \frac{1}{aL}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{1}{[Dq^2+|ϵ_n|+|ϵ_m|]}}\theta (ϵ_nϵ_m).`$ (29)
We have kept non-zero external momentum $`q`$ so that we can include magnetic field effects later. For now we consider the case of zero magnetic field, so we set $`q=0`$, and the two contributions to $`t\mathrm{\Lambda }`$ are identical. It follows that $`\mathrm{\Lambda }(ϵ_n,ϵ_m)`$ does not depend upon the relative sign of $`ϵ_n`$ and $`ϵ_m`$ and has the form

$$\mathrm{\Lambda }(ϵ_n,ϵ_m)=I_1(n+m+1),$$
(30)
where $`I_1`$ is defined in Eqn. (23). Eqn. (27) can then be written in matrix form as

$$\widehat{\mathrm{\Gamma }}=|\gamma |\widehat{1}+t\widehat{\mathrm{\Lambda }}[|\gamma |\widehat{1}+t\widehat{\mathrm{\Lambda }}]\widehat{ϵ}^1\widehat{\mathrm{\Gamma }},$$
(31)
where $`\widehat{\mathrm{\Gamma }}_{nm}=\mathrm{\Gamma }(ϵ_n,ϵ_m)`$, $`\widehat{1}_{nm}=1`$, $`\widehat{\mathrm{\Lambda }}_{nm}=\mathrm{\Lambda }(ϵ_n,ϵ_m)`$ and $`\widehat{ϵ}_{nm}=(n+1/2)\delta _{nm}`$. This has the solution

$$\widehat{\mathrm{\Gamma }}=\widehat{ϵ}^{1/2}(\widehat{I}|\gamma |\widehat{\mathrm{\Pi }})^1\widehat{ϵ}^{1/2}(|\gamma |\widehat{1}+t\widehat{\mathrm{\Lambda }}),$$
(32)
where

$$\widehat{\mathrm{\Pi }}=\widehat{ϵ}^{1/2}[\widehat{1}|\gamma |^1t\widehat{\mathrm{\Lambda }}]\widehat{ϵ}^{1/2},$$
(33)
and $`\widehat{I}_{nm}=\delta _{nm}`$ is the identity matrix. It follows that when the matrix $`\widehat{\mathrm{\Pi }}`$ has an eigenvalue equal to $`1/|\gamma |`$, the pair amplitude diverges, and we have found $`T_c`$. In other words, we need to solve the eigenvalue equation

$$\left[|\gamma |^1\widehat{I}\widehat{\mathrm{\Pi }}(T_c)\right]|\psi =0.$$
(34)
We solve Eqn. (34) by diagonalizing the matrix $`[|\gamma |^1\widehat{I}\widehat{\mathrm{\Pi }}(T)]`$ numerically, and decreasing the temperature until its lowest eigenvalue equals zero. Since we only need the lowest eigenvalue, and the matrix has such a simple structure, the Lanczos method is very efficient here. The only minor complication is that the matrix $`\widehat{\mathrm{\Pi }}(T)`$ depends upon $`T`$ both through the dependence of $`\widehat{\mathrm{\Lambda }}`$ upon $`T`$ and through its rank $`M=1/2\pi T\tau `$. We therefore start at the mean-field value of $`M`$, which we will call $`M_0`$, and decrease the temperature by increasing $`M`$ successively by one. We diagonalize the matrix $`[|\gamma |^1\widehat{I}\widehat{\mathrm{\Pi }}]`$ for each value of $`M`$ until the lowest eigenvalue changes sign. At this point we have found $`T_c`$ for the given problem, and $`T_c/T_{c0}=M_0/M`$. This is obviously the method of choice, since no approximation is involved.

We can now compare theory and experiment: there are three sets of experimental data that we know about – those of Graybeal et al (1987), Sharifi et al (1993), and Xiong et al (1997). The data of Graybeal et al consists of $`T_c`$ vs $`a`$ for wires of constant thickness (and hence constant resistance per square). Their data shows that the relative suppression of $`T_c`$ varies as $`1/a^2`$ with a proportionality constant independent of $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$. We have performed calculations which predict a dependence roughly of the form $`1/a`$, but with a proportionality constant dependent on $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$. We will therefore not consider this data further. The data of Sharifi et al and Xiong et al are both of the form of $`T_c(R_{\text{ }\text{ }\text{ }\text{ }\text{ }},a)`$ vs $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$ for wires of different width, but only the latter shows a consistent trend with decreasing width, so it is only this data that we shall attempt to fit. The fact that this is the most recent data is reassuring, since it suggests that more recent fabrication techniques have produced better quality wires.

In Fig. 3 we compare the data of Xiong et al with our theoretical predictions. In doing this we encounter two difficulties. The first is that experiment and theory do not agree for some of the 2D films used. We assume that this is due to some additional physical processes that will affect samples in the same manner independent of width. We therefore multiply our theoretical predictions by $`T_c(\text{2D,expt})/T_c(\text{2D,theory})`$ to factor out this effect. The second difficulty is that wires of different widths do not appear to have the same value of $`T_c`$ in the clean limit, $`T_c(0,a)`$. We therefore assume a different $`T_c(0,a)`$ for each wire. Upon making these two assumptions, we see that the agreement between theory and experiment is good. However plotting the data in this form emphasises the suppression of transition temperature due to disorder, whereas we are interested in the additional suppression due to the decrease in width. For this reason we need to replot the data to emphasise this latter point.

In Fig. 4 we replot Fig. 3 in terms of the additional relative suppression, $`\mathrm{\Theta }(R_{\text{ }\text{ }\text{ }\text{ }\text{ }},a)`$, which is defined by

$$\mathrm{\Theta }(R_{\text{ }\text{ }\text{ }\text{ }\text{ }},a)=\frac{T_c(R_{\text{ }\text{ }\text{ }\text{ }\text{ }},a)/T_c(0,a)}{T_c(R_{\text{ }\text{ }\text{ }\text{ }\text{ }},2D)/T_c(0,2D)}.$$
(35)
This quantity focuses on the basic physical phenomenon we are interested in – the additional suppression of transition temperature caused by the reduced dimensionality in the quasi-1D wires. It has the additional benefit of factoring out the two difficulties discussed in the previous paragraph. When plotted in this form, the agreement between theory and experiment no longer appears to be as good. This is because our theory, and indeed any theory, predicts a smooth dimensional crossover whereas the experimental data shows an abrupt crossover between the 2D film and the $`1000\AA `$ wire, which lie on one curve, and the $`580\AA ,350\AA ,250\AA `$ wires, which lie on a different curve. The fact that the data for the three narrowest wires fall on top of each other might be due to them having the same effective dimensionless width $`a/L_T`$. If we allowed ourselves the additional freedom of fitting this parameter separately for each wire, we would obtain good agreement between theory and experiment, but we do not believe this is justified. We have used the value of the thermal length quoted in Sharifi et al (viz $`L_T=400\AA `$ at $`T=4K`$, which corresponds to $`L_T=300\AA `$ at $`T_{c0}=6.9K`$).

The basic problem is clearly that there is not enough experimental data available – there are only 19 points in total, and each of these corresponds to a different wire. To carry out a more systematic comparison between theory and experiment, one would need wires of many different resistances (thicknesses) for each value of width. This would, however, be very costly in terms of experimental effort. A less costly method for generating more data would be to apply a magnetic field to a given wire. We therefore extend our theory to include the effect of such a field in the next section.

## IV Effect of an Applied Magnetic Field

We will now consider the effect of an applied perpendicular magnetic field. There are two effects to consider here: the orbital effect, caused by the magnetic field scrambling the relative phase of the electrons in a Cooper pair, and the Zeeman effect, caused by the magnetic field lifting the degeneracy of up and down spins. Let us first consider the relative importance of the two effects. The orbital effect will destroy superconductivity when the magnetic dephasing rate, $`1/\tau _H`$, is of the order of $`k_BT_c`$, which gives $`H4T`$ for a Pb wire of width $`250\AA `$. The Zeeman splitting will destroy superconductivity when its magnitude $`\mu _BH`$ is of the order of $`k_BT_c`$, which corresponds to $`H10T`$ for Pb. Since these fields have the same order of magnitude, we will need to consider both effects.

We first examine the orbital effect. The magnetic field causes the eigenvalues of the Cooperon operator to change from the zero-field values

$$Dq^2=D\left[\left(\frac{2\pi n}{L}\right)^2+\left(\frac{\pi p}{a}\right)^2\right],n=0,\pm 1,\pm 2\mathrm{},p=0,1,2\mathrm{}$$
(36)
to the eigenvalues, $`\lambda `$, of the minimally coupled Cooperon operator

$$D(i+2e𝐀)^2\psi (x,y)=D\left[\left(i\frac{}{x}+2eHy\right)^2+\left(i\frac{}{y}\right)^2\right]\psi (x,y)=\lambda \psi (x,y),$$
(37)
where we have chosen the gauge $`𝐀=(Hy,0,0)`$. The eigenfunctions $`\psi (x,y)`$ must satisfy zero-gradient boundary conditions on the sides of the wire, $`y=\pm a/2`$, and periodic boundary conditions at the ends of the wire, $`x=0,L`$. The solution has the separable form

$$\psi (x,y)=\varphi _{q_n,p}(y)e^{iqx},q_n=\frac{2\pi n}{L},n=0,\pm 1,\pm 2\mathrm{}$$
(38)
where $`\varphi _{q_n,p}`$ satisfies the equation

$$D\left[\frac{d^2}{dy^2}+(q+2eHy)^2\right]\varphi _{q_n,p}(y)=\lambda _{q_n,p}\varphi _{q_n,p}(y),\frac{d\varphi _{q_n,p}}{dy}(\pm a/2)=0,$$
(39)
and $`p=0,1,2\mathrm{}`$ labels the eigenvalues in ascending order.

To solve this eigenvalue problem numerically, we convert it into a matrix diagonalization problem by working in a truncated function space spanned by eigenfunctions of the zero-field problem. The lowest eigenvalue, $`\lambda _{0,0}`$, equals $`2/\tau _H`$, where $`1/\tau _H`$ is the magnetic dephasing rate which appears in the Abrikosov-Gor’kov equation for $`T_c(H)`$. There are two asymptotic limits where $`1/\tau _H`$ is easily determined: the 2D limit, where the problem becomes that of a simple harmonic oscillator, with $`1/\tau _H=DeH`$; and the 1D limit, where the eigenfunction is simply a constant, and hence (using perturbation theory) $`1/\tau _H=De^2H^2a^2/6`$. The 2D to 1D crossover occurs when the width, $`a`$, is of the order of the magnetic length, $`L_H=(\mathrm{}/eH)^{1/2}`$, which occurs at the field $`H_0=h/ea^2`$. At larger magnetic fields the cyclotron orbit of the electron lies entirely within the wire, and hence we are in the 2D limit; at smaller magnetic fields the electron hits the sides of the wire before completing a cyclotron orbit, and hence we are in the 1D limit. To summarise,

$$\frac{1}{\tau _H}=\{\begin{array}{cc}DeH\hfill & HH_0\text{ (}aL_H\text{)\hspace{1em}2D}\hfill \\ De^2H^2a^2/6\hfill & HH_0\text{ (}aL_H\text{)\hspace{1em}1D}\hfill \end{array}$$
(40)
These eigenvalues then enter into our equations for the pair propagator in two ways: (i) the term $`Dq^2`$ in the $`q`$-dependent pair propagator is replaced by the lowest eigenvalue $`\lambda _{0,0}`$; (ii) sums over functions of $`Dq^2`$ in diagrams which involve Cooperonic ladders (e.g. diagrams 2 and 4 in Fig. 1(c)) are replaced by sums over $`\lambda _{q_n,p}`$. Obviously sums over functions of $`Dq^2`$ in diagrams which involve diffusion ladders (e.g. diagrams 1 and 3 in Fig. 1(c)) remain unchanged. It is numerically very expensive to perform the sums over all the eigenvalues, so we would like to see how important it is to use these eigenvalues in any sums. To answer this question we again use perturbation theory as our testing ground. We perform the Cooperonic sum in two different ways – both correctly, and by replacing it by the corresponding diffusion sum. In both cases we still make the replacement of $`Dq^2`$ by $`\lambda _{0,0}=2/\tau _H`$. The perturbative equations are then

$`{\displaystyle \frac{\delta P_1}{N(0)}}`$ $`=`$ $`{\displaystyle \frac{t}{2}}{\displaystyle \underset{m=1}{\overset{M}{}}}\left\{2\left[\psi \left({\displaystyle \frac{1}{2}}+m+\alpha \right)\psi \left({\displaystyle \frac{1}{2}}+\alpha \right)\right]I_2(m)+\left[\psi ^{}\left({\displaystyle \frac{1}{2}}+\alpha \right)\psi ^{}\left({\displaystyle \frac{1}{2}}+m+\alpha \right)\right]I_1(m)\right\}`$ (41)
$`{\displaystyle \frac{\delta P_2}{N(0)}}`$ $`=`$ $`+{\displaystyle \frac{t}{2}}{\displaystyle \underset{m=1}{\overset{M}{}}}\left[\psi ^{}\left({\displaystyle \frac{1}{2}}+\alpha \right)\psi ^{}\left({\displaystyle \frac{1}{2}}+m+\alpha \right)\right]C_1(m)`$ (42)
$`{\displaystyle \frac{\delta P_3}{N(0)}}`$ $`=`$ $`t{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{1}{m+2\alpha }}\left[\psi \left({\displaystyle \frac{1}{2}}+m+\alpha \right)\psi \left({\displaystyle \frac{1}{2}}+\alpha \right)\right](I_1(m)+2\alpha I_2(m))`$ (43)
$`{\displaystyle \frac{\delta P_4}{N(0)}}`$ $`=`$ $`t{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{1}{m+2\alpha }}\left[\psi \left({\displaystyle \frac{1}{2}}+m+\alpha \right)\psi \left({\displaystyle \frac{1}{2}}+\alpha \right)\right]C_1(m)`$ (44)
where $`\alpha =1/2\pi T\tau _H`$. The diffusive sums $`I_1(m)`$ and $`I_2(m)`$ are defined as in Eqn. (23) except that here we perform the momentum sums in the longitudinal direction exactly to give

$`I_1(m)`$ $`=`$ $`{\displaystyle \frac{4\pi D}{aL}}{\displaystyle \underset{q}{}}{\displaystyle \frac{1}{Dq^2+2\pi Tm}}={\displaystyle \frac{L_T^2}{\pi La}}{\displaystyle \underset{n=n_0}{\overset{n_0}{}}}{\displaystyle \underset{k=0}{\overset{k_0}{}}}{\displaystyle \frac{1}{\left[m+\left(\frac{kL_T}{2a}\right)^2+\left(\frac{nL_T}{L}\right)^2\right]}}`$ (45)
$`I_2(m)`$ $`=`$ $`{\displaystyle \frac{4\pi D}{aL}}{\displaystyle \underset{q}{}}{\displaystyle \frac{2\pi T}{(Dq^2+2\pi Tm)^2}}={\displaystyle \frac{L_T^2}{\pi La}}{\displaystyle \underset{n=n_0}{\overset{n_0}{}}}{\displaystyle \underset{k=0}{\overset{k_0}{}}}{\displaystyle \frac{1}{\left[m+\left(\frac{kL_T}{2a}\right)^2+\left(\frac{nL_T}{L}\right)^2\right]^2}}`$ (46)
and the Cooperonic sum, $`C_1`$, which is now distinct from $`I_1`$, is given by

$$C_1(m)=\frac{4\pi D}{aL}\underset{q}{}\frac{1}{Dq^2+2\pi Tm}=\frac{L_T^2}{\pi La}\underset{n=n_0}{\overset{n_0}{}}\underset{k=0}{\overset{k_0}{}}\frac{1}{\left[m+\lambda _{nk}/2\pi T\right]}.$$
(47)
We found essentially no difference between using the correct $`C_1`$, or replacing it by $`I_1`$, over a wide range of magnetic fields. It follows that the only replacement we need to make is of $`Dq^2`$ by $`\lambda _{0,0}`$, which greatly reduces the numerical cost.

Let us next consider the effect of the Zeeman splitting. This leads to changes in both the mean field and perturbative terms, and we will show that only the former effect is important. To do this we will evaluate the mean-field pair polarization function and the first-order perturbation correction, $`\delta P`$, both with and without Zeeman splitting. We find that the latter is essentially unaffected. To treat the Zeeman effect, we need formulas for the Cooperon impurity ladder and Hikami boxes, where we note that the boxes in diagrams 1 and 3 of Fig. 1(c) will give different contributions. The non-zero elements of the spin-dependent Cooperon ladder, $`C_{\alpha \beta \gamma \delta }`$, are

$$C_{}=\frac{1}{2\pi N(0)\tau ^2}\frac{1}{Dq^2+|\omega |+|\omega ^{}|2i\mathrm{\Delta }},C_{}=\frac{1}{2\pi N(0)\tau ^2}\frac{1}{Dq^2+|\omega |+|\omega ^{}|+2i\mathrm{\Delta }}=C_{}^{}$$
(48)
where $`\mathrm{\Delta }=\mu _BH`$ is the Zeeman splitting. The appropriate spin combination which occurs in an s-wave pair propagator is $`C_{}+C_{}=2\mathrm{}eC_{}`$. This leads to the mean field formula for $`T_c`$,

$$\mathrm{ln}\left(\frac{T_c}{T_{c0}}\right)=\psi \left(\frac{1}{2}\right)\mathrm{}e\left[\psi \left(\frac{1}{2}+\kappa \right)\right]$$
(49)
where $`\kappa =\alpha i\beta `$ with $`\beta =\mathrm{\Delta }/2\pi T_c=\mu _BH/2\pi T_c`$. To calculate the first-order pair polarization terms $`\delta P_{14}`$, we note that: (a) the diffusion ladders on Coulomb vertices are unchanged; (b) Cooperon ladders pick up an extra $`2i\mathrm{\Delta }`$ as described above; (c) the Hikami box in $`\delta P_1`$ picks up an extra $`2i\mathrm{\Delta }`$, whilst the box in $`\delta P_3`$ is unaffected. This leads to the formulas

$`{\displaystyle \frac{\delta P_1}{N(0)}}`$ $`=`$ $`{\displaystyle \frac{t}{2}}\mathrm{}e{\displaystyle \underset{m=1}{\overset{M}{}}}\left\{2\left[\psi \left({\displaystyle \frac{1}{2}}+m+\kappa \right)\psi \left({\displaystyle \frac{1}{2}}+\kappa \right)\right]I_2(m)+\left[\psi ^{}\left({\displaystyle \frac{1}{2}}+\kappa \right)\psi ^{}\left({\displaystyle \frac{1}{2}}+m+\kappa \right)\right]I_1(m)\right\}`$ (50)
$`{\displaystyle \frac{\delta P_2}{N(0)}}`$ $`=`$ $`+{\displaystyle \frac{t}{2}}\mathrm{}e{\displaystyle \underset{m=1}{\overset{M}{}}}\left[\psi ^{}\left({\displaystyle \frac{1}{2}}+\kappa \right)\psi ^{}\left({\displaystyle \frac{1}{2}}+m+\kappa \right)\right]\overline{C}_1(m)`$ (51)
$`{\displaystyle \frac{\delta P_3}{N(0)}}`$ $`=`$ $`t\mathrm{}e{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{1}{m+2\kappa }}\left[\psi \left({\displaystyle \frac{1}{2}}+m+\kappa \right)\psi \left({\displaystyle \frac{1}{2}}+\kappa \right)\right](I_1(m)+2\alpha I_2(m))`$ (52)
$`{\displaystyle \frac{\delta P_4}{N(0)}}`$ $`=`$ $`t\mathrm{}e{\displaystyle \underset{m=1}{\overset{M}{}}}{\displaystyle \frac{1}{m+2\kappa }}\left[\psi \left({\displaystyle \frac{1}{2}}+m+\kappa \right)\psi \left({\displaystyle \frac{1}{2}}+\kappa \right)\right]\overline{C}_1(m)`$ (53)
where $`I_1(m)`$ and $`I_2(m)`$ are the same as before, whilst

$$\overline{C}_1(m)=\frac{4\pi D}{aL}\underset{q}{}\frac{1}{Dq^2+2\pi Tm2i\mathrm{\Delta }}=\frac{L_T^2}{\pi aL}\underset{n=n_0}{\overset{n_0}{}}\underset{k=0}{\overset{k_0}{}}\frac{1}{\left[m+2\kappa +\left(\frac{kL_T}{2a}\right)^2+\left(\frac{nL_T}{L}\right)^2\right]}$$
(54)
Calculating the $`\delta P_{14}`$ as above, or with $`\beta `$ set to zero, we find essentially no difference over a wide range of magnetic fields until we get to widths of order $`10\AA `$. It follows that we only need consider the Zeeman splitting at the mean field level. Even this is a small effect, only noticeable for wires narrower than about $`150\AA `$.

Let us finally perform calculations of $`T_c`$ using the resummation technique. From the discussions above, we now know that to include magnetic field effects we need only replace $`Dq^2`$ in Cooperons by $`2/\tau _H`$, which must be calculated numerically, and include the Zeeman splitting at the mean field level. The Cooperon denominator $`Dq^2+2|ϵ_m|`$ then becomes $`2|ϵ_m|+2/\tau _H2i\mathrm{\Delta }`$, and hence $`\widehat{ϵ}_{nm}=(n+1/2+\kappa )\delta _{nm}`$, which is now complex. Because the two spin configurations which occur in an s-wave pair propagator are complex conjugates of each other, we end up taking the real part of the complex matrix equation. The equation for $`\mathrm{\Lambda }_{nm}`$ becomes

$$\mathrm{\Lambda }_{nm}=\{\begin{array}{cc}I_1(n+m+1+2\alpha )\hfill & \text{ }ϵ_nϵ_m>0\hfill \\ I_1(n+m+1+2\alpha )+2I_2(n+m+1+2\alpha )\hfill & \text{ }ϵ_nϵ_m<0\hfill \end{array}$$
(55)
where we have not included the Zeeman effect as we have shown this to be much smaller in perturbation theory.

We are now in a position to calculate the transition temperature, $`T_c`$, as a function of resistance per square, $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$, wire width, $`a`$, and applied magnetic field, $`H`$. This obviously leads to a wide choice of ways to present results. Experimentally one would probably vary magnetic field on a sample of fixed width and thickness, and measure $`T_c(H)`$. In Fig. 5 we therefore plot $`T_c(H)`$ for wires of different resistances per square, $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}=0\mathrm{\Omega },1000\mathrm{\Omega },2000\mathrm{\Omega }`$, and widths $`a=500\AA ,600\AA ,750\AA ,2D`$ film. There is far too much data on this plot, but some general trends are apparent. Three physical processes are occuring here: increasing $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$ increases the effective Coulomb repulsion which strongly suppresses $`T_c`$; decreasing $`a`$ means that a larger magnetic field is required to fit a flux quantum inside the wire, and hence $`H_{c2}`$ increases at fixed $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$; decreasing $`a`$ also causes a crossover from two to one-dimensional behavior, and so increases fluctuation effects, leading to a small reduction of $`T_c`$ at fixed $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$. The last two effects mean that the $`T_c(H)`$ curves for wires with the same $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$ and different $`a`$ cross each other at non-zero field.

Another effect of the dimensional crossover is to change the shape of the $`T_c(H)`$ curve. This occurs because the magnetic pair-breaking rate, $`1/\tau _H`$, is linear in magnetic field in the two-dimensional limit, but quadratic in magnetic field in the one-dimensional limit. This leads to two limiting forms for the dimensionless $`T_c(H)`$ curve,

$$\mathrm{ln}\left(\frac{T_c(H)}{T_c(0)}\right)=\psi \left(\frac{1}{2}\right)\psi \left(\frac{1}{2}+e^{\psi (1/2)}\left[\frac{H}{H_{c2}}\right]^p\left[\frac{T_c(0)}{T_c(H)}\right]\right),$$
(56)
where $`p=1`$ in the 2D limit and $`p=2`$ in the 1D limit. We can move from the 2D limit to the 1D limit either by decreasing $`a`$ at fixed $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$, as shown in Fig. 6, or by increasing $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$ at fixed $`a`$, as shown in Fig. 7. Considering first the width dependence, we see in Fig. 6 that the dimensional crossover occurs when $`aL_T`$. For the wider wires the $`T_c(H)`$ curve joins the two-dimensional curve at the value of the magnetic field such that the magnetic length, $`L_H`$, is of order the width, $`a`$; for narrower wires this value of $`H`$ exceeds the two-dimensional value of $`H_{c2}`$, and the curves never meet. For ultranarrow wires, $`a<250\AA `$, the Zeeman effect becomes important. This causes the $`T_c(H)`$ curves to change shape again, moving slowly back in the direction of the 2D limit. Now let us consider the resistance dependence of the $`T_c(H)`$ curve. For narrow wires, $`a250\AA `$, all $`T_c(H)`$ curves fall on the 1D curve independent of $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$. For wide wires, $`a10000\AA `$, all $`T_c(H)`$ curves fall on the 2D curve independent of $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$. For wires of intermediate width, the $`T_c(H)`$ curves can be changed from the 1D curve to the 2D curve by increasing $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$, as shown in Fig. 7. Finally for ultranarrow wires, $`a50\AA `$, the $`T_c(H)`$ curves can again be tuned by varying $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$, in this case the tuning being between one-dimensional orbit and Zeeman curves.

## V Discussion and Conclusions

In this paper we have considered the effect of localization and interaction on transition temperature in ultranarrow superconducting wires. Our main conclusion is that, taking into account the relatively small amount of experimental data, theory and experiment are consistent with each other. The agreement is not perfect: the experimental data seems to show a very sudden crossover from 2D to 1D behavior when the wire width becomes of the order of the thermal length, whilst theory predicts a smooth crossover. Indeed, when plotted as the additional relative suppression of transition temperature due to finite width, the experimental data falls on two curves – one for the 2D film and $`1000\AA `$ wire; the other for the $`580\AA `$, $`350\AA `$ and $`250\AA `$ wires. These conclusions are based on a single experiment, as this is the only experiment which shows any consistent trend with wire width. More experimental data is obviously needed to investigate this further. Given the difficulty in fabricating these wires, we suggest that the magnetic field dependence of transition temperature be measured in future experiments. We have therefore made predictions for magnetic field dependence of transition temperature in wires of varying width and resistance in the hope that such measurements will be made in the future. Our main prediction is of a change in the shape of the reduced $`T_c(H)`$ curves due to the dimensional crossover. This may be tuned either by varying wire width or resistance.

Our work is based on the resummation technique developed by Oreg and Finkel’stein, which is a very powerful method of going beyond perturbation theory, and hence giving believable results at strong disorder. The suppression of the transition temperature is very sensitive to boundary conditions, due to the different relative weights of the zero momentum mode. We have used the physically correct zero-gradient boundary conditions in our work. The unexpected strong dependence of results on boundary conditions led us to worry about all details of the calculation. We therefore also investigated the effects of using a featureless Coulomb interaction instead of the exact one, ignoring self-energy contributions, and solving the OF matrix equation exactly instead of using an approximate method. We found that it is legitimate to use a featureless Coulomb interaction and ignore self-energy contributions, but that our matrix equations should be solved by exact diagonalization (which in any case is no more difficult).

In our calculations of the magnetic field dependence of transition temperature, we find that the orbital effect is the most important, and can be treated by replacing the $`Dq^2`$ terms in the momentum dependent pair propagator, $`L(q,0)`$, by twice the magnetic dephasing rate, $`2/\tau _H`$. The latter is the lowest eigenvalue of the minimally coupled Cooperon operator, and must be found numerically. The Zeeman splitting is found to have a non-negligible effect only at the mean field level, and even this effect is small. The perturbative term only becomes important in very narrow wires (of order $`10\AA `$), which cannot currently be fabricated.

ACKNOWLEDGEMENTS

We thank I. Aleiner, A. Clerk, A.M. Finkel’stein, I.V. Lerner, Y. Oreg and M. Steiner for helpful discussions, and F. Sharifi and P. Xiong for making their data available to us. RAS acknowledges support from the Nuffield Foundation and the UK EPSRC under Grant No GR/M98975. BSH acknowledges support from the UK EPSRC. VA is supported by the US National Science Foundation under Grant No DMR-9805613.

## A Details of Some Momentum Sums

In this appendix we derive expressions for the quantitities $`I_1(m)`$ and $`I_2(m)`$ defined in Eqn. (23) as

$`I_1`$ $``$ $`{\displaystyle \frac{4\pi D}{aL}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{1}{Dq^2+|\omega _m|}}`$ (A1)
$`I_2`$ $``$ $`{\displaystyle \frac{4\pi D}{aL}}{\displaystyle \underset{q^{}}{}}{\displaystyle \frac{2\pi T}{(Dq^2+|\omega _m|)^2}}`$ (A2)
In two dimensions the sum becomes an integral, and cutting off in a circularly symmetric manner,

$$\frac{1}{aL}\underset{q^{}}{}=_0^{Dq^2=1/\tau }\frac{2\pi q^{}dq^{}}{(2\pi )^2}=_0^{1/\tau }\frac{d(Dq^2)}{4\pi D},$$
(A3)
we obtain the results

$$I_1(m)=\mathrm{ln}\left(\frac{M}{m}\right),I_2(m)=\frac{1}{m}.$$
(A4)
In quasi-1D wires we replace the integral in the transverse direction by a sum, so that $`I_1(m)`$ becomes

$`I_1(m)`$ $`=`$ $`4\pi D{\displaystyle \frac{1}{a}}{\displaystyle \underset{q_a^{}}{}}{\displaystyle \frac{1}{2\pi }}{\displaystyle _0^{Dq_L^2=1/\tau }}{\displaystyle \frac{dq_L^{}}{Dq_L^2+Dq_a^2+2\pi Tm}}`$ (A5)
$`=`$ $`{\displaystyle \frac{4}{a}}{\displaystyle \underset{q_a^{}}{}}{\displaystyle \frac{1}{\sqrt{q_a^2+{\displaystyle \frac{2\pi Tm}{D}}}}}\mathrm{tan}^1\left[{\displaystyle \frac{1/\sqrt{D\tau }}{\sqrt{q_a^2+{\displaystyle \frac{2\pi Tm}{D}}}}}\right],`$ (A6)
and $`I_2(m)`$ can be found by differentiating the above with respect to $`m`$ and dividing by $`2\pi T`$. Substituting in the two types of boundary condition listed in Table 1 yields

$`I_1(m)=\{\begin{array}{cc}{\displaystyle \frac{2}{\pi }}{\displaystyle \underset{p=p_0}{\overset{p_0}{}}}(p^2+mA^2)^{\frac{1}{2}}\mathrm{tan}^1\left[{\displaystyle \frac{p_0}{\sqrt{p^2+mA^2}}}\right]\hfill & \text{ }\text{PBC}\hfill \\ {\displaystyle \frac{4}{\pi }}{\displaystyle \underset{p=0}{\overset{2p_0}{}}}(p^2+4mA^2)^{\frac{1}{2}}\mathrm{tan}^1\left[{\displaystyle \frac{2p_0}{\sqrt{p^2+4mA^2}}}\right]\hfill & \text{ }\text{ZGBC}\hfill \end{array}`$ (A7)
and

$`I_2(m)=\{\begin{array}{cc}{\displaystyle \frac{A^2}{\pi }}{\displaystyle \underset{p=p_0}{\overset{p_0}{}}}(p^2+mA^2)^{\frac{3}{2}}\left(\mathrm{tan}^1\left[{\displaystyle \frac{p_0}{\sqrt{p^2+mA^2}}}\right]+{\displaystyle \frac{p_0\sqrt{p^2+mA^2}}{p^2+(m+M)A^2}}\right)\hfill & \text{ }\mathrm{PBC}\hfill \\ {\displaystyle \frac{8A^2}{\pi }}{\displaystyle \underset{p=0}{\overset{2p_0}{}}}(p^2+4mA^2)^{\frac{3}{2}}\left(\mathrm{tan}^1\left[{\displaystyle \frac{2p_0}{\sqrt{p^2+4mA^2}}}\right]+{\displaystyle \frac{2p_0\sqrt{p^2+4mA^2}}{p^2+4(m+M)A^2}}\right)\hfill & \text{ }\mathrm{ZGBC}\hfill \end{array}`$ (A8)
where $`A=a/L_T`$ is a dimensionless width, the ratio of the wire width to the thermal length.

TABLE I. Details of the boundary conditions discussed in the text. The upper cut-off in momentum is given by $`Dq_a^2=1/\tau `$, so that $`p_0=a/2\pi \sqrt{D\tau }=\sqrt{M}(a/L_T)`$ where $`L_T=\sqrt{2\pi D/T}`$ is the thermal length and $`M=1/2\pi T\tau `$.

FIG. 1. The diagrams needed to calculate the first-order perturbative correction to the pair-propagator $`L(q,0)`$. (a) Definition of pair propagator, $`L(q,0)`$, in terms of the pair polarization function, $`P(q,0)`$, and the 4-point BCS interaction $`|\gamma |`$. (b) Zeroth-order (mean-field) polarization function, $`P_0(q,0)`$, in a dirty superconductor. The solid lines are electron Green functions; the dashed lines are impurity lines. (c) First-order perturbative correction to pair polarization function, $`\delta P(q,0)`$. The wiggly line is the screened Coulomb interaction.

FIG. 2. Diagrammatic equation for the pair amplitude matrix $`\mathrm{\Gamma }(ϵ_n,ϵ_l)`$. Block $`\gamma `$ is the BCS interaction. Block $`t\mathrm{\Lambda }`$ is the correction to the effective interaction caused by the interplay of Coulomb interaction and disorder. Block $`\mathrm{\Sigma }`$ is the correction to the bare Cooperon ladder caused by Coulomb interaction and disorder, which yields a negligible contribution and is thus ignored in our numerical calculations.

FIG. 3. Comparison of theory (solid lines) and experiment (symbols) for transition temperature, $`T_c`$, as a function of resistance per square, $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$, for various wire widths, $`a=250\AA ,350\AA ,580\AA ,1000\AA `$. The experimental data is that of Xiong et al. The theory used is the one we have argued to be correct – resummation technique solved by exact diagonalization using zero-gradient boundary conditions. The theory plots are, from top to bottom, for $`1000\AA `$, $`580\AA `$ $`350\AA `$ and $`250\AA `$ wires. Rescaling has been applied to take into account the fact that theory and experiment do not agree with each other in two dimensions, and that the transition temperature at zero resistance appears to depend upon width. The thermal length, $`L_T`$, has been taken to be $`300\AA `$ at $`T_{c0}=6.9K`$ throughout – this is the value quoted by Sharifi et al. Plotting the data in this form focuses on the suppression of transition temperature due to disorder, and we see that the agreement between theory and experiment is good, particularly for the $`1000\AA `$ and $`250\AA `$ wires.

FIG. 4. Comparison of theory (solid lines) and experiment (symbols) for the additional relative suppression of transition temperature, $`\mathrm{\Theta }(R_{\text{ }\text{ }\text{ }\text{ }\text{ }},a)`$, as a function of resistance per square, $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$, for various wire widths, $`a=250\AA ,350\AA ,580\AA ,1000\AA `$. Plotting the data in this form focuses on the suppression of transition temperature due to width, and we see that the agreement between theory and experiment is not particularly good. This may have been anticipated from the experimental data which shows an abrupt crossover between the 2D film and the $`1000\AA `$ wire, which lie on one curve, and the $`580\AA ,350\AA ,250\AA `$ wires, which lie on a different curve. Our theory, and indeed any theory, predicts a smooth crossover.

FIG. 5. Predicted transition temperature, $`T_c`$, as a function of applied perpendicular magnetic field, $`H`$, for wires of various widths, $`a=500\AA ,600\AA ,750\AA ,2D`$ film, and resistances per square, $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}=0\mathrm{\Omega },1000\mathrm{\Omega },2000\mathrm{\Omega }`$. Three physical effects can be seen in this figure: increasing $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$ reduces $`T_c`$ dramatically due to the increased Coulomb repulsion; decreasing $`a`$ increases the upper critical field, $`H_{c2}`$, as a larger field is needed to fit a flux quantum into the wire; decreasing $`a`$ at fixed $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$ causes a small additional suppression of $`T_c`$ at zero field due to increased fluctuation effects caused by the crossover from two to one-dimensional behavior. The last two effects mean that the $`T_c(H)`$ curves for wires with the same $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$ and different $`a`$ cross each other at non-zero field. Finally the shape of the $`T_c(H)`$ curve changes as width $`a`$ is altered due to the dimensional crossover. There is obviously a large amount of data in this figure, and to analyse the individual trends we will need to re-plot some of this data in a different form.

FIG. 6. Predicted dependence of the shape of the $`T_c(H)`$ curve on wire width, $`a`$, at a fixed resistance per square, $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}=1000\mathrm{\Omega }`$. We plot using dimensionless variables $`T_c(H)/T_c(H=0)`$ vs $`H/H_{c2}`$ in order to compare the curve shape. As width is reduced we see a crossover from the two-dimensional limit, where the magnetic pair-breaking rate $`1/\tau _HH`$, to the one-dimensional limit where $`1/\tau _HH^2`$. The crossover occurs when the wire width is of the order of the thermal length, $`L_T`$. For the wider wires the $`T_c(H)`$ curve joins the two-dimensional curve at the value of the magnetic field such that the magnetic length, $`L_H`$, is of order the width, $`a`$; for narrower wires this value of $`H`$ exceeds the two-dimensional value of $`H_{c2}`$, and the curves never meet. For ultranarrow wires, $`a<250\AA `$, the Zeeman effect becomes important. This causes the $`T_c(H)`$ curves to change shape again, moving slowly back in the direction of the 2D limit.

FIG. 7. Predicted dependence of the shape of the $`T_c(H)`$ curve on resistance per square, $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$, at a fixed wire width, $`a=1000\AA `$. We see a crossover from two-dimensional to one-dimensioanl behavior as $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$ is increased. This disorder-tuned crossover occurs only for wires of intermediate width; a similar plot for $`a=10000\AA `$ shows all data falling on the two-dimensional curve independent of $`R_{\text{ }\text{ }\text{ }\text{ }\text{ }}`$, whilst that for $`a=250\AA `$ shows all data falling on the one-dimensional curve.


text:
# Measurement of Light-Cone Wave Functions by Diffractive Dissociation

## Introduction

### The Pion Light-Cone Wave Function

The internal momentum distribution of valence quarks in hadrons are fundamental to QCD. . They are generated from the valence light-cone wave functions integrated over $`k_t<Q^2`$, where $`k_t`$ is the intrinsic transverse momentum of the valence constituents and $`Q^2`$, the total momentum transfer squared. Even though these amplitudes were calculated about 20 years ago, there have been no direct measurements until those reported here.

The pion wave function can be expanded in terms of Fock states:

$$\mathrm{\Psi }=a_1|q\overline{q}+a_2|q\overline{q}g+a_3|q\overline{q}gg+\mathrm{}.$$
(1)
Two functions have been proposed to describe the momentum distribution amplitude for the quark and antiquark in the $`|q\overline{q}`$ configuration. The asymptotic function was calculated using perturbative QCD (pQCD) methods , and is the solution to the pQCD evolution equation for very large $`Q^2`$ ($`Q^2\mathrm{}`$):

$$\varphi _{as}(u)=\sqrt{3}u(1u).$$
(2)
$`u`$ is the fraction of the longitudinal momentum of the pion carried by the quark in the infinite momentum frame. It should not be confused with $`x_{Bj}`$ which is not specific to a single Fock state. The antiquark carries a fraction ($`1u`$). Using QCD sum rules, Chernyak and Zhitnitsky  proposed a function that is expected to be correct for low Q2:

$$\varphi _{cz}(u)=5\sqrt{3}u(1u)(12u)^2.$$
(3)
As can be seen from eqns. 2 and 3 and from Fig. 1, there is a large difference between the two functions. Measurements of form factors are insensitive to the wave function as these quantities are derived by integrating over the wave function and interpretation of the results is model dependent . In this work we describe an experimental study that maps the momentum distribution of the valence $`|q\overline{q}`$ in the pion. This provides the first direct measurement of the pion light-cone wave function (squared). The concept of the measurement is the following: a high energy pion dissociates diffractively on a heavy nuclear target. This is a coherent process in which the quark and antiquark break apart and hadronize into two jets. If in this fragmentation process the quark momentum is transferred to the jet, measurement of the jet momentum gives the quark (and antiquark) momentum. Thus: $`u_{measured}=\frac{p_{jet1}}{p_{jet1}+p_{jet2}}.`$ From simple kinematics and assuming that the masses of the jets are small compared with the mass of the di-jets, the virtuality and mass-squared of the di-jets are given by: $`Q^2M_{DJ}^2=\frac{k_t^2}{u(1u)}`$ where $`k_t`$ is the transverse momentum of each jet. By studying the momentum distribution for various $`k_t`$ bins, one can observe changes in the apparent fractions of asymptotic and Chernyak-Zhitnitsky (CZ) contributions to the pion wave function.

The basic assumption that the momentum carried by the dissociating $`q\overline{q}`$ is transferred to the di-jets was examined by Monte Carlo (MC) simulations of the asymptotic and (CZ) wave functions (squared). The MC samples were allowed to hadronize through the LUND PYTHIA-JETSET model  and then passed through simulation of the experimental apparatus (described in the next section) to simulate the effect of unmeasured neutrals and other experimental distortions. In Fig. 1 the initial distributions at the quark level are compared with the final distributions of the detected di-jets. As can be seen, the qualitative features of the distributions are retained. The results of this analysis come from comparing the observed $`u`$-distribution to a combination of the distributions shown, as examples, on the right of Figure 1.

### The Color Transparency Effect

The Color-Transparency (CT) phenomenon is derived from the prediction that color fields cancel for physically small color singlet systems of quarks and gluons . This effect of color neutrality (or color screening) is expected to lead to the suppression of initial and final state interactions for a small sized system or point-like configuration (PLC) formed in a large angle hard process . Observation of CT requires that a PLC is formed and that the energies are high enough so that expansion of the PLC does not occur  (the frozen approximation). Under conditions of $`k_t>`$1.5 GeV/c, which translates to $`Q^210(\mathrm{GeV}/\mathrm{c})^2`$ and $`r0.1`$ fm, observation of these effects can be expected. Bertsch et al.  proposed that the small $`|q\overline{q}`$ component will be filtered by the nucleus. Frankfurt et al.  show that for $`k_t>1.5GeV/c`$ the interaction with the nucleus is completely coherent and $`\sigma (|q\overline{q}N`$di-jets $`N`$) is small. This leads to an $`A^2`$ dependence of the forward amplitude squared.

## Experimental Results

Fermilab experiment E791 recorded $`2\times 10^{10}`$ events from interactions of a 500 GeV/c $`\pi ^{}`$ beam with carbon and platinum targets. Details of the experiment are given in . Only about 10% of the E791 data was used for the analysis presented here. The data were analysed by selecting events in which 90% of the beam momentum was carried by charged particles. Jets were identifies using the JADE jet-finding algorithm . To insure clean selection of two-jet events, a minimum $`k_t`$ of 1.25 GeV/c was required and their relative azimuthal angle, which for pure di-jets should be $`180^{}`$ was required to be within $`20^{}`$ of this value.

Diffractive di-jets were identified through the $`e^{bq_t^2}`$ dependence of their yield ($`q_t^2`$ is the square of the transverse momentum transferred to the nucleus and $`b=\frac{<R^2>}{3}`$ where R is the nuclear radius). Figure 2 shows the $`q_t^2`$ distributions of di-jet events from platinum and carbon. The different slopes in the low $`q_t^2`$ coherent region reflect the different nuclear radii. Events in this region come from diffractive dissociation of the pion.

### The Pion Wave Function

For measurement of the wave function we used data from the platinum target as it has a sharp diffractive distribution and low background. We used events with $`q_t^2<0.015`$ GeV/c2. For these events, the value of $`u`$ was computed from the measured longitudinal momentum of each jet. A background, estimated from the $`u`$ distribution for events with larger $`q_t^2`$ was subtracted. The analysis was carried out in two windows of $`k_t`$: $`1.25\mathrm{GeV}/\mathrm{c}\mathrm{k}_\mathrm{t}1.5\mathrm{GeV}/\mathrm{c}`$ and $`1.5\mathrm{GeV}/\mathrm{c}\mathrm{k}_\mathrm{t}2.5\mathrm{GeV}/\mathrm{c}`$. The resulting $`u`$ distributions are shown in Fig. 3. In order to get a measure of the correspondence between the experimental results and the calculated light-cone wave functions, we fit the results with a linear combination of squares of the two wave functions (right side of Fig. 1). This assumes an incoherent combination of the two wave functions and that the evolution of the CZ function is slow (as stated in ).

The results of the fits are given in Fig. 3 and in Table 1 in terms of the coefficients $`a_{as}`$ and $`a_{cz}`$ representing the contributions of the asymptotic and CZ functions, respectively. The results for the higher $`k_t`$ window show that the asymptotic wave function describes the data very well. Hence, for $`k_t>`$1.5 GeV/c, which translates to $`Q^210(\mathrm{GeV}/\mathrm{c})^2`$, the pQCD approach that led to construction of the asymptotic wave function is reasonable. The distribution in the lower window is consistent with a significant contribution from the CZ wave function or may indicate contributions due to other non-perturbative effects.

The $`k_t`$ dependence of diffractive di-jets is another observable that can show how well the perturbative calculations describe the data. As shown in  assuming interaction via two gluon exchange and $`\varphi _{as}`$ would lead to $`\frac{d\sigma }{dk_t}k_t^6`$. The results, corrected for experimental acceptance, are shown in figure 4(a) fitted by $`k_t^n`$ for $`k_t>1.25`$ GeV with n = $`9.2\pm 0.4(stat)\pm 0.3(sys)`$ and $`\chi ^2/dof`$ = 1.0. This slope is significantly larger than expected. However, the region above $`k_t`$ 1.8 GeV/c can be fitted (Fig. 4(b)) with n = $`6.5\pm 2.0`$ with $`\chi ^2/dof`$ = 0.8, consistent with the predictions. This would support the evaluation of the light-cone wave function at large $`k_t`$ as due to one gluon exchange. The lower $`k_t`$ region can be fitted with the non-perturbative Gaussian function: $`\psi e^{\beta k_t^2}`$ , with $`\beta =1.78\pm 0.05(stat)\pm 0.1(sys)`$ and $`\chi ^2/dof`$ = 1.1. Model-dependent values in the range of 0.9 - 4.0 were used . These results are consistent with the measurements of the wave function that indicated noticeable non-perturbative effects up to $`k_t`$ 1.5 GeV/c.

### The Color Transparency Effect

To study the CT effect we measure the A-dependence of the diffractive di-jet yield. The coherence length is estimated using $`2p_{lab}`$ = 1000 GeV/c and $`M_J`$ 5 $`\mathrm{GeV}/\mathrm{c}^2`$. The result is $`l_c`$ 10 fm, larger than the platinum nuclear radius. In order to correct for experimental acceptance, we generate MC simulations of diffractive di-jets using the asymptotic wave function and di-jet masses of 4,5, and 6 GeV/c. The simulated coherent $`q_t^2`$ distributions of the di-jets represent the nuclear form factors of carbon (R=2.44 fm) and platinum (R=5.27 fm)  and the incoherent dissociation is simulated according to the nucleon radius (R=0.8 fm)  truncated at $`q_t^2<`$ 0.015. A combination of these simulations was used to fit the data (Fig. 2). We derive the numbers of producted di-jet events in the data for each target in three $`k_t`$ bins by integrating over the diffractive terms in the fits. Using the resulting yields and the known target thicknesses, we determine the ratio of cross sections for diffractive dissociation on platinum and carbon (the two targets were subjected to essentially the same beam flux). The exponents $`\alpha `$ are then calculated using the cross section dependence $`\sigma A^\alpha `$. The results are listed in Table 2 and compared with CT theoretical predictions . The results are consistent with those expected from color-transparency calculations and clearly inconsistent with $`\alpha `$ values for incoherent scattering observed in other hadronic interactions.

## The Photon Light-Cone Wave Function

The photon l.c. wave function can be described in a way similar to that of the pion except that it has two major components: the electromagnetic and the hadronic. Being a gauge field capable of piont-like coupling it has also a point bare-photon component. Consequently, the photon light-cone wave function in the pQCD regime can be expanded in terms of Fock states:

$`\psi _\gamma =a|\gamma _p+b|l^+l^{}+c|l^+l^{}\gamma +(othere.m.)`$
$`+d|q\overline{q}+e|q\overline{q}g+(otherhadronic)+\mathrm{}.`$ (4)
where $`|\gamma _p`$ describes the point bare-photon and $`|l^+l^{}`$ stands for $`|e^+e^{},|\mu ^+\mu ^{}`$ etc. The photon can be described as containing $`|q\overline{q}`$ components even in the non-pQCD regime of the vector-meson dominance model , in particular the $`\rho ,\omega ,\varphi `$ and $`J/\psi `$ mesons. There is also a $`|q\overline{q}`$ component with large relative transverse momentum that is connected to the point-like (direct) photon. The wave function of the photon is very rich: it can be studied for real photons, for virtual photons of various virtualities, for transverse and longitudinal photons and the hadronic component may be decomposed according to the quarks flavor. Strange quarks will reflect the strange content of the photon and charmed quarks will reflect the smaller sized charmed components.

The point-bare photon does not have internal structure, it can only Compton-scatter and will not be studied here. The other electromagnetic Fock states begin with $`|l^+l^{}`$ and continue to more complex systems. The interaction of these Fock states with the target is expected to be purely electromagnetic, as shown in figure 5(a). These were studied extensively through measurements of Bethe-Heitler leptoproduction  and of $`F_{2,QED}^\gamma `$ . Exclusive measurements of the $`|l^+l^{}`$ wave function will complement these inclusive measurements, may provide more tests of QED and will help define the analysis tools for the more complex hadronic structure.

The hadronic components interact according to the diagram of Fig. 5(b). The distribution amplitudes for these components, integrated over $`k_t`$ were calculated by Petrov et al.  using the instanton model and are shown in figure 6. Balitsky et al.  predict for real photons a distribution amplitude identical to the pion asymptotic function  (see Fig. 1).

While a variety of cross sections, form factors etc. depend on the light-cone wave function they are usually not sensitive to its structure. This is because it normally enters the calculations in an integrated form. Such are also the photon inclusive structure functions $`F_2^\gamma `$ . While this is an advantage when one needs calculations that do not depend on the internal structure of the photon, it is a disadvantage if we want to understand it. A differential measurement which will be sensitive to the $`u`$ and $`k_t`$ dependence of $`\varphi (u,Q^2)`$ will test the most fundamental description of the photon internal structure. Such measurements will make it possible to compare with the theoretical predictions and determine the regime of their validity.

The experimental program is presently in preparation and is based on diffractive dissociation of a real or virtual photon into di-leptons or di-jets. The electromagnetic component will be studied through measurements of the $`u`$ distribution for pure diffractive $`\mu ^+\mu ^{}`$ ($`e^+e^{}`$) elastic photoproduction and the $`k_t`$ distribution of these events. Similar measurements of diffractive hadronic di-jets will be used to study the hadronic component.

I would like to acknowledge the efforts of the E791 collaboration, of which I am a member, for the data presented in this work and the data analysis performed by my graduate student R. Weiss-Babai.


text:
# QCD corrections to top quark spin correlations at hadron colliders Report-no. DESY 00-120

## Discussion

A.P. Contogouris, Univ. of Athens
I am worried not by the method of your calculation, but by the importance of your numerical results. In particular the channel $`ggt\overline{t}(g)`$ can change them very much, and I am glad that you have already set up its calculation. Of course, $`q\overline{q}t\overline{t}(g)`$ is required to be complete.
A. Brandenburg
I agree. Note that the plots show results at the parton level for $`q\overline{q}t\overline{t}(g)`$.