text:
# Contents

## 1  Introduction

Channelling of a beam of charged particles in a bent monocrystal, is going to become a working tool for the next generation of accelerators . Therefore there is a need for a theoretical tool which describes a whole set of experimental data on channelling in the GeV range, and which also simulates the processes important for the future applications.

The output should be the distribution of exiting primary and secondary particles, the energy loss in crystal, and any other interesting quantity related to channelling. Since these processes are sensitive to the orientation, the simulation should track every particle through the crystal lattice computing the probability of any process as a function of the co-ordinates.

Historically, one of the independent discoveries of channelling, in the early sixties, was done in a Monte Carlo simulation of low-energy ($``$MeV) ions propagating in crystals . The very low ranges and the thin crystals used allowed the study of binary collisions of the incident ion with the atoms of the crystal. At GeV energy , crystals of a few centimetres in length are used, therefore tracking with binary collisions would take a considerable time. Instead, an approach with continuous potential introduced by J. Lindhard  can be used. In this approach one considers collisions of the incoming particle with the atomic strings or planes instead of separate atoms, if the particle is sufficiently aligned with respect to the crystallographic axis or plane. The validity of doing this improves with the increase of the particle velocity .

Besides the motion in the potential one must take into account the scattering. This makes it necessary to use either kinetic equations , or computer simulation  to transport particles through a crystal. The general feature of the methods described in Refs.  is the use of the diffusion approach which omits the single scattering acts. However, in Monte Carlo methods it is easy to include the single collisions with nuclei and electrons. Besides their influence on channelling, these close-encounter processes are the source of secondary particles emitted from the crystal. Moreover, such collisions with electrons provide interesting peculiarities of the energy loss spectra in aligned crystals. Here we describe the Monte Carlo program CATCH for the simulation of planar channelling, which does not use the diffusion approximation.

To simulate the nuclear collisions we use the LUND  routines; the CATCH program serves as a frame to provide the orientational dependence of these processes. Another useful feature of the CATCH program is the various imperfections of the crystal surface, incorporated in the simulation because of their essential role in the crystal-assisted beam extraction. Version 1.4 of CATCH allows variable curvature of the crystal, both longitudinal and transversal, in order to set, for instance, a twist of crystal faces.

## 2 Continuous potential

For the potential of the atomic plane we use the Molière approximation; details can be found e.g. in the review by D. Gemmel . The static-lattice potential is modified to take into account the thermal vibrations of the lattice atoms; this is done by integration over the Gaussian distribution of the atom displacement. Bending of the crystal has no effect on this potential. However, it causes a centrifugal force in the non-inertial frame related to the atomic planes. To solve the equation of motion in the potential $`U(x)`$ of the bent crystal, as a first approximation to the transport of a particle,

$$pv\frac{d^2x}{dz^2}=\frac{dU(x)}{dx}\frac{pv}{R(z)},$$
(1)
($`x`$ being the transversal, $`z`$ the longitudinal coordinate, $`pv`$ the particle longitudinal momentum and velocity product, $`R(z)`$ the local radius of curvature), we use the fast form of the Verlet algorithm :

$$x_{i+1}x_i=(\theta _i+0.5f_i\delta z)\delta z,$$
(2)
$$\theta _{i+1}\theta _i=0.5(f_{i+1}+f_i)\delta z$$
(3)
with $`\theta `$ for $`dx/dz`$, $`f`$ for the ‘force’, and $`\delta z`$ for the step. It was chosen over the other second order algorithms for non-linear equations of motion, such as Euler-Cromer’s and Beeman’s, owing to the better conservation of the transverse energy shown in the potential motion.

Figure 1 shows the simulated phase trajectories of protons on the plane $`(x,\theta )`$ in Si (111). Scattering was also included; as a result, there are no channelled protons near the atomic planes. With version 1.4 of CATCH one can specify in the input cards the planar geometry equivalent either to the Si(110) geometry (equidistant planes), or to the Si(111) one. It is easy to ‘build’ any other geometry by editing the program. The bending curvature and the crystal plane orientation can be arbitrary functions of spacial coordinates. Any data measured (or assumed) for the real crystal shape can be implemented in the simulation.

## 3  Scattering

Beam bending by a crystal is due to the trapping of some particles in the potential well $`U(x)`$, where they then follow the direction of the atomic planes. This simple picture is disturbed by scattering processes which could cause (as result of one or many acts) the trapped particle to come to a free state (feed out, or de channelling process), and an initially free particle to be trapped in the channelled state (feed in, or volume capture).

### 3.1  Scattering on electrons

Feed out is mostly due to scattering on electrons , because the channelled particles keep far from the nuclei. The mean energy loss in this scattering can be written as follows :

$$\frac{dE}{dz}=\frac{D}{\beta ^2}(0.5\mathrm{ln}\frac{2m_ec^2\beta ^2\gamma ^2}{I}\beta ^2\frac{\delta }{2}+0.5\rho _e(x)\mathrm{ln}\frac{T_{max}}{I}),$$
(4)
with $`D=4\pi N_Ar_e^2m_ec^2z^2\frac{Z}{A}\rho `$, $`z`$ for the charge of the incident particle (in units of $`e`$), $`\rho `$ for the crystal density, $`Z`$ and $`A`$ for atomic number and weight, and the other notation being standard .

The second part in the brackets is due to single collisions and depends on the local density $`\rho _e(x)`$ (normalized on the amorphous one) of electrons. The angle of scattering in soft acts can be computed as a random Gaussian with r.m.s. value

$$\theta _{rms}^2=\frac{m_e}{p^2}(\delta E)_{soft},$$
(5)
where $`(\delta E)_{soft}`$ is the soft acts, contribution in Eq. (4). The probability of the hard collision contributing to the second part of Eq. (4) is computed at every step. The energy transfer $`T`$ in such an act is generated according to the distribution function $`P(T)`$:

$$P(T)=\frac{D\rho _e(x)}{2\beta ^2}\frac{1}{T^2}.$$
(6)
The transverse momentum transfer $`q`$ is equal to

$$q=\sqrt{2m_eT+(T/c)^2}.$$
(7)
Its projections are used to modify the angles $`\theta _x`$ and $`\theta _y`$ of the particle.

### 3.2  Scattering on nuclei

The multiple Coulomb scattering on nuclei is computed by the approximation Kitagawa–Ohtsuki :

$$\theta _{nucl,sc}^2=\theta _{sc}^2_{amorph}\rho _n(x)$$
(8)
i.e. the mean angle of scattering squared is proportional to the local density of nuclei $`\rho _n(x)`$. The density function is Gaussian with r.m.s. value $`u`$ being the amplitude thermal vibration of the atom.

The probability of nuclear collision, proportional to $`\rho _n(x)`$, is checked at every step. If such a collision succeeds (then the flag ISTATU is set to –1), the LUND routine responsible for the event generation may be called. This scheme may be used for the description of any close-encounter process of interest; then one should define the collision length (ADLI in the input cards) properly. In version 1.4 of CATCH the secondary particles are normally not tracked through the crystal lattice, but their exit angles can be modified due to multiple scattering in the rest part of the crystal. In principle, a full-scale tracking of secondary particles (and even their products) is possible via the construction of the appropriate main program (description below) with loops and trees.

## 4  Crystal imperfection

The complete description of the program usage follows below. Here we show how the crystal imperfection is implemented.

### 4.1  Variable curvature and twist

If KRADV = 0 is set in the input cards, the transverse curvature is set to 0. The longitudinal curvature is 0 for the longitudinal coordinate $`z`$ between 0 and STRAIT (length of the forthcoming straight part), and is constant equal to 1/RADIUS for $`z`$ between STRAIT and STRAIT + DLINA (hence DLINA means the length of the bent part ).

If KRADV = 1 is set in the input cards, the curvature at any point must be defined by the user in function CURVAT. The total length of the crystal is STRAIT + DLINA in any case.

If KRADVY = 0 is set in the input cards, the orientation of the crystal planes at both entry and exit faces is constant. If KRADVY = 1 is set in the input cards, these orientations are computed according to the functions TCFROY and TCENDY, for the front face and for the back face respectively. These functions must be supplied by the user.

### 4.2  Surface effects

In the case of beam extraction from the accelerator, extremely small impact parameters are possible, making the surface effects essential. The particle entering the crystal very close to its edge, can suffer from various additional factors:

* miscut angle (between the atomic planes and the surface),
* roughness (i.e. non flatness) of the surface,
* possible amorphous layer,
* bent surface.

Therefore one must pay particular attention to the near-surface tracking, where the particule is entering and leaving the cristal materi al (due to the roughness, holes and bend), both coherent and non-coherent scattering in this peculiar region, bending in short channels, and so on. The surface effects mentioned above are simulated in the current version of CATCH. The roughness is expressed by a periodical function, $`a\mathrm{sin}(z/\lambda )`$, where $`a`$ is the amplitude of ‘bumps’ and $`\lambda `$ is their periodicity. The ‘rough’ crystalline material can be superimposed by a uniform amorphous layer. The position of surfaces is computed at every step in accordance with bending with the variable curvature. To set up these imperfections, simply give non-zero values to the corresponding parameters in the input cards.

## 5  Usage of routines

Figure 2 shows an example of usage of the CATCH routines. The subroutine XPREP asks for a filename with a problem description, and performs some preliminary procedures. The subroutine XCATCH takes the entry values of $`x`$, $`\theta _x`$, $`y`$, $`\theta _y`$ and $`E`$ and returns the exit values. The entry/exit values of X, Y, XP, YP and ENERGY are double precision . The units in CALL…(..) are radians and metres; the energy is returned modified according to the E loss in crystal; hitting/missing the crystal is checked, the flag for that is ISTATU (returned) :

* ISTATU = 0, crystal missed,
* ISTATU = 1, or –1, crystal hit,
* ISTATU $`<`$ 0, there was nuclear interaction.

In the case of interaction one may call the necessary routines to produce the secondary particles and to get their exit parameters modified for multiple scattering. The XPOST can be called to save the histograms, to fit the data, etc. Figure 3 shows an example of the input file. The input data is listed below with comments when necessary.

* title
* KGEOM, geometry type: either 110 (equidistant planes) or 111 (planes like those in the Diamond(111) lattice)
* KTRACE, trace level:

+ if $``$ 0, no output during tracking
    + if $`=`$ 1, some output at every step of tracking
    + if $``$ 2, all current parameters are printed at every step
* KRADV, switch on/off $`R=R(z)`$ option; with KRADV = 1 one must supply FUNCTION CURVAT, evaluating the curvature $`1/R`$ for every point; with KRADV = 0 the curvature is 1/RADIUS.
* KRADVY, switch on/off face twist; with KRADVY=1 one must supply functions TCFROY and TCENDY, evaluating the orientation of the planes for any point at the crystal faces; with KRADVY = 0 that orientation is constant.
* I1, I2, initial numbers for pseudo-random generator RAN(I1,I2)
* DZ, step size, microns. Try different DZ to be sure the result is independent of the step size. A step of the order of 1 $`\mu `$m is suggested for 450 GeV/$`c`$ protons in Si(110). The reasonable value of DZ should be scaled like $`\sqrt{E}`$ (like the oscillation period in channelling)
* ENERGY, GeV
* PMASS, mass of the beam particle, (GeV)
* RADIUS, in cm, for the bent part of the crystal. When KRADV = 1 is set, the RADIUS is used to define the histogram windows only.
* DLINA, length of the bent part of the crystal (cm). The whole length of the crystal is STRAIT + DLINA. When KRADV = 0 is set, the forthcoming part of the length STRAIT is assumed to be straight. The following part of the th DLINA is assumed to be bent with R = RADIUS.
* STRAIT, length of the forthcoming straight part, (cm)
* HDL, interplanar half-spacing, for (110) or (111)L, in Å
* HDS, interplanar half-spacing for (111)S, in Å. Must be 0 for KGEOM=110
* U, thermal vibration amplitude of the lattice atoms,(Å)
* DENS, crystal density, (gram/cm3)
* RDLI, radia ion length, (cm)
* ADLI, absorption length, (cm)
* ZN, charge of the nucleus
* AN, atomic weight
* X0XTL, the crystal entry-face left-edge $`x`$-coordinate, (cm)
* X1XTL, the crystal entry-face right-edge $`x`$-coordinate, (cm)
* Y0XTL, the crystal entry-face left-edge $`y`$-coordinate, (cm)
* Y1XTL, the crystal entry-face right-edge $`y`$-coordinate, (cm)
* THXTAL, crystal angle $`\theta _x`$, ($`\mu `$rad)
* THYTAL, crystal angle $`\theta _y`$, ($`\mu `$rad)
* SKINTH, amorphous skin thickness, (cm)
Some amorphous skin with uniform thickness SKINTH can be superimposed over the surface.
* ROUGH, (cm) the roughness (= non-flatness) amplitude of every surface. It is described by some periodical function, therefore you should supply in the input both its amplitude ROUGH and period ROUGPD. It is assumed that material under the rough surface is crystalline.
* ROUGPD, surface roughness period, (cm)
* SKIMOZ, ($`\mu `$rad) There is the same skin, SKINTH + ROUGH, at the entry face. If you wish, you can also set some angular spread SKIMOZ of the crystal lattice in the entry-face skin (in this skin only!). The SKIMOZ means the sigma of Gaussian distribution.
* DISCUT, miscut angle, ($`\mu `$rad) Note that its sign is essential; when your impact is near the left edge of the crystal, the negative sign (DISCUT$`<`$0) is preferable.

The position of all surfaces in space, as a function of $`\theta _x`$ and $`\theta _y`$ misalignment of the crystal (THXTAL and THYTAL in the input cards), as a function of bending with variable radius, as well as a function of the surface ’rough structure’, is computed at every step.

All variables read from the input cards are placed into COMMON blocks, and hence you can vary them from the main program during the execution (if you wish to scan or to optimize). Every particle (both channelled and non-channelled) is tracked in the crystal (X, Y are changed) and can leak out through any surface. Near the rough surface it is even possible to leak out and be caught again many times, i.e. the particle traverses sequentially the crystal bumps and vacuum (or amorphous skin) between. Hence, for particles touching any surface one has both Coulomb multiple and coherent collisions.

## 6  Acknowledgements

The hospitality of the CERN SL Division, in particular of W. Scandale, and discussions with users of CATCH at Protvino, Pisa, and CERN have fruitfully contributed to this work.


format: text

text:
# A Simulation of the ’t Hooft Model at finite 𝑁_𝑐 with the Overlap Dirac OperatorPresented by F. Berruto. Work supported in part by grant DE-FG02-91ER40676 and by INFN.

## 1 The ’t Hooft Model with Overlap Fermions

The ’t Hooft model action in the overlap regularization reads

$`S`$ $`=`$ $`\beta {\displaystyle \underset{x,\mu <\nu }{}}\left[1{\displaystyle \frac{1}{2N_c}}\text{Re}\text{Tr}U_{\mu \nu }(x)\right]`$
$`+`$ $`{\displaystyle \underset{x,y}{}}\overline{\psi }(x)\left[(1{\displaystyle \frac{ma}{2}})D_{xy}+m\delta _{xy}\right]\psi (y)`$
where $`U_{\mu \nu }`$ is the usual Wilson plaquette, $`a`$ is the lattice spacing and $`\psi `$ are $`N_f`$ fermion fields with mass $`m`$ and $`\beta =2N_c/(ag)^2`$. $`D=(1+V)/a`$ is the massless Neuberger-Dirac operator, $`V`$ being the unitary component of $`D_W1/a`$, where $`D_W`$ is the usual Wilson-Dirac operator.

For $`m=0`$ the action (1) is invariant under a continuous symmetry $`\delta \psi =\gamma _5(1aD)\psi ,\delta \overline{\psi }=\overline{\psi }\gamma _5`$, i.e. a form of chiral symmetry that holds at finite lattice spacing. The $`N_f=2`$ models exhibit an isotriplet of particles ($`\pi `$), whose masses vary with the quark mass as

$`M_\pi ^2`$ $`=`$ $`2\sqrt{{\displaystyle \frac{C_t\pi }{3}}}m+\mathrm{},U(N_c)\text{ models}`$ (2)
$`M_\pi ^2`$ $`=`$ $`{\displaystyle \frac{9}{\pi }}(2^7C_t)^{\frac{1}{3}}({\displaystyle \frac{e^\gamma }{\pi }})^{\frac{4}{3}}m^{\frac{4}{3}}+\mathrm{},SU(N_c)\text{ models}`$ (3)
Eq.(2) has been derived for $`N_c\mathrm{}`$ while Eq.(3) also involves a semiclassical WKB approximation.

In the $`U(N_c)`$ models in the chiral limit $`M_\eta ^{}^2=N_fg^2/\pi `$. The Witten-Veneziano formula for the $`U(N_c)`$ models gives

$$M_\eta ^{}^2=\frac{4N_f}{f_\pi ^2}\frac{Q^2}{V}$$
(4)
## 2 The meson spectrum

We performed extensive simulations with the exact Neuberger-Dirac operator on a $`18\times 18`$ lattice at $`C_t=4/3`$, i.e. $`\beta =6,13.5,24`$ for $`N_c=2,3,4`$ respectively. We computed $`M_\pi `$ for $`m=0.04,0.05,0.06,0.07,0.08,0.1`$ so that the relation $`N_xM_\pi >4`$ holds. The effects of dynamical fermions have been included by weighting the observables with the appropriate power of the fermion determinant. The smallness of the lattice warrants this procedure (cfr. ), which would lead to unacceptable variance on larger systems.

We computed the meson propagators projected over zero momentum from the fermion propagators and extracted the meson masses in a standard manner. The errors have been estimated with the jacknife analysis.

In table (1) we report the pion masses for the $`U(N_c)`$ and $`SU(N_c)`$ models with two flavors of dynamical fermions and in the quenched approximation.

The quenched results get closer to the unquenched two flavors results when $`N_c`$ gets larger, as one would expect since the large $`N_c`$ limit corresponds to discarding the internal quark loops. For the $`U(N_c)`$ case Eq.(2) motivated us to do the fit

$$\frac{M_\pi ^2}{g^2N_c}=A+B\frac{m}{g\sqrt{N_c}}$$
(5)
We obtained $`A=0.012(17),B=1.57(19)`$ for $`N_c=2`$, $`A=0.006(4),B=1.53(4)`$ for $`N_c=3`$ and $`A=0.006(8),B=1.66(8)`$ for $`N_c=4`$.

For the $`SU(N_c)`$ models, following Eq.(3), we performed the fit

$$\frac{M_\pi }{g\sqrt{N_c}}=A+B(\frac{m}{g\sqrt{N_c}})^{\frac{2}{3}}$$
(6)
We obtained $`A=0.022(26),B=1.524(120)`$ for $`N_c=2`$, $`A=0.045(43),B=1.548(166)`$ for $`N_c=3`$ and $`A=0.009(40),B=1.788(11)`$ for $`N_c=4`$. In the $`SU(N_c)`$ case we performed also the fit $`M_\pi /g\sqrt{N_c}=C(m/g\sqrt{N_c})^\gamma `$ and we obtained $`\gamma `$ compatible with $`2/3`$ for $`N_c=2,3,4`$. Figures(1,2) give evidence of universality in $`N_c`$ for the pion masses.

In the $`U(N_c)`$ models we computed also $`M_\eta ^{}`$. The $`\eta ^{}`$ propagator is given by the difference between the connected and disconnected terms in the correlator and, because of cancellations, the errors are larger than in the $`\pi `$ case. Therefore, we resorted the method proposed in  to compute the $`\eta ^{}`$ mass. The method exploits the quenched two-loop disconnected $`\mathrm{\Gamma }_q^{2\text{loop}}`$ and quenched one-loop connected $`\mathrm{\Gamma }_q^{1\text{loop}}`$ contributions to the $`\eta ^{}`$ propagator:

$$M_\eta ^{}^2=2M_\pi \underset{t\mathrm{}}{lim}\frac{\mathrm{\Gamma }_q^{2\text{loop}}(t)}{|t|\mathrm{\Gamma }_q^{1\text{loop}}(t)}$$
(7)
As one can see from table (2) $`M_\eta ^{}^2`$ computed using Eq.(7) and $`N_f=1`$ compares remarkably well with the analytic result.

An alternative calculation of $`M_\eta ^{}`$ can be done by using the Witten-Veneziano formula. In the quenched approximation we computed $`f_\pi =0.774(3),0.956(3),1.115(5)`$ for $`N_c=2,3,4`$ respectively. Since in our calculations we diagonalize the full Neuberger-Dirac operator, we could determine the number of zero modes, and thus the topological charge, for each individual gauge configuration. We thus obtained $`Q^2/V=0.0258(16),0.0298(23),0.0319(37)`$ for $`N_c=2,3,4`$ respectively. (Notice that with increasing $`N_c`$ these values get closer and closer to the analytic value $`Q^2/V=C_t/4\pi ^2=0.0337`$. Building on the trade-off between spatial and internal degrees of freedom that prompted the introduction of the $`N_c\mathrm{}`$ single plaquette model , one could heuristically argue that the theory should be less affected by finite size effects increasing $`N_c`$. A rigorous justification of this argument, however, would require a detailed analysis of discretization and finite size effects that goes beyond the scope of our investigation.) From $`f_\pi `$ and $`Q^2/V`$, we use Eq.(4) to calculate $`M_\eta ^{}`$. The results, also in table (2), are consistent for $`N_c=3`$ and 4 with those obtained with the method of Ref.  and with the analytical formula.


text:
# A fractional-spin phase in the power-law Kondo model

## I Introduction

The low-temperature physics of the standard Kondo model, describing a single magnetic impurity embedded in a metal, is by now well understood . For antiferromagnetic coupling between the impurity and the conduction electron spins the effective interaction grows with decreasing temperature. The low-energy behavior is completely determined by a single energy scale, the Kondo temperature $`T_\mathrm{K}`$, and the impurity spin is fully quenched in the low-temperature limit, $`TT_\mathrm{K}`$. In contrast, for ferromagnetic coupling the effective interaction decreases upon lowering the temperature, leaving the impurity spin essentially decoupled from the environment.

The standard Kondo picture has to be revised if the conduction band density of states (DOS) is not constant near the Fermi level. This is the case in systems with a power-law DOS $`\rho (ϵ)|ϵ|^r`$. Exponents $`r>0`$ lead to a vanishing DOS at the Fermi level – such a pseudogap DOS arises in one-dimensional interacting systems, in certain zero-gap semiconductors, and in systems with long-range order where the order parameter has nodes at the Fermi surface, e.g., $`p`$\- and $`d`$-wave superconductors ($`r=2`$ and 1). The pseudogap Kondo problem has attracted a lot of attention during recent years . These studies show the existence of a zero-temperature boundary phase transition at a critical antiferromagnetic Kondo coupling, $`J_c`$, below which the impurity spin is unscreened even at lowest temperatures. A comprehensive discussion of possible fixed points and their thermodynamic properties has been given by Gonzalez-Buxton and Ingersent  based on the numerical renormalization group (NRG) approach.

In this work we investigate a power-law Kondo model for exponents $`1<r<0`$ which has not been explored before. Negative $`r`$ implies a diverging low-energy DOS, corresponding to a critical or van-Hove singularity at the Fermi level.

The Kondo Hamiltonian for a spin-1/2 impurity can be written as $`H=H_{\mathrm{band}}+H_{\mathrm{int}}`$, with

$$H_{\mathrm{int}}=J𝐒𝐬_0+Vc_{0\sigma }^{}c_{0\sigma }$$
(1)
and $`H_{\mathrm{band}}=_{𝐤\alpha }ϵ_𝐤c_{𝐤\alpha }^{}c_{𝐤\alpha }`$ in standard notation, $`𝐬_0=_{\mathrm{𝐤𝐤}^{}\alpha \beta }c_{𝐤\alpha }^{}\sigma _{\alpha \beta }c_{𝐤^{}\beta }`$ is the conduction band spin operator at the impurity site $`𝐫_0=0`$, and $`𝐒`$ denotes the impurity spin operator. For simplicity we will use a conduction band with a symmetric density of states, $`\rho (ϵ)=\rho _0|ϵ/D|^r`$ for $`|ϵ|<D=1`$, $`\rho _0=(1+r)/(2D)`$. Particle-hole asymmetry is introduced via the potential scattering term $`V`$.

Our findings (Fig. 1) can be summarized as follows: Antiferromagnetic Kondo interactions in the particle-hole symmetric model ($`V=0`$) always flow to a strong-coupling fixed point with a singlet ground state (SSC fixed point) and complete screening of the impurity. The ferromagnetic sector has very rich behavior: large interactions flow to a ferromagnetic strong-coupling fixed point with a triplet ground state (TSC). Even more interestingly, for $`V=0`$ there exists a stable intermediate-coupling fixed point (FS). Its most remarkable property is a Curie response corresponding to a fractional spin,

$$\chi _{\mathrm{imp}}(T0)=\frac{𝒞(r)}{T},$$
(2)
where $`\chi _{\mathrm{imp}}`$ is the total impurity contribution to the uniform susceptibility, and $`𝒞`$ is a universal, irrational number which depends only on the DOS exponent $`r`$. (The $`V0`$ behavior is quite complicated as well and will be discussed below.)

Such fractional-spin states at intermediate coupling have so far only been found at infrared unstable (critical) fixed points, namely in the pseudogap Kondo model ($`r>0`$)  and for impurities coupled to quantum-critical magnets . We note that “exotic” states corresponding to stable intermediate-coupling fixed points also occur in multichannel Kondo models ; here, however, the leading term in the impurity susceptibility does not follow a Curie law.

## II Weak-coupling analysis

We start our considerations with the standard weak-coupling renormalization group treatment (poor man’s scaling ), here modified for the power-law Kondo model.

Integrating out a shell of high-energy conduction electrons gives the one-loop renormalization of the Kondo interaction of order $`J^2`$, rescaling of energies and band cut-off $`\mathrm{\Lambda }`$ leads to an additional renormalization of all couplings proportional to $`r`$ . Expressed in $`\beta `$ functions for the dimensionless running couplings $`j=\rho _0J`$ and $`v=\rho _0V`$ we have

$`{\displaystyle \frac{dj}{d\mathrm{ln}\mathrm{\Lambda }}}=j(rj)+𝒪(j^3),{\displaystyle \frac{dv}{d\mathrm{ln}\mathrm{\Lambda }}}=rv.`$ (3)
For $`r=0`$ the potential scattering term $`v`$ remains unchanged whereas positive (negative) Kondo coupling $`j`$ flows to $`\mathrm{}`$ (zero) upon decreasing the cut-off. In the pseudogap case, $`r>0`$, small $`j`$ of either sign renormalizes to zero. For small positive $`r`$ one can deduce a critical fixed point at $`J_c=r/\rho _0`$, $`V_c=0`$, corresponding to the transition between a local moment phase and an antiferromagnetic strong coupling phase in the pseudogap model .

Now we turn to $`r<0`$: here small $`j`$ grows for both signs of the Kondo coupling. In analogy to the above, on the ferromagnetic side for small negative $`r`$, we can predict a fixed point at $`J^{}=|r|/\rho _0`$ which is now infrared stable (for $`V=0`$). The properties of this novel intermediate-coupling fixed point are perturbatively accessible in a double expansion in $`j`$ and $`r`$ – as known from the theory of critical phenomena the properties are determined by the universal value of the running coupling constant $`j`$, and are completely different from the properties of usual stable phases where the renormalized couplings are either zero or infinite. In particular, local spin correlations show a power law decay with some non-trivial, $`r`$-dependent exponent, the ground-state entropy has a universal, non-zero value, and the total impurity-induced susceptibility shows fractional Curie behavior (2). In principle the calculation of these properties can be done in renormalized perturbation theory up to high orders; this does, however, not provide information about the fixed point structure for finite negative $`r`$ nor about the flow of large initial couplings $`J`$.

## III Numerical results

To verify the above picture and investigate the strong-coupling behavior, we have performed extensive studies of the power-law Kondo model using the NRG technique . From the flow of the NRG energy levels we have deduced a number of (meta)stable fixed points, with properties listed in the following:

(LM) The symmetric local-moment fixed point corresponds to the system with a decoupled impurity, $`J=V=0`$. The $`T=0`$ limit of the susceptibility is $`T\chi _{\mathrm{imp}}=1/4`$, the entropy $`S_{\mathrm{imp}}=\mathrm{ln}2`$. This fixed point is unstable both w.r.t. to finite $`J`$ and $`V`$, which follows from (3). The instability against potential scattering can be easily understood: finite $`V`$ leads to a pole in the scattering T matrix which creates a local (anti)bound state, leading to the ALM behavior below.

(ALM) An additional asymmetric local-moment fixed point exists, corresponding to $`J=0`$ and $`|V|=\mathrm{}`$; here the conduction electron site next to the impurity is either doubly occupied or empty. The impurity thermodynamic properties are similar to the LM fixed point.

(SSC) The singlet strong-coupling fixed point is reached for positive $`J`$ and small $`|V|`$, it corresponds to $`J=\mathrm{}`$, $`V=0`$ and displays a fully quenched spin with $`T\chi _{\mathrm{imp}}=0`$, $`S_{\mathrm{imp}}=0`$. (The leading term in $`\chi _{\mathrm{imp}}`$ has the form $`T^r`$, i.e., the impurity susceptibility itself vanishes in the zero-temperature limit.) In contrast to the power-law model with $`r>0`$, this fixed point is stable w.r.t. particle-hole symmetry breaking.

(TSC) A triplet strong-coupling fixed point is reached for large negative $`J`$, its properties follow from $`J=\mathrm{}`$, $`V=0`$ as $`T\chi _{\mathrm{imp}}=2/3`$, $`S_{\mathrm{imp}}=\mathrm{ln}3`$. This fixed point is also stable against finite $`V`$.

(FS) The new, fractional-spin fixed point exists for small $`r<0`$ and is reached for small $`J<0`$, $`V=0`$. It features universal, $`r`$-dependent values of $`T\chi _{\mathrm{imp}}`$ and $`S_{\mathrm{imp}}`$; it is unstable against non-zero potential scattering.

From the NRG results we can construct the flow diagrams shown in Fig. 1, and we can obtain various thermodynamic quantities as function of $`T`$ . In Fig. 2 we display the temperature dependence of $`T\chi _{\mathrm{imp}}`$ for different Kondo couplings $`J`$ and a fixed DOS exponent $`r=0.2`$, which puts the model in the intervall $`\overline{r}^{}<r<0`$, with the flow diagram depicted in Fig. 1a. The particle-hole symmetric case $`V=0`$ is shown in Fig. 2a. Positive values of $`J`$ lead to a fully quenched spin in the low-temperature limit. The “Kondo” temperature characterizing the flow to strong coupling depends in a power-law fashion on $`J`$, $`T_KD(\rho _0J)^\alpha `$, with the leading term being $`\alpha =1/r`$. The behavior is in contrast to the exponential dependence $`T_K(J)`$ in the metallic $`r=0`$ case, reflecting the fact that for $`r<0`$ antiferromagnetic $`J`$ is a relevant rather than a marginally relevant perturbation of the local-moment fixed point.

Turning to negative values of $`J`$, we see that all small ferromagnetic values of the initial coupling yield a flow to the predicted intermediate-coupling fixed point (FS), characterized by a universal Curie behavior of the impurity susceptibility, here $`T\chi _{\mathrm{imp}}0.30`$. Interestingly, the basin of attraction of this fixed point does not extend to arbitrarily large $`|J|`$: large ferromagnetic couplings flow to a different stable fixed point with $`T\chi _{\mathrm{imp}}=2/3`$ – this is the triplet strong-coupling (TSC) fixed point. The boundary quantum phase transition between the fractional-spin phase and the triplet phase appears to be continuous, and consequently there is a critical particle-hole symmetric fixed point (SCR, thick line in Fig. 2a), separating the fractional-spin and triplet strong-coupling regimes.

Switching on potential scattering, Fig. 2b, we find that the local-moment and fractional-spin fixed points are unstable w.r.t. particle-hole symmetry breaking, and both yield flows towards the asymmetric local moment (ALM) fixed point. (Note, however, that $`T\chi _{\mathrm{imp}}=1/4`$ for both the symmetric and asymmetric local-moment fixed points.) In contrast, both strong-coupling fixed points as well as the critical fixed point (SCR) are stable w.r.t. finite $`V`$ (Fig. 1a). On the antiferromagnetic side, finite $`V`$ suppresses Kondo screening, leading to another non-trivial transition between the singlet strong-coupling (SSC) and asymmetric local-moment (ALM) phases, with an associated particle-hole asymmetric critical fixed point at positive $`J`$ (AF-CR, here $`T\chi _{\mathrm{imp}}0.18`$ in Fig. 2b).

Our numerical investigation shows that the described fixed point structure on the ferromagnetic side changes for more negative $`r`$ values. With decreasing $`r<0`$, first the symmetric critical fixed point between the triplet strong-coupling (TSC) and fractional-spin (FS) phases becomes unstable w.r.t. broken particle-hole symmetry, and a second critical fixed point (ACR) with finite asymmetry appears. This means that at a value $`r=\overline{r}^{}`$ the critical fixed point splits, and the transition between triplet strong-coupling (TSC) and asymmetric local-moment phases (ALM) is now governed by the new asymmetric critical point ACR. The resulting renormalization group flow is shown in Fig. 1b. (Strictly speaking, each described fixed point at finite $`|V|`$ represents a pair of fixed points, which are trivially related by a particle-hole transformation.)

Upon further decreasing $`r`$, the ferromagnetic particle-hole symmetric critical fixed point SCR moves to smaller couplings, towards the stable FS fixed point. At a certain value $`r=\overline{r}_{\mathrm{max}}`$ both fixed points meet and disappear (!). This implies that for $`r<\overline{r}_{\mathrm{max}}`$ and $`V=0`$ any negative value of $`J`$ flows to strong coupling, i.e., leads to a triplet state between the impurity and the conduction electrons. At finite particle-hole asymmetry, there is still a transition between a local moment and ferromagnetic strong coupling phase, and the critical behavior is governed by a single critical fixed point at finite asymmetry, ACR, see Fig. 1c. The above changes do not influence the structure of the flow diagram on the antiferromagnetic side.

The Curie part of the impurity susceptibility, $`T\chi _{\mathrm{imp}}`$, of the various intermediate-coupling fixed points is shown in Fig. 3. This illustrates nicely the splitting of the unstable ferromagnetic fixed point at $`r=\overline{r}^{}`$ and the collapse of the two fixed points at $`r=\overline{r}_{\mathrm{max}}`$.

The ferromagnetic side of the flow diagrams in Fig. 1 shows some superficial similarity to the antiferromagnetic $`r>0`$ situation : Also in this case, two critical fixed points coexist over a certain range of DOS exponents, $`r^{}<r<r_{\mathrm{max}}`$. However, there are crucial differences between the present $`r<0`$, $`J<0`$ case and the $`r>0`$, $`J>0`$ regime studied in Ref. : First, the novel stable intermediate-coupling fractional-spin phase (FS) has no counterpart for $`r>0`$. Second, in the $`r<0`$ case the weak (strong) coupling fixed points are unstable (stable) w.r.t. to particle-hole symmetry breaking; for $`r>0`$ this is reversed.

Fig. 4 summarizes the phase diagram for the particle-hole symmetric power-law Kondo model for both signs of the exponent $`r`$ and both signs of $`J`$. The line through the origin describes the fixed points which can be deduced from the weak-coupling equations (3) – note that the line of critical fixed points for $`r>0`$ has the same initial slope as the line of stable fixed points for $`r<0`$.

In the vicinity of each critical point one can define an energy scale $`T^{}`$, which vanishes at the transition, and defines the crossover energy above which quantum-critical behavior is observed . We note, however, that on the ferromagnetic side simple one-parameter scaling (as function of $`T/T^{}`$ and $`\omega /T^{}`$) can only be expected for $`V=0`$ and $`\overline{r}_{\mathrm{max}}<r<0`$ and for $`V0`$ and $`1<r<\overline{r}^{}`$. In contrast, in the asymmetric case for $`\overline{r}^{}<r<0`$ the $`J<0`$ critical behavior is governed by the symmetric multicritical point (SCR, Fig. 1a), leading to two distinct energy scales vanishing at the transition. (This situation is similar to the behavior for small positive $`r`$ .) We have numerically checked the scaling behavior of $`\chi _{\mathrm{imp}}`$ and of the T matrix, details will appear elsewhere.

## IV Conclusions and outlook

In this paper, we have examined a novel part of the phase diagram of the power-law Kondo model, namely the behavior occuring for a bath density of states diverging as $`|\omega |^r`$ with exponent $`1<r<0`$. We have found rich physics for both antiferromagnetic and ferromagnetic Kondo coupling. Particularly interesting is a stable phase with fractional ground state spin, i.e., a Curie susceptibility (2) with a non-trivial universal coefficient $`𝒞(r)`$, which also shows a non-zero ground state entropy. This phase occurs for ferromagnetic Kondo coupling and corresponds to an intermediate-coupling fixed point in the renormalization group sense. It cannot be understood in a weak- or strong-coupling picture in terms of the coupling between impurity and conduction electron degrees of freedom. Instead, it shares many properties of critical fixed points, namely the absence of quasiparticle excitations, scale invariance of the low-energy correlations, and non-trivial power laws in response functions like the dynamical susceptibility. A preliminary analysis of the NRG levels at the fractional-spin fixed point shows that it cannot be described in terms of non-interacting particles – in contrast to, e.g., the two-channel Kondo fixed point which allows for a Majorana fermion description .

In closing, we mention that the discussed power-law Kondo model is not only of purely theoretical interest: Besides describing “real” magnetic impurities – where a power law DOS can occur due to band structure effects such as van-Hove singularities – effective single-impurity models arise in the context of lattice models in the limit of infinite dimensions (dynamical mean-field theory (DMFT) ). Here, the DOS of the effective embedding medium is determined from a self-consistency condition involving both the impurity spectrum and the bare DOS of the conduction band. The interplay of these quantities can e.g. generate a power-law DOS for the effective medium at a critical point of the lattice model. Particularly interesting in this context is the so-called extended DMFT scheme  where the influence of spin fluctuations can drive the local impurity problem critical, which in turn will generate a power-law spectral density of the effective fermionic medium in extended DMFT.

Furthermore, we note that the low-energy physics of Kondo models is related to that of dissipative two-level systems, i.e., a Kondo model with a diverging power-law DOS corresponds to a spin-boson model with a certain non-ohmic bath. Therefore, applications in mesoscopics, like two-level systems coupled to a noisy environment, or in glass physics may be found.

These prospects, together with the possibility of non-Fermi-liquid physics in lattice realizations of the power-law Kondo physics, e.g., in extended DMFT, will be investigated in the future.

###### Acknowledgements.

We thank M. Glossop, S. Kehrein, and Th. Pruschke for valuable discussions. This research was supported by the DFG through SFB 484.


text:
# Quasiparticle Electronic structure of Copper in the GW approximation

## Abstract

We show that the results of photoemission and inverse photoemission experiments on bulk copper can be quantitatively described within band-structure theory, with no evidence of effects beyond the single-quasiparticle approximation. The well known discrepancies between the experimental bandstructure and the Kohn-Sham eigenvalues of Density Functional Theory are almost completely corrected by self–energy effects. Exchange–correlation contributions to the self–energy arising from 3s and 3p core levels are shown to be crucial.

Experimental techniques for the determination of the electronic bandstructure of solids have made considerable progress in recent years . From the theoretical point of view, state–of–the–art calculations within many–body perturbation theory allow to obtain band energies in a rigorous way, i.e. as the poles of the one–particle Green’s function $`G`$ . To obtain G, one needs the electron self-energy $`\mathrm{\Sigma }`$, usually evaluated according to the so-called GW approximation, derived by Hedin in 1965 , which is based on an expansion in terms of the dynamically screened Coulomb interaction. However, due to the high complexity and large computational requirements of ab–initio calculations of $`\mathrm{\Sigma }`$, the experimental bandstructures are often compared with the results of (simpler) calculations performed within Density Functional Theory (DFT), in the Local Density Approximation (LDA) or in the Generalized Gradient Approximation (GGA). The consequences of this approach (which can be set on a firm ground by considering the exchange–correlation potential of the LDA or GGA as an approximation to the self–energy operator) must however be considered with great care, particularly when the system under study differs from those –semiconductors and insulators– for which the approximation is usually made. A recent example is given by Ref. , where the measured Cu bandstructure is compared with state–of–the–art DFT-GGA results, finding discrepancies which are both large and dependent on the considered band and k–point. As a consequence of this result, the authors of Ref.  state that it is intriguing to find such pronounced “deviations from the band theory”, copper being much less correlated than metals with an open 3d shell (e.g., Ni).

In the present work, we demonstrate that copper does not, however, deviate from quasiparticle band theory. It only deviates from the simple “rigid shift” behavior of the self–energy commonly found in the case of semiconductors and insulators. For these, the application of the GW method to compute self–energy corrections ontop of ab–initio DFT results has become a quite well-established and standard technique, giving energy levels generally in good agreement with experiments, even for complicated systems like reconstructed surfaces and clusters. The gaps between empty and filled states generally increase by a substantial amount with respect to those obtained in the Kohn-Sham (KS) formulation of DFT, reaching agreement with experimental results. The KS–DFT underestimation of the experimental gap is generally found to be weakly dependent on the particular band or k–point, although some semiconductor surface gaps show stronger QP corrections than bulk ones . This is the basis for the introduction of the so–called “scissors operator”, often invoked in order to correct the discrepancies by rigidly shifting upwards the DFT empty bands, hence avoiding explicit self–energy calculations. A priori, there is no reason why the “scissors” approach should also work for metals, where no gap exists between filled and empty states.

Unlike semiconductors, the case of metals has received limited attention so far. The band width (i.e., the energy range of filled states) of simple metals has often been compared with that calculated for the homogeneous electron gas (jellium), finding discrepancies of a few tenths of eV. On the other hand, the band structure of transition and noble metals cannot be approximated by that of jellium: valence electrons of d character play an important role both in the electronic structure and in the optical properties of these metals. Previous estimates of many–body corrections to the 3d Cu bandstructure, based on the self–energy of the homogeneous electron gas, yielded only partial agreement with the experimental data, suggesting that in the case of non–free electron metals screening effects (accounted for by the dielectric function) should be included in a more realistic way. This is precisely the aim of the present work, where a full, first–principles GW calculation of the quasiparticle (QP) bandstructure of bulk copper is performed. The resulting QP energies correct in a non trivial way the DFT-LDA results, and are in remarkably good agreement with direct and inverse photoemission experiments.

Among transition metals, full quasiparticle calculations have been carried out so far only for Ni, while we are not aware of similar calculations for noble metals. In the case of Ni, $`GW`$ yields a good description of photoemission data, except for the $`6`$ eV satellite, which is due to strong short–range correlations within the partially–filled d–shell. For copper, we expect band–theory to work better than for transition metals, since d–shells are completely filled. From the computational point of view, the very presence of strongly localized d–bands makes the ab–initio calculations based on plane waves (PW) and norm–conserving pseudopotentials (NCPP) for noble metals much heavier than those for semiconductors.

Let us first recall the relations between Kohn–Sham eigenvalues, the poles of the one–particle Green’s function G, and experimental band energies measured in photoemission (PE) or inverse photoemission (IPE). In DFT, KS eigenvalues cannot be identified with electron addition or removal energies, since there is no equivalent of Koopman’s theorem. Instead, the experimental bandstructure should be compared with the poles of the Green’s Function $`G(𝐫,𝐫^{};\omega )`$. The latter are determined by an equation of the form:

$$\begin{array}{c}\left[\frac{\mathrm{}^2}{2m}\mathrm{}_𝐫+V_{external}\left(𝐫\right)+V_{Hartree}\left(𝐫\right)\right]\psi _{n𝐤}(𝐫,\omega )+\hfill \\ \hfill 𝑑𝐫^{}\mathrm{\Sigma }(𝐫,𝐫^{};\omega )\psi _{n𝐤}(𝐫^{},\omega )=E_{n𝐤}\left(\omega \right)\psi _{n𝐤}(𝐫,\omega ),\end{array}$$
(1)
containing the non-local, non-hermitian and frequency dependent self-energy operator $`\mathrm{\Sigma }`$. The poles of G are the QP energies $`ϵ_{n𝐤}^{QP}`$, the solutions of $`ϵ_{n𝐤}^{QP}=E_{n𝐤}\left(ϵ_{n𝐤}^{QP}\right)`$. Following Ref. , we separate the static, bare–exchange part $`\mathrm{\Sigma }_x(𝐫,𝐫^{})`$ from $`\mathrm{\Sigma }_c(𝐫,𝐫^{};\omega )`$, the energy–dependent correlation contribution. Eq. (1) is formally similar to Kohn-Sham equations which are solved in the determination of the ground-state properties, where the local and energy independent exchange-correlation potential $`V_{xc}\left(𝐫\right)`$ has been substituted by $`\mathrm{\Sigma }`$. Hence, KS eigenvalues can be considered as a zeroth-order approximation to the true QP energies, if the exchange-correlation (xc) potential of the DFT is seen as an approximation to the true self-energy operator $`\mathrm{\Sigma }`$. This suggests the possibility to look for a first-order, perturbative solution of Eq. (1) with respect to $`\left(\mathrm{\Sigma }V_{xcLDA}\right)`$, where $`V_{xcLDA}`$ is the xc potential of the LDA. We have checked that the off–diagonal elements of $`\left(\mathrm{\Sigma }V_{xcLDA}\right)`$ are two orders of magnitude smaller than the diagonal ones. To determinate the QP energies $`ϵ_{nk}^{QP}`$ starting from the DFT-LDA values we have used an iterative procedure.

Our bandstructure calculation hence starts with a DFT-LDA calculation of the ground-state properties, performed using norm-conserving pseudopotentials (PPs) and a plane waves basis. The use of soft (Martins-Troullier) PPs allows us to work at full convergence with a reasonable kinetic energy cutoff (60 Rydbergs if the 3s and 3p atomic states are frozen into the core, 160 Rydbergs when they are explicitly included into the valence). Details can be found in Ref. , where the non-trivial discrepancies between the LDA Kohn-Sham eigenvalues and the experimental bandstructure have also been evidentiated. In particular, we find substantial differences with respect to the experiment for both the d–bands width (3.70 instead of 3.17 eV) and position (more than 0.5 eV upshifted in the DFT), in agreement with previous calculations .

QP energies are then computed by evaluating $`\mathrm{\Sigma }`$ as GW , so that the dynamically screened potential $`W(\omega )`$ (convolution of the inverse dielectric function $`ϵ^1\left(\omega \right)`$ with the bare Coulomb potential) is needed. Most GW calculations on semiconductor systems use a Plasmon-Pole Approximation (PPA) for $`W(\omega )`$ , based on the observation that the Fourier components $`ϵ_{𝐆,𝐆^{}}^1(𝐪;\omega )`$ of the inverse dielectric function are generally peaked functions of $`\omega `$, and can be approximated by a single pole. Since the evaluation of $`\mathrm{\Sigma }_c`$ involves an integration over the energy, the fine details of the $`\omega `$-dependence are not critical, and the PPA turns out to work reasonably well for most applications. However, in the case of Cu the use of a PPA becomes more critical. The presence of flat d–bands 2 eV below the Fermi level implies the presence of strong transitions in $`ϵ_{𝐆,𝐆^{}}^1(𝐪;\omega )`$ spread over a large energy range. For small values of $`𝐆`$ and $`𝐆^{}`$, the behaviour of $`ϵ_{𝐆,𝐆^{}}^1(𝐪;\omega )`$ is often very different from that of a single-pole function, leading to instabilities when determining the Plasmon-Pole parameters. Hence, we have found it more convenient to avoid the PPA. Instead, we explicitly compute $`ϵ_{𝐆,𝐆^{}}^1(𝐪;\omega )`$ over a grid of about 200 frequencies from zero to $``$130 eV, and perform the energy integral numerically. Another characteristic of metallic systems which leads to additional difficulties in practical GW calculations is the discontinuity of the band occupation crossing the Fermi surface. The problem already shows up in the evaluation of the bare exchange matrix elements, which read:

$$\begin{array}{c}\mathrm{\Sigma }_x^{n𝐤}n𝐤|\mathrm{\Sigma }_x(𝐫,𝐫^{})|n𝐤=\hfill \\ \hfill \underset{n_1}{}_{BZ}\frac{d𝐪}{\left(2\pi \right)^3}\underset{𝐆}{}4\pi \\ \hfill \left|\frac{n𝐤|e^{i\left(𝐪+𝐆\right)𝐫}|n_1𝐤𝐪}{𝐪+𝐆}\right|^2f_{n_1\left(𝐤𝐪\right)}.\end{array}$$
(2)
Here $`0f_{n_𝐤}2`$ is the occupation of band $`n`$ at point $`𝐤`$, and its discontinuity at the Fermi Surface slows down the convergency of any numerical evaluation of Eq. (2) on a finite $`𝐪`$-mesh, no matter how $`n𝐤|e^{i\left(𝐪+𝐆\right)𝐫}|n_1𝐤𝐪`$ is smooth. To overcome this difficulty, we have first computed the integral of the term $`f_{n_1\left(𝐤𝐪\right)}/\left|𝐪+𝐆\right|^2`$ over small regular volumes centered around each $`𝐪`$-point used to evaluate Eq. (2), and then used these values in the finite $`𝐪`$-sum. A similar approach has been used in Ref. . This procedure yields results for $`\mathrm{\Sigma }_x^{n𝐤}V_{xcLDA}^{n𝐤}`$ computed with 19 $`𝐪`$-points in the IBZ converged within 0.1 eV, while barely evaluating Eq. (2) over the same set of $`𝐪`$-points yields an error of about 1.5 eV.

Normally (e.g., in GW calculations for semiconductors), the calculation of G and W to correct the DFT valence bandstructure can be performed by including only valence states, and fully neglecting the core states which have been frozen in the pseudopotential approach. In copper, however, when $`\mathrm{\Sigma }`$ is computed neglecting the 3s and 3p atomic core states (which in the solid create two flat bands, at about 112 and 70 eV, respectively, below the Fermi level), the resulting QP corrections on the d-bands are clearly unphysical: GW corrections move the highest occupied d–bands above the DFT-LDA Fermi level. On the other hand, the situation for s/p states (e.g. for the state $`L_2^{}`$) is much more reasonable, with correlation and exchange parts of the self-energy which cancel largely each other (as in the case of semiconductors), and negative QP corrections of order of the eV.

The solution of this puzzling situation is provided by the role of the above–mentioned 3s and 3p states, which, despite being well separated in energy from the 3d ones, have a large spatial overlap with the latter. As a consequence, non-negligible contributions to the self-energy are expected from the exchange contributions between 3d and 3s/3p states. These contributions can be included in our calculation in a straightforward way, by starting with a pseudopotential which puts the whole $`3^{rd}`$ atomic shell into the valence (see Ref.). Indicating with $`G_0^{core}`$ the non-interacting Green’s function containing also the 3s/3p levels, the self-energy computed as $`G_0^{core}W`$ will yield the desired dynamical exchange between 3d bands and core bands composed of a bare exchange (static) contribution, $`\mathrm{\Sigma }_x^{core}`$, plus a correlation part, $`\mathrm{\Sigma }_c^{core}`$. Naturally, a larger plane-waves cutoff (160 Ry) is needed for convergence. As a result, we find large bare–exchange contributions from core levels, leading to very different values of $`\mathrm{\Sigma }_x^{n𝐤}V_{xcLDA}^{n𝐤}`$ for s/p and d states, as illustrated in Fig. 1. The role of core levels in the calculation of the bare exchange contributions (whose importance was already addressed for transition metals by Aryasetiawan and Gunnarsson , but estimated to be of the order of 1 eV) is hence crucial in the case of copper . Moreover their effect is unexpected on the basis of DFT–LDA calculations, where the band structure does not change appreciably even if the 3s/3p orbitals are fully included in the valence . The effect of the core states is expected to be smaller on the correlation part of $`\mathrm{\Sigma }\left(\omega \right)`$, since the screened and bare Coulomb interaction, whose difference is involved in $`\mathrm{\Sigma }_c^{core}`$, are very similar to each other in the core region. As a result, $`\mathrm{\Sigma }_c^{core}`$ is small, of the order of $`0.5`$ eV for d–bands, and close to zero for s/p–bands. The inclusion of the correlation part $`\mathrm{\Sigma }_c^{core}`$ does not change the widths of d-bands, that remain in good agreement with the experimental results, but improves considerably their absolute positions.

In Tab. I we compare our calculations with selected experimental data reported in Ref. . In Fig. 2 the full theoretical band structure is compared with experimental data: the agreement is remarkably good and the fact that the GW corrections cannot be reproduced by any rigid shift of the LDA bands clearly appears. For instance, at the L point the shifts change sign for different valence bands, with QP corrections ranging from $`0.29`$ eV to $`0.61`$ eV, to $`0.57`$ eV for the three d–bands. As a conseguence, the $`L_1L_2^{}`$ gap decreases after inclusion of self–energy corrections (see also Tab. I, last line), in contrast with the usual behaviour occurring in semiconductors and insulators. The other gaps at $`L`$ and $`X`$, instead, increase. All of them come to good agreement with experimental data. Unoccupied (conduction) bands at points L, X and K are also obtained in very good agreement with the experimental data. These findings demonstrate that the strong deviations from the single quasiparticle band theory in copper, as suggested in Ref. , do not occur. QP band–theory in the $`GW`$ approximation, instead, remarkably deviates from the DFT–LDA + scissors–operator approach, due to the interplay of the different localization and correlation properties of d and s/p states.

In conclusion, we have carried out the first ab–initio quasi-particle calculation for a noble metal, copper. In order to do that, we had to overcome several problems, namely, the large number of plane waves needed to describe d orbitals, the calculation of k-space sums close to the Fermi surface, and the inclusion of the exchange interaction between valence and core levels. The single–quasiparticle bandstructure computed within Hedin’s $`GW`$ approximation for the electron self–energy turns out to be in very good agreement with experiments. The GW method, originally devised to describe the long-range charge oscillations , is hence shown to yield a good description also of copper, a system characterized by localized orbitals and short-range correlation effects. The corrections to the LDA Kohn–Sham eigenvalues for the d-bands are highly non–trivial, since even their sign turns out to be k-point and band–dependent, ruling out the validity of any “scissors operator” approximation. Furthermore, exchange effects between d–electrons and 3s, 3p core states play a key role and cannot be neglected in the calculation of QP corrections. Finally, a full screening calculation of the self-energy has been carried out, avoiding the single-plasmon pole approximation. This allows a realistic treatment also of the imaginary part of the self energy (and therefore the calculation of QP lifetimes), which will be studied in a forthcoming work.

This work has been supported by the INFM PRA project “1MESS”, MURST-COFIN 99 and by the EU through the NANOPHASE Research Training Network (Contract No. HPRN-CT-2000-00167). We thank Angel Rubio, Lucia Reining and Michele Cini for useful discussions, and Conor Hogan for a critical reading.


text:
# HIGH ENERGY CONSTRAINTS ON LORENTZ SYMMETRY VIOLATIONS

## 1 Introduction

There are several reasons to suspect that Lorentz invariance may be only a low energy symmetry. This possibility is suggested by the divergences of local quantum field theory, as well as by tentative results in various approaches to quantum gravity. Moreover, the unboundedness of the boost parameter makes experimental verification of exact Lorentz symmetry impossible in principle.

One can study the possibility of Lorentz violation, without a particular fundamental theory in hand, by considering its manifestation in dispersion relations for matter. It is natural to assume that such dispersion relations $`E^2(p)`$ can be characterized by an expansion with integral powers of momentum,

$$E^2=p^2+m^2+\left[Ap^2+Bp^3/K_0+Cp^4/K_0^2+O(p^5)\right].$$
(1)
Here $`A,B,C`$ are dimensionless coefficients which might be positive as well as negative and $`K_0`$ is the “quantum gravity” scale (often identified with the inverse Planck length, $`K_0=2\pi /L_{\mathrm{Pl}}`$). \[Throughout this paper $`p`$ denotes the absolute value of the 3-momentum vector $`𝐩`$, and we use units with the low energy speed of light in vacuum equal to unity.\]

Different approaches to quantum gravity suggest different leading order Lorentz violating terms. The terms with coefficients $`A,B,C`$ have mostly been considered so far. Since the $`p^2`$ term is not suppressed by the Planck scale, it might be thought to be largest, however observations severely limit $`A`$ to be much less than unity (see e.g.  and references therein). The higher order terms are naturally small, with coefficients of order unity.

Our strategy here is to take a purely phenomenological stance and consider the constraints that high energy observations impose on dispersion relations of the form

$$E_a^2p_a^2+m_a^2+\eta _ap_a^n/K_0^{n2},$$
(2)
where $`a`$ labels different fields and $`n3`$. In the absence of a fundamental theory one has no reason to expect any particular relation between the coefficients $`\eta _a`$ for different particles, except perhaps that they should all be of the same order of magnitude.

Dispersion relations of the kind (2) produce kinematic constraints from energy-momentum conservation that differ from the usual Lorentz invariant case. As a result reactions can take place that are normally forbidden, and thresholds for reactions are modified. Observational consequences may seem out of reach because of the Planck scale suppression of the Lorentz violating terms (assuming that, as generally expected, $`K_0`$ is of order the Planck scale). However this is not so. One can expect deviations from standard kinematics when the the last two terms of (2) are of comparable magnitude.111In general of course one must look at the specific reaction in order to estimate the energies at which deviations from standard behavior can be expected. For example in the case of the reaction $`\gamma \gamma e^+e^{}`$ it will be the electron mass that sets the scale. Assuming $`\eta `$ is of order unity this yields the condition $`p_{\mathrm{dev}}(m^2K_0^{n2})^{1/n}`$, which is $`(m/m_e)^{2/3}\times 10`$ TeV for $`n=3`$ and $`(m/m_e)^{1/2}\times 10^4`$ TeV for $`n=4`$. Although these energies are currently not achievable in particle accelerators (except in the case of the massive neutrinos which however are weakly coupled) they are in the range of current astrophysical observations.

In fact, it has been suggested by several authors  (see also  and references therein) that we may already be observing deviations from Lorentz invariance via the existence of two puzzles in modern astrophysics: the missing GZK cut-off on cosmic ray protons with ultra high energy greater than $`7\times 10^{19}`$ eV, and the apparent overabundance of gamma rays above 10 TeV from the BL Lac system Mkr 501 . Here we shall mostly not consider the constraints imposed by asking Lorentz violation to explain these puzzles. Instead we restrict our attention to constraints imposed by consistency with known phenomena (or lack thereof). (See also  for a similar discussion.)

## 2 Observational constraints

To find the strongest observational constraints without assuming a priori relations between the coefficients $`\eta _a`$ of (2) we focus here on processes involving just two fundamental particles, photons and electrons. We also restrict to the case $`n=3`$, since this should be most tightly constrained. If it can be ruled out, one can then move on to the $`n=4`$ case.

The modified dispersion relations for photons and electrons in general allow two normally forbidden interactions: photon decay, $`\gamma e^+e^{}`$, and vacuum Čerenkov radiation, $`e^{}e^{}\gamma `$. If allowed these processes happen very rapidly. In addition the threshold for photon annihilation, $`\gamma \gamma e^+e^{}`$, is shifted. We consider constraints that follow from three observations: (i) Gamma rays up to $`50`$ TeV of cosmological origin arrive on earth , so photon decay does not occur up to this energy. (ii) Electrons of energy $`100`$ TeV are believed to produce observed X-ray synchrotron radiation coming from supernova remnants. Assuming these electrons are actually present, vacuum Čerenkov radiation must not occur up to that energy. (iii) Cosmic gamma rays below 10 TeV are believed to be absorbed in a manner consistent with photon annihilation off the IR background with the standard threshold. Observation (iii) is not model independent, so the corresponding constraint is tentative and subject to future verification.

To derive the observational constraints one needs to determine the threshold for each the process, i.e. the lowest energy for which the process occurs. The details concerning our determination of the thresholds are reported elsewhere . Assuming monotonicity of $`E(p)`$ (for the relevant momenta $`pK_0`$) we have shown that all thresholds for processes with two particle final states occur when the final momenta are parallel. Moreover for two particle initial states the incoming momenta are antiparallel. This geometry has been assumed in previous works but to our knowledge it was not shown to be necessary. (In fact it is not necessary if $`E(p)`$ is not monotonic.)

To eliminate the subscript $`a`$ we introduce $`\xi :=\eta _\gamma `$ and $`\eta :=\eta _e`$. The constraints will restrict the allowed region of the $`\eta `$$`\xi `$ plane. For the rest of the paper we assume $`K_0`$ is the Planck energy and we use units with $`K_0=1`$.

### 2.1 Photon decay

Photon decay is allowed only above a broken line in the $`\eta `$$`\xi `$ plane given by $`\xi =\eta /2`$ in the quadrant $`\xi ,\eta >0`$ and by $`\xi =\eta `$ in the quadrant $`\xi ,\eta <0`$. Above this line, the threshold is given by

$`k_{\mathrm{th}}`$ $`=`$ $`\left({\displaystyle \frac{8m^2}{2\xi \eta }}\right)^{1/3}\text{for }\xi 0,`$ (3)
$`k_{\mathrm{th}}`$ $`=`$ $`\left({\displaystyle \frac{8m^2\eta }{(\xi \eta )^2}}\right)^{1/3}\text{for }\eta <\xi <0\text{.}`$ (4)
The first relation (3) arises when the electron and positron momenta are equal at threshold. In standard Lorentz invariant kinematics such “equipartition” of momentum always holds at threshold for pair production. In all previous work on Lorentz violating dispersion it has been assumed to hold. Surprisingly, however, in the present case the threshold may occur with an asymmetric distribution of momentum. The second relation (4) applies in those cases. The constraint we impose is that the threshold is above 50 TeV, the energy of highest observed gamma rays from the Crab nebula .

### 2.2 Photon annihilation

The threshold relations for a gamma ray to annihilate with an IR background photon of energy $`ϵ`$ take approximately the same form as for photon decay, with the replacement $`\xi \xi ^{^{}}`$, where $`\xi ^{^{}}\xi +4ϵ/k_{\mathrm{th}}^2`$. (Here we have used the fact that $`ϵ`$ is much smaller than any other scale in the problem.) The two different relations arise for the same reason as in the case of photon decay, and they correspond respectively to cubic and quartic polynomial equations for $`k_{\mathrm{th}}`$.

For the observational consequences it is important to recognize that the threshold shifts are much more significant at higher energies than at lower energies. To exhibit this dependence, it is simplest to fix a gamma ray energy $`k`$ and to solve for the corresponding soft photon threshold energy $`ϵ_{\mathrm{th}}`$. Taking the ratio with the usual threshold $`ϵ_{\mathrm{th},0}`$, we find a dependence on $`k`$ at least as strong as $`k^{3/2}`$. Introducing $`k_{10}:=k/(10\text{TeV})`$, we have

$`{\displaystyle \frac{ϵ_{\mathrm{th}}}{ϵ_{\mathrm{th},0}}}`$ $`=`$ $`1+0.05(\eta 2\xi )k_{10}^3\text{for }\xi ^{}0,`$ (5)
$`{\displaystyle \frac{ϵ_{\mathrm{th}}}{ϵ_{\mathrm{th},0}}}`$ $`=`$ $`0.1(\eta \xi )k_{10}^3+\sqrt{0.2\eta k_{10}^3}\text{for }\eta <\xi ^{}<0\text{.}`$ (6)
High energy TeV gamma rays from the blazars Markarian 421 and Markarian 501 have been detected out to 17 TeV and 24 TeV respectively. Although the sources are not well understood, and the intergalactic IR background is also not fully known, detailed modeling shows that the data are consistent with some absorption by photon annihilation off the IR background (see e.g.  and references therein). However, while the inferred source spectrum for Mrk 501 is consistent with expectations for energies less than around 10 TeV, above this energy there are far more photons than expected according to some IR background models .

The uncertainty in the blazar and IR background models currently precludes sharp constraints from photon annihilation. Instead, we just determine the range of parameters $`\xi ,\eta `$ for which the threshold $`k_{\mathrm{th}}`$ lies between 10 TeV and 20 TeV for an IR photon with which a 10 TeV photon would normally be at threshold. We choose this range since (i) lowering the threshold would make the overabundance problem worse, and (ii) raising this threshold by a factor of two might explain the potential overabundance of photons over 10 TeV . If the overabundance puzzle exists and is to be resolved in this way, the effect of the threshold shifts must be enhanced for photons above 10 TeV relative to those below 10 TeV (since the latter seem to be well accounted for). This enhancement could arise partly from the shape of the IR background spectrum, but the shift itself is also enhanced, as seen in the energy dependence of equations (5,6). The hypothesis  that the potential overabundance is due to Lorentz violating dispersion with $`n=3`$ therefore appears to be consistent with current observations.

### 2.3 Vacuum Čerenkov radiation

An electron can emit Čerenkov radiation in the vacuum if $`\eta >0`$ or if $`\eta <0`$ and $`\xi <\eta `$. In this case there are also two threshold relations, depending on whether the threshold occurs with emission of a zero-energy photon or with emission of a finite energy photon. These two cases correspond to the two following relations, respectively:

$`p_{\mathrm{th}}`$ $`=`$ $`\left({\displaystyle \frac{m^2}{2\eta }}\right)^{1/3}\text{for }\eta >0\text{ and }\xi 3\eta ,`$ (7)
$`p_{\mathrm{th}}`$ $`=`$ $`\left({\displaystyle \frac{4m^2\left(\xi +\eta \right)}{\left(\xi \eta \right)^2}}\right)^{1/3}\text{for }\xi <3\eta <0\text{ or }\xi <\eta 0\text{.}`$ (8)
The vacuum Čerenkov process is extremely efficient, leading to an energy loss rate that goes like $`E^2`$ well above threshold. Thus any electron known to propagate must lie below the threshold. Electrons of energy $`100`$ TeV are believed to produce observed X-ray synchrotron radiation coming from supernova remnants. Thus for example in the region of the parameter plane where (7) holds we obtain the constraint $`\eta <m^2/2p_{\mathrm{th}}^310^3`$.

## 3 Combined constraints

Putting together all the constraints and potential constraints we obtain the allowed region in the $`\eta `$$`\xi `$ plane (see Figure 1). The photon decay and Čerenkov constraints exclude the horizontally and vertically shaded regions, respectively. The allowed region lies in the lower left quadrant, except for an exceedingly small sliver near the origin with $`0<\eta \stackrel{<}{}\mathrm{\hspace{0.33em}10}^3`$ and a small triangular region ($`0.16\stackrel{<}{}\eta <0`$, $`0<\xi \stackrel{<}{}\mathrm{\hspace{0.33em}0.08}`$) in the upper left quadrant. The range of the photon annihilation threshold discussed in subsection 2.2 falls between the two roughly parallel diagonal lines. This intersects the otherwise allowed region in a finite, narrow wedge where $`\xi `$ and $`\eta `$ are negative and of order unity (apart from a minuscule invisible region near the origin with $`\eta >0`$).

If future observations of the blazar fluxes and the IR background confirm agreement with standard Lorentz invariant kinematics, the region allowed by the photon annihilation constraint will be squeezed toward the upper line ($`k_{\mathrm{th}}k_\mathrm{s}`$). If the overabundance of gamma rays from Mrk 501 is indeed due to Lorentz violation of the sort we are considering, the region will be squeezed toward the lower line ($`k_{\mathrm{th}}2k_\mathrm{s}`$) or may even shift toward lower values of $`\xi `$ and $`\eta `$ if a yet larger threshold shift is indicated.

## 4 Conclusions

We have seen that a conservative interpretation of observations puts strong constraints on the possibility of Planck scale cubic modifications to the electron and photon dispersion relations. The allowed region includes $`\xi =\eta =1`$, which has been a focus of previous work . The negative quadrant has most of the allowed parameter range. It is interesting to note that in this quadrant all group velocities are less than the low energy speed of light.

To completely rule out the cubic case would require new observations. Finding higher energy electrons would not help much, while finding higher energy undecayed photons would squeeze the allowed region onto the line $`\xi =\eta `$. To further shrink the allowed segment of this line would require observations allowing the usual threshold for photon annihilation to be confirmed to higher precision. Perhaps other processes could be used as well. If a priori relations among the coefficients in the dispersion relations for different particles are hypothesized, stronger constraints can of course be obtained.

This research was supported in part by NSF grant PHY-9800967.

## References


text:
# Ultrahigh energy neutrinos, small 𝒙 and unitarity

## I Introduction

Ultrahigh energy neutrinos are predicted from a number of sources. One source is from cosmic ray interactions with the microwave background radiation Protheroe:1996ft , producing charged pions which decay into neutrinos. Another possible source is decaying cosmic strings or extremely massive relics Bhattacharjee:2000qc , which ultimately contribute to a cosmic neutrino flux. Detection of ultrahigh energy neutrinos may shed light on the observation of air shower events with energies in excess of $`10^{11}`$ GeV, reveal aspects of grand unification or yield some insight into the sources of the highest energy cosmic rays.

A number of detectors are able or will be able to detect neutrino induced showers. For example, the Auger experiment should be able to detect neutrino induced horizontal air showers initiated by neutrinos with energies above $`10^9`$ GeV. The proposed OWL/EUSO satellite experiments should be able to detect upward air showers produced by $`\nu _\tau \tau `$ just below the Earth’s surface. The event rates predictions depend on the ultrahigh energy neutrino cross section, extrapolated beyond the measured regime, as well as on the predicted neutrino fluxes.

Recent discussions by Dicus, Kretzer, Repko, and Schmidt Dicus:2001kb  about the implications of perturbative unitarity have refocused attention on the ultrahigh energy extrapolation of the neutrino-nucleon cross section. In the next section we review the cross section evaluation including the extrapolation of the parton distribution function to small parton momentum fraction $`x`$. We examine to what extent the cross section may be sensitive to the presence of saturation effects in the evolution of the parton distributions. In the following section, we outline the unitarity argument and comment on what can and cannot be learned by relating the neutrino-nucleon forward scattering amplitude to the total neutrino-nucleon cross section.

## II Neutrino-nucleon cross section

The expression for the neutrino-nucleon charged current cross section, for $`\nu _l(k)N(p)l(k^{})X(p^{})`$, is

$$\frac{d^2\sigma }{dxdQ^2}=\frac{G_F^2}{\pi }\left(\frac{M_W^2}{Q^2+M_W^2}\right)^2\left[q(x,Q)+(1y)^2\overline{q}(x,Q)\right],$$
(1)
in terms of $`Q^2=(kk^{})^2`$, $`x=Q^2/(2M\nu )`$ and $`y=\nu /E`$, for neutrino energy $`E`$ and lepton energy transfer $`\nu =EE^{}`$ defined in the nucleon rest frame. The nucleon mass is $`M`$, and the center of mass energy squared is $`S=2ME`$. In Eq. (1) we have followed Refs. Dicus:2001kb ; Gandhi:1998ri  to introduce effective quark and antiquark densities that contain the contributions from the various flavors as well as the appropriate electroweak mixing angles. The expression for the neutral current reaction $`\nu _l(k)N(p)\nu _l(k^{})X(p^{})`$ can be cast into an identical form, with the obvious replacement $`M_WM_Z`$ and with different effective quark densities. In what follows, we will only consider the neutrino-nucleon cross section. At high energies, the antineutrino cross section is expected to be very similar. We will neglect perturbative QCD corrections to the cross section which were found to be small in Refs. Dicus:2001kb ; Gluck:1999js . Finally, we also neglect the contributions to the cross section arising from charm quarks in the initial state; these can be sizable at high energies, but are unimportant for our more qualitative purposes.

In Eq. (1), increasing $`Q^2`$ has two effects: as $`Q^2`$ rises, the cross section decreases due to the $`W`$ propagator, but the contributions of the quark and antiquark distribution functions $`q(x,Q)`$ and $`\overline{q}(x,Q)`$ increase due to QCD evolution Andreev:1979cp . The propagator dominates and effectively cuts off the growth in $`Q^2`$ at $`Q^2M_W^2`$. As a consequence, the typical $`x`$ value as a function of incident neutrino energy $`E`$ is $`xM_W^2/(2ME)`$, so ultrahigh neutrino energies translate to small parton $`x`$. For the highest neutrino energies considered, $`E10^{12}`$ GeV, $`x10^8`$. HERA measurements of structure functions hera  extend to $`x10^6`$; however, such low values of $`x`$ are measured at $`Q^2<0.1`$ GeV2. For $`QM_W`$, the structure functions are measured down to $`x10^3`$ in the D0 experiment’s analysis of inclusive jets Babukhadia:2001bd . Small $`x`$ extrapolations of the parton distribution functions are therefore necessary to extend the predictions for the neutrino-nucleon cross section above $`E10^7`$ GeV.

Perturbative QCD governs the small $`x`$ extrapolations. The sea quark distributions dominate the cross section at high energies. Sea quarks are produced by gluon splitting $`gq\overline{q}`$, so the gluon distribution $`g(x,Q)`$ dictates the eventual quark and antiquark distributions at small $`x`$. The gluon distribution is parameterized as $`xg(x,Q_0)x^\lambda `$ for $`x1`$ at a reference scale $`Q_0`$. Approximate small-$`x`$ DGLAP evolution dglap , for $`\lambda `$ close to 0.5, yields a gluon distribution function of the same form, at a larger value of $`Q`$: $`xg(x,Q)x^\lambda `$ Ellis:1994rb . As a practical matter, $`\lambda `$ was determined at $`Q=M_W`$ to extrapolate the parton distribution functions Gandhi:1998ri , for example, those by the CTEQ collaboration Lai:1997mg , below $`x=10^5`$. Glück, Kretzer, and Reya Gluck:1999js  have checked that the full DGLAP evolution of the Glück, Reya, and Vogt Gluck:1998xa  distribution functions yields only some 20% difference at $`x=10^8`$ compared with the power law extrapolation of the CTEQ densities. Kwiecinski, Martin, and Stasto Kwiecinski:1999yf  have performed a BFKL-type bfkl  evolution, yielding results in substantial agreement at the highest energies considered (1012 GeV). The resulting total neutrino-nucleon cross sections can be parameterized by power laws for $`10^7`$ GeV$`<E<10^{12}`$ GeV. For example, the charged current neutrino-nucleon cross section, using the CTEQ4 parton distribution functions, scales as $`\sigma =5.5\times 10^{36}(E/\mathrm{GeV})^{0.36}`$ cm2 Gandhi:1998ri .

Ultimately, the growth of the parton distribution functions – and hence that of the cross section – predicted by both DGLAP and BFKL evolution will have to slow down when the gluon densities become large enough that $`ggg`$ recombination processes become important sat . The regime in $`x,Q^2`$ in which this happens is referred to as “saturation region” and can be estimated. It can be roughly characterized by the condition gbw

$$Q_s^2(x)=(1\mathrm{G}\mathrm{e}\mathrm{V}^2)\left(\frac{x_0}{x}\right)^\lambda $$
(2)
with $`\lambda `$ and $`x_0`$ obtained from a fit to HERA data: $`\lambda =0.288`$, $`x_0=3.04\times 10^4`$. $`Q_s`$ is the “saturation scale” – saturation should roughly occur when $`Q<Q_s(x)`$ at a given $`x`$. Clearly, as $`x`$ decreases, saturation effects are expected to become relevant already at larger and larger $`Q^2`$.

In Figs. 1 and 2 we examine whether the neutrino-nucleon cross section at ultrahigh energy, $`E=10^{12}`$ GeV, might be sensitive to such saturation effects. Figure 1 shows $`\mathrm{log}_{10}(d^2\sigma /d\mathrm{log}_{10}(x)d\mathrm{log}_{10}(Q^2))`$, evaluated using the parton distribution functions of Gluck:1998xa , as a function of $`\mathrm{log}_{10}x`$ and $`\mathrm{log}_{10}Q^2`$. It is evident that scales $`QM_{W,Z}`$ dominate the cross section. Note that the parton distributions of Gluck:1998xa  can be used also down to rather small $`Q^21`$ GeV2. This is convenient because, in order to obtain the total neutrino-nucleon cross section one needs to integrate the expression in Eq. (1) over the range $`0<Q^2<xS`$. For $`Q<1`$ GeV, we have frozen the scale in the parton distribution functions to 1 GeV. Fortunately, as can be seen from Fig. 1, the region of $`Q^2<1`$ GeV2 contributes only very little to the total cross section.

The line in the $`x,Q^2`$ plane corresponding to the saturation condition, Eq. (2), is also shown in Fig. 1. In Fig. 2 we follow Ref. Kwiecinski:1999yf  to look at a projection onto the $`x,Q^2`$ plane. The contours are the lines at which the cross section $`d^2\sigma /dxdQ^2`$ has fallen to $`10^{310.5n}`$ cm2, where $`n=1,\mathrm{},6`$. From these figures, it becomes evident that, even at the highest neutrino energies, contributions to the cross section resulting from the regime sensitive to gluon recombination effects are marginal. For our example at $`E=10^{12}`$ GeV, the “saturation region” contributes far less than 1% to the total cross section.

## III Unitarity considerations

Dicus, Kretzer, Repko, and Schmidt in Ref. Dicus:2001kb  have brought unitarity considerations to the fore. A restatement of the optical theorem relates the total neutrino-nucleon cross section to the neutrino-nucleon forward elastic scattering amplitude. The latter can be written in terms of the differential elastic cross section, evaluated at Mandelstam variable $`t=0`$, which, with some approximations, yields:

$$\frac{d\sigma _{\mathrm{el}}}{dt}|_{t=0}\frac{1}{16\pi }\sigma _{\mathrm{tot}}^2.$$
(3)
One can view this as a lower bound on the forward scattering elastic cross section, or as an upper bound on the total cross section. In Ref. Dicus:2001kb , the authors observe that the inequality is saturated at a relatively low energy by using the lowest order, $`G_F^2`$, contribution for the elastic cross section on the left, and the $`G_F^4`$ contribution that comes from the inclusive cross section (1) on the right. Specifically, using the leading term in $`G_F`$ for the elastic differential cross section, they conclude that

$$\sigma _{\mathrm{tot}}9.3\times 10^{33}\mathrm{cm}^2,$$
(4)
which already is violated for $`E2\times 10^8\mathrm{GeV}`$. From this they deduce that at yet higher energies, where the right-hand side of Eq. (3) increases, while the left is constant (at $`𝒪(G_F^2)`$), previously neglected terms that are higher order in the weak coupling $`g`$, in particular, $`g^6`$ or $`g^8`$ terms, must become important. They go on to suggest that this signals a breakdown in perturbation theory in the weak coupling, $`g`$. This is a striking implication indeed, especially given the small size of the cross section in Eq. (4).

There is another, and we believe more natural, interpretation of the equality when $`E2\times 10^8\mathrm{GeV}`$. First, we observe that the forward elastic cross section receives two qualitatively different and quantum mechanically incoherent contributions. The first of these describes the coherent elastic scattering of the entire nucleon through weak vector boson exchange, which begins at tree level, that is, at $`G_F^2`$ in the cross section. The second is the contribution of high-$`Q^2`$ virtual states that results from the incoherent scattering of partons. The latter, not the former, is related independently by the optical theorem to the inelastic cross section on the right-hand side of Eq. (3), and will saturate that inequality identically at order $`G_F^4`$, regardless of its size, just as at order $`G_F^2`$ the forward cross section is identically equal to the corresponding contribution from the square of the real part, which has been neglected on the right of Eq. (3).

That being said, we may still ask whether the dominance of the partonic part of the cross section, higher-order by $`g^2`$ compared to the elastic part, might not be a sign of large contributions from yet higher orders in the weak coupling. Integrating the factorized form Eq. (1), over $`x`$ and $`Q^2`$, however, shows that at very high energy the square of the total cross section behaves as $`G_F^2`$ times $`[g^2(S/M_W^2)^\lambda ]^2`$. This is to be compared to $`G_F^2`$ on the left-hand side of Eq. (3). The factor $`g^2`$ is the default size of a higher-order electroweak correction. The factor $`(S/M_W^2)^\lambda `$ is due to the large number of partons of size $`1/M_W`$ at $`xM_W^2/S`$. For higher orders in $`g^2`$ to contribute at a similar level, they would have to come accompanied by a similar large counting factor. At the leading power in $`1/M_W`$, which is given by Eq. (1), this cannot happen, simply because $`q(x)`$ and $`\overline{q}(x)`$ already count the partons. It would still be possible if more partons are involved in the hard scattering, but this involves going to higher twist, that is, to explicit suppression by additional powers of $`1/M_W`$, which would have to be compensated for by higher-twist multi-parton matrix elements. While such contributions are not very well-known even at low momentum transfers, there is no experimental indication of such large scales implicit within the nucleon.

The forgoing arguments, of course, assume that the unaided QCD extrapolations described above are equal to the task of so many orders of magnitude. We have shown above the self-consistency of these extrapolations, and that they do not, by themselves, lead to problems with unitarity, or give evidence of a breakdown in perturbation theory in the weak coupling sigl . The very fact of the self-consistency of the QCD extrapolations shows that ultra high energy neutrinos offer an exploration of the strong interactions, as well as of cosmic dynamics, into unprecedented length scales.

###### Acknowledgements.

We thank Duane Dicus, Stefan Kretzer, and Jamal Jalilian-Marian for discussions. This work was supported in part by NSF Grants No. PHY-9802403 (M.H.R.), PHY-0098527 (G.S.), DOE Contracts DE-FG02-95ER40906 and DE-FG03-93ER40792 (I.S.). M.S. thanks SUNY Stony Brook and RIKEN, and G.S. and M.S. thank Brookhaven National Laboratory for hospitality and support. W.V. is grateful to RIKEN, Brookhaven National Laboratory and the U.S. Department of Energy (contract number DE-AC02-98CH10886) for providing the facilities essential for the completion of his work.