text:
# Optimization and Parallelization of a force field for silicon using OpenMP

## I Program survey

Title of program: siliconiap

Computer hardware and operating system: Any shared memory computer running under Unix or Linux

Programming Language: Fortran90 with OpenMP compiler directives

Memory requirements: roughly 150 words per atom

No. of bits in a word: 64

No. of processors used: tested on up to 4 processors

Has the code been vectorized or parallelized: Parallelized withe OpenMP

No. of bytes in distributed program, including test data, etc: 50 000

Distribution format: Compressed tar file

Keywords: Silicon, Interatomic potential, Force field, Molecular dynamics

Nature of physical problem: condensed matter physics

Method of solution: Interatomic potential

Restrictions on the complexity of the problem: None

Typical running time: 30 $`\mu `$sec per step and per atom on a Compaq DEC Alpha

Unusual features of the program: None

## II Introduction

Due to its technological importance, silicon is one of the most studied materials. For small system sizes ab-initio density functional calculations  are the preferred approach. Unfortunately this kind of calculation becomes unfeasible for larger systems required to study problems such as interfaces or extended defects. For this type of calculations one resorts to force fields which are several orders of magnitude faster. Recent progress in the development of force fields has demonstrated that they can be a reliable tool for such studies. A highly accurate silicon force field has been developed by Lenosky and coworkers . Its transferability has been demonstrated by extensive tests containing both bulk and cluster systems . Its accuracy is in part due to the fact that second nearest neighbor interactions are included. This makes it unfortunately somewhat slower than force fields containing only nearest neighbor interactions. In the following a highly optimized parallel implementation of this force field will be presented that allows large scale calculations with this force field. The parallelization is achieved by using OpenMP, an emerging industry standard for medium size shared memory parallel computers.

Molecular dynamics calculations  have also been parallelized on distributed memory supercomputers . This approach is considerably more complex than the one presented here. Since few researches have access to massively parallel supercomputers and are willing to overcome the complexities of doing molecular dynamics on such machines, medium scale parallelization  of molecular dynamics has an important place in practice.

## III Calling the subroutine

User friendliness was one of the major design goals in the development of this routine. Using Fortran90 made it possible to hide all the complexities in an object oriented fashion from the user. The calling sequence is just

```
call lenosky(nat,alat,rxyz,fxyz,ener,coord,ener_var,coord_var,count)

```

On input the user has to specify the number of atoms, $`nat`$, the vector $`alat`$ containing the 3 lattice constant of the orthorhombic periodic volume and the atomic positions $`rxyz`$. The program then returns the total energy, $`ener`$, the forces, $`fxyz`$, the average coordination number, the variation of the energy per atom and of the coordination number as well as an counter that is increased in each call. In particular the user has not to supply any Verlet list.

Since the calculation of the forces is typically much more expensive than the update of the atomic positions in molecular dynamics or geometry optimizations, we expect that the subroutine will be called in most cases from within a serial program. In case the user is on a shared memory machine the subroutine will then nevertheless be executed in parallel if the program is compiled with the appropriate OpenMP options.

In addition the subroutine can of course also be used on a serial machine. In this case all the parallelization directives are considered by the compiler to be comments.

## IV Calculation of the Verlet list

The Verlet list gives all the atoms that are contained within the potential cutoff distance $`cut`$ of any given atom. Typically the Verlet list consists of two integer arrays. The first array, called $`lsta`$ in this work, points to the first/last neighbor position in the second array $`lstb`$ that contains the numbering of the atoms that are neighbors. A straightforward implementation for a non-periodic system containing $`nat`$ atoms is shown below. In this simple case the search through all atoms is sequential with respect to their numbering and it is redundant to give both the starting positions $`lsta(1,iat)`$ and the ending position $`lsta(2,iat)`$, since $`lsta(1,iat)=lsta(2,iat1)+1`$. But in the more complicated linear scaling algorithm to be presented below, both will be needed.

```
        indc=0
      do 10 iat=1,nat
c       starting position
        lsta(1,iat)=indc+1
        do 20 jat=1,nat
        if (jat.ne.iat) then
          xrel1= rxyz(1,jat)-rxyz(1,iat)
          xrel2= rxyz(2,jat)-rxyz(2,iat)
          xrel3= rxyz(3,jat)-rxyz(3,iat)
          rr2=xrel1**2 + xrel2**2 + xrel3**2
           if ( rr2 .le. cut**2 ) then
           indc=indc+1
c          nearest neighbor numbers
           lstb(indc)=jat
           endif
        endif
20      continue
c       ending position
        lsta(2,iat)=indc
10      continue

```

This straightforward implementations has a quadratic scaling with respect to the numbers of atoms. Due to this scaling the calculation of the Verlet list starts to dominate the linear scaling calculation of the energies and forces for system sizes of more than 10 000 atoms. It is therefore good practice to calculate the Verlet list with a modified algorithm that has linear scaling  as well.

To do this one first subdivides the system into boxes that have a side length that is equal to or larger than $`cut`$ and then finds all the atoms that are contained in each box. The CPU time for this first step is less than 1 percent of the entire Verlet list calculation. Hence this part was not parallelized. It could significantly affect the parallel performance according to Amdahls law  only if more than 50 processors are used. The largest SMP machines at our disposal had however only 4 processors.

To implement periodic boundary conditions all atoms within a distance $`cut`$ of the boundary of the periodic volume are replicated on the opposite part as shown in Figure 1. This part is equally well less than 1 percent of the CPU time for the Verlet list for a 8000 atom system. Being a surface term it becomes even smaller for larger systems. Consequently it wasn’t parallelized either.

After these two preparing steps one has to search only among all the atoms in the reference cell containing the atom for which one wants to find its neighbors as well as all the atoms in the cells neighboring this reference cell (26 cells in 3 dimensions). This implies that starting and ending positions $`lsta`$ for the atoms $`1`$ to $`nat`$ are not calculated in a sequential way, necessitating, as mentioned before, separate starting and ending positions in the array $`lsta`$. The corresponding parallel code is shown below. The indices $`l1,l2,l3`$ refer to the cells, $`icell(l1,l2,l3,0)`$ contains the number of atoms in cell $`l1,l2,l3`$ and $`icell(l1,l2,l3,)`$ their numbering. The array $`rel`$ saves the relative positions and distances that will again be needed in the loop calculating the forces and energies. Each thread has its own starting position $`iammyspace+1`$ in the shared memory space $`lstb`$ and these starting positions are uniformly distributed. This approach allows the different threads to work independently. The resulting speedup is much higher than the one that one would obtain by calculating in the parallel version an array $`lstb`$ that is identical to the one from the serial version. If there are on the average more neighbors than expected (24 by default) the allocated space becomes too small. In this case the array $`lstb`$ is deallocated and a new larger version is allocated. This check for sufficient memory requires some minimal amount of coordination among the processors and is implemented by a critical section. If a reallocation is necessary, a message is written into an file to alert the user of the inefficiency due to the need of a second calculation of the Verlet list.

```
        allocate(lsta(2,nat))
        nnbrx=24
2345    nnbrx=3*nnbrx/2
        allocate(lstb(nnbrx*nat),rel(5,nnbrx*nat))

        indlstx=0

!$omp parallel  &
!$omp private(iat,cut2,iam,ii,indlst,l1,l2,l3,myspace,npr) &
!$omp shared (indlstx,nat,nn,nnbrx,ncx,ll1,ll2,ll3,icell,lsta,lstb,lay, &
!$omp rel,rxyz,cut,myspaceout)

        npr=1
!$       npr=omp_get_num_threads()
        iam=0
!$       iam=omp_get_thread_num()

        cut2=cut**2
        myspace=(nat*nnbrx)/npr
        if (iam.eq.0) myspaceout=myspace
! Verlet list, relative positions
        indlst=0
      do 6000,l3=0,ll3-1
      do 6000,l2=0,ll2-1
      do 6000,l1=0,ll1-1
      do 6600,ii=1,icell(0,l1,l2,l3)
        iat=icell(ii,l1,l2,l3)
        if ( ((iat-1)*npr)/nat .eq. iam) then
        lsta(1,iat)=iam*myspace+indlst+1
        call sublstiat(iat,nn,ncx,ll1,ll2,ll3,l1,l2,l3,myspace, &
             rxyz,icell,lstb(iam*myspace+1),lay,rel(1,iam*myspace+1),cut2,indlst)
        lsta(2,iat)=iam*myspace+indlst
        endif

6600    continue
6000    continue
!$omp critical
        indlstx=max(indlstx,indlst)
!$omp end critical
!$omp end parallel

           if (indlstx.gt.myspaceout) then
               write(10,*) count,’NNBRX too  small’, nnbrx
               deallocate(lstb,rel)
               goto 2345
           endif

        subroutine sublstiat(iat,nn,ncx,ll1,ll2,ll3,l1,l2,l3,myspace, &
                   rxyz,icell,lstb,lay,rel,cut2,indlst)
        implicit real*8 (a-h,o-z)
        dimension rxyz(3,nn),lay(nn),icell(0:ncx,-1:ll1,-1:ll2,-1:ll3), &
                  lstb(0:myspace-1),rel(5,0:myspace-1)

        do 6363,k3=l3-1,l3+1
        do 6363,k2=l2-1,l2+1
        do 6363,k1=l1-1,l1+1
        do 6363,jj=1,icell(0,k1,k2,k3)
          jat=icell(jj,k1,k2,k3)
          if (jat.eq.iat) goto 6363
          xrel= rxyz(1,iat)-rxyz(1,jat)
          yrel= rxyz(2,iat)-rxyz(2,jat)
          zrel= rxyz(3,iat)-rxyz(3,jat)
          rr2=xrel**2 + yrel**2 + zrel**2
          if ( rr2 .le. cut2 ) then
           indlst=min(indlst,myspace-1)
           lstb(indlst)=lay(jat)
           tt=sqrt(rr2)
           tti=1.d0/tt
           rel(1,indlst)=xrel*tti
           rel(2,indlst)=yrel*tti
           rel(3,indlst)=zrel*tti
           rel(4,indlst)=tt
           rel(5,indlst)=tti
           indlst= indlst+1
          endif
6363        continue

        return
        end

```

## V Calculation of the energies and forces

The computationally most important part taking some 80 percent of the CPU time is the calculation of the energies and forces. The energy expression for the Lenosky force field is given by

$$E=\underset{i,j}{}\varphi (r_{ij})+\underset{i}{}U\left(\underset{j}{}\rho (r_{ij})+\underset{j,k}{}f(r_{ij})f(r_{ik})g(cos(\theta _{jik}))\right)$$
(1)
All the functions ($`\varphi `$, $`U`$, $`\rho `$, $`f`$, $`g`$) in this energy expression are given by cubic splines. The subroutine for evaluating the cubic spline is listed below. The case that the argument is outside the cubic spline interval $`[tmin,tmax]`$ is rare and unimportant for performance considerations. The important cubic spline case is characterized by many dependencies. In the case of such dependencies the latency of the functional unit pipeline comes into play and reduces the attainable speed . A latency of some 20 cycles comes from the first two statements ( tt=(x-tmin)\*hi ; klo=tt) alone, requiring arithmetic operations and a floating point to integer conversion. For this reason the calculation of tt was taken out of the (most likely occurring) else block to overlap its evaluation with the evaluation of the if clauses. To further speed up the evaluation of the splines the structure of the energy expression 1 was exploited. In the computationally most important loop over $`j`$ and $`k`$ two splines ($`f(r_{ik})`$ and $`g(cos(\theta jik))`$ ) have to be evaluated. Inlining by hand the subroutine splint for both evaluations and calculating alternatingly one step of the first spline evaluation and one step of the second spline evaluation introduces two independent streams. This reduces the effect of latencies and boosts speed. Compilers are not able to do these complex type of optimizations. The best performance after these optimizations was obtained with low level compiler optimization flags (-O3 -qarch=pwr3 -qtune-pwr3 on IBM Power3, -O2 on the Compaq DEC Alpha, -O2 -xW on Intel Pentium4 )

```
        subroutine splint(ya,y2a,tmin,tmax,hsixth,h2sixth,hi,n,x,y,yp)
        implicit real*8 (a-h,o-z)
        dimension y2a(0:n-1),ya(0:n-1)

! interpolate if the argument is outside the cubic spline interval [tmin,tmax]
        tt=(x-tmin)*hi
        if (x.lt.tmin) then
          yp=hi*(ya(1)-ya(0)) -  &
          ( y2a(1)+2.d0*y2a(0) )*hsixth
          y=ya(0) + (x-tmin)*yp
        else if (x.gt.tmax) then
          yp=hi*(ya(n-1)-ya(n-2)) +  &
          ( 2.d0*y2a(n-1)+y2a(n-2) )*hsixth
          y=ya(n-1) + (x-tmax)*yp
! otherwise evaluate cubic spline
        else
          klo=tt
          khi=klo+1
          ya_klo=ya(klo)
          y2a_klo=y2a(klo)
          b=tt-klo
          a=1.d0-b
          ya_khi=ya(khi)
          y2a_khi=y2a(khi)
          b2=b*b
          y=a*ya_klo
          yp=ya_khi-ya_klo
          a2=a*a
          cof1=a2-1.d0
          cof2=b2-1.d0
          y=y+b*ya_khi
          yp=hi*yp
          cof3=3.d0*b2
          cof4=3.d0*a2
          cof1=a*cof1
          cof2=b*cof2
          cof3=cof3-1.d0
          cof4=cof4-1.d0
          yt1=cof1*y2a_klo
          yt2=cof2*y2a_khi
          ypt1=cof3*y2a_khi
          ypt2=cof4*y2a_klo
          y=y + (yt1+yt2)*h2sixth
          yp=yp + ( ypt1 - ypt2 )*hsixth
        endif
      return
      end

```

The final single processor performance for the entire subroutine is 460 Mflops on a Compaq DEC Alpha at 833 MHz, 300 Mflops on a IBM Power3 at 350 MHz and 550 Mflops on a Pentium 4. In order to obtain a high parallel speedup in this central part of the subroutine the threads are completely decoupled. This was done by introducing private copies for each thread to accumulate the energies $`tener`$ and forces $`txyz`$. The global energy and force are summed up in an additional loop at the end of the parallel region in a critical section.

```
!$omp parallel  &
!$omp private(iam,npr,iat,iat1,iat2,lot,istop,tcoord,tcoord2, &
!$omp tener,tener2,txyz,f2ij,f3ij,f3ik,npjx,npjkx) &
!$omp shared (nat,nnbrx,lsta,lstb,rel,ener,ener2,fxyz,coord,coord2,istopg)

        npr=1
!$       npr=omp_get_num_threads()
        iam=0
!$       iam=omp_get_thread_num()

        npjx=300 ; npjkx=3000
        istopg=0

        if (npr.ne.1) then
! PARALLEL CASE
! create temporary private scalars for reduction sum on energies and
!        temporary private array for reduction sum on forces
!$omp critical
        allocate(txyz(3,nat),f2ij(3,npjx),f3ij(3,npjkx),f3ik(3,npjkx))
!$omp end critical
       if (iam.eq.0) then
        ener=0.d0
        ener2=0.d0
        coord=0.d0
        coord2=0.d0
        do 121,iat=1,nat
        fxyz(1,iat)=0.d0
        fxyz(2,iat)=0.d0
121     fxyz(3,iat)=0.d0
       endif

        lot=nat/npr+.999999999999d0
        iat1=iam*lot+1
        iat2=min((iam+1)*lot,nat)
        call subfeniat(iat1,iat2,nat,lsta,lstb,rel,tener,tener2, &
               tcoord,tcoord2,nnbrx,txyz,f2ij,npjx,f3ij,npjkx,f3ik,istop)
!$omp critical
        ener=ener+tener
        ener2=ener2+tener2
        coord=coord+tcoord
        coord2=coord2+tcoord2
        istopg=istopg+istop
        do 8093,iat=1,nat
        fxyz(1,iat)=fxyz(1,iat)+txyz(1,iat)
        fxyz(2,iat)=fxyz(2,iat)+txyz(2,iat)
        fxyz(3,iat)=fxyz(3,iat)+txyz(3,iat)
8093    continue
!$omp end critical
        deallocate(txyz,f2ij,f3ij,f3ik)

        else
! SERIAL CASE
        iat1=1
        iat2=nat
        allocate(f2ij(3,npjx),f3ij(3,npjkx),f3ik(3,npjkx))
        call subfeniat(iat1,iat2,nat,lsta,lstb,rel,ener,ener2, &
               coord,coord2,nnbrx,fxyz,f2ij,npjx,f3ij,npjkx,f3ik,istop)
        deallocate(f2ij,f3ij,f3ik)

        endif
!$omp end parallel

        if (istopg.gt.0) stop ’DIMENSION ERROR (see WARNING above)’
        ener_var=ener2/nat-(ener/nat)**2
        coord=coord/nat
        coord_var=coord2/nat-coord**2

        deallocate(rxyz,icell,lay,lsta,lstb,rel)

        end

        subroutine subfeniat(iat1,iat2,nat,lsta,lstb,rel,tener,tener2, &
               tcoord,tcoord2,nnbrx,txyz,f2ij,npjx,f3ij,npjkx,f3ik,istop)
        implicit real*8 (a-h,o-z)
        dimension lsta(2,nat),lstb(nnbrx*nat),rel(5,nnbrx*nat),txyz(3,nat)
        dimension f2ij(3,npjx),f3ij(3,npjkx),f3ik(3,npjkx)

     initialize data ........

! create temporary private scalars for reduction sum on energies and
        tener=0.d0
        tener2=0.d0
        tcoord=0.d0
        tcoord2=0.d0
        istop=0
        do 121,iat=1,nat
        txyz(1,iat)=0.d0
        txyz(2,iat)=0.d0
121     txyz(3,iat)=0.d0

! calculation of forces, energy

        do 1000,iat=iat1,iat2

        dens2=0.d0
        dens3=0.d0
        jcnt=0
        jkcnt=0
        coord_iat=0.d0
        ener_iat=0.d0
        do 2000,jbr=lsta(1,iat),lsta(2,iat)
        jat=lstb(jbr)
        jcnt=jcnt+1
        if (jcnt.gt.npjx) then
            write(6,*) ’WARNING: enlarge npjx’
            istop=1
        endif

        fxij=rel(1,jbr)
        fyij=rel(2,jbr)
        fzij=rel(3,jbr)
        rij=rel(4,jbr)
        sij=rel(5,jbr)

! coordination number calculated with soft cutoff between first and
! second nearest neighbor
        if (rij.le.2.36d0) then
        coord_iat=coord_iat+1.d0
        else if (rij.ge.3.83d0) then
        else
        x=(rij-2.36d0)*(1.d0/(3.83d0-2.36d0))
        coord_iat=coord_iat+(2*x+1.d0)*(x-1.d0)**2
        endif

! pairpotential term
        call splint(cof_phi,dof_phi,tmin_phi,tmax_phi, &
           hsixth_phi,h2sixth_phi,hi_phi,10,rij,e_phi,ep_phi)
        ener_iat=ener_iat+(e_phi*.5d0)
        txyz(1,iat)=txyz(1,iat)-fxij*(ep_phi*.5d0)
        txyz(2,iat)=txyz(2,iat)-fyij*(ep_phi*.5d0)
        txyz(3,iat)=txyz(3,iat)-fzij*(ep_phi*.5d0)
        txyz(1,jat)=txyz(1,jat)+fxij*(ep_phi*.5d0)
        txyz(2,jat)=txyz(2,jat)+fyij*(ep_phi*.5d0)
        txyz(3,jat)=txyz(3,jat)+fzij*(ep_phi*.5d0)

! 2 body embedding term
        call splint(cof_rho,dof_rho,tmin_rho,tmax_rho, &
             hsixth_rho,h2sixth_rho,hi_rho,11,rij,rho,rhop)
        dens2=dens2+rho
        f2ij(1,jcnt)=fxij*rhop
        f2ij(2,jcnt)=fyij*rhop
        f2ij(3,jcnt)=fzij*rhop

! 3 body embedding term
        call splint(cof_fff,dof_fff,tmin_fff,tmax_fff, &
             hsixth_fff,h2sixth_fff,hi_fff,10,rij,fij,fijp)

        do 3000,kbr=lsta(1,iat),lsta(2,iat)
        kat=lstb(kbr)
        if (kat.lt.jat) then
        jkcnt=jkcnt+1
        if (jkcnt.gt.npjkx) then
            write(6,*) ’WARNING: enlarge npjkx’
            istop=1
        endif

! begin optimized version
        rik=rel(4,kbr)
      if (rik.gt.tmax_fff) then
        fikp=0.d0 ; fik=0.d0
        gjik=0.d0 ;  gjikp=0.d0 ; sik=0.d0
        costheta=0.d0 ; fxik=0.d0 ; fyik=0.d0 ; fzik=0.d0
      else if (rik.lt.tmin_fff) then
        fxik=rel(1,kbr)
        fyik=rel(2,kbr)
        fzik=rel(3,kbr)
        costheta=fxij*fxik+fyij*fyik+fzij*fzik
        sik=rel(5,kbr)
        fikp=hi_fff*(cof_fff(1)-cof_fff(0)) -  &
             ( dof_fff(1)+2.d0*dof_fff(0) )*hsixth_fff
        fik=cof_fff(0) + (rik-tmin_fff)*fikp
        tt_ggg=(costheta-tmin_ggg)*hi_ggg
        if (costheta.gt.tmax_ggg) then
           gjikp=hi_ggg*(cof_ggg(8-1)-cof_ggg(8-2)) + &
                ( 2.d0*dof_ggg(8-1)+dof_ggg(8-2) )*hsixth_ggg
                 gjik=cof_ggg(8-1) + (costheta-tmax_ggg)*gjikp
        else
           klo_ggg=tt_ggg
           khi_ggg=klo_ggg+1
           cof_ggg_klo=cof_ggg(klo_ggg)
           dof_ggg_klo=dof_ggg(klo_ggg)
           b_ggg=tt_ggg-klo_ggg
           a_ggg=1.d0-b_ggg
           cof_ggg_khi=cof_ggg(khi_ggg)
           dof_ggg_khi=dof_ggg(khi_ggg)
           b2_ggg=b_ggg*b_ggg
           gjik=a_ggg*cof_ggg_klo
           gjikp=cof_ggg_khi-cof_ggg_klo
           a2_ggg=a_ggg*a_ggg
           cof1_ggg=a2_ggg-1.d0
           cof2_ggg=b2_ggg-1.d0
           gjik=gjik+b_ggg*cof_ggg_khi
           gjikp=hi_ggg*gjikp
           cof3_ggg=3.d0*b2_ggg
           cof4_ggg=3.d0*a2_ggg
           cof1_ggg=a_ggg*cof1_ggg
           cof2_ggg=b_ggg*cof2_ggg
           cof3_ggg=cof3_ggg-1.d0
           cof4_ggg=cof4_ggg-1.d0
           yt1_ggg=cof1_ggg*dof_ggg_klo
           yt2_ggg=cof2_ggg*dof_ggg_khi
           ypt1_ggg=cof3_ggg*dof_ggg_khi
           ypt2_ggg=cof4_ggg*dof_ggg_klo
           gjik=gjik + (yt1_ggg+yt2_ggg)*h2sixth_ggg
           gjikp=gjikp + ( ypt1_ggg - ypt2_ggg )*hsixth_ggg
        endif
      else
        fxik=rel(1,kbr)
        tt_fff=rik-tmin_fff
        costheta=fxij*fxik
        fyik=rel(2,kbr)
        tt_fff=tt_fff*hi_fff
        costheta=costheta+fyij*fyik
        fzik=rel(3,kbr)
        klo_fff=tt_fff
        costheta=costheta+fzij*fzik
        sik=rel(5,kbr)
        tt_ggg=(costheta-tmin_ggg)*hi_ggg
        if (costheta.gt.tmax_ggg) then
          gjikp=hi_ggg*(cof_ggg(8-1)-cof_ggg(8-2)) + &
                ( 2.d0*dof_ggg(8-1)+dof_ggg(8-2) )*hsixth_ggg
                gjik=cof_ggg(8-1) + (costheta-tmax_ggg)*gjikp
          khi_fff=klo_fff+1
          cof_fff_klo=cof_fff(klo_fff)
          dof_fff_klo=dof_fff(klo_fff)
          b_fff=tt_fff-klo_fff
          a_fff=1.d0-b_fff
          cof_fff_khi=cof_fff(khi_fff)
          dof_fff_khi=dof_fff(khi_fff)
          b2_fff=b_fff*b_fff
          fik=a_fff*cof_fff_klo
          fikp=cof_fff_khi-cof_fff_klo
          a2_fff=a_fff*a_fff
          cof1_fff=a2_fff-1.d0
          cof2_fff=b2_fff-1.d0
          fik=fik+b_fff*cof_fff_khi
          fikp=hi_fff*fikp
          cof3_fff=3.d0*b2_fff
          cof4_fff=3.d0*a2_fff
          cof1_fff=a_fff*cof1_fff
          cof2_fff=b_fff*cof2_fff
          cof3_fff=cof3_fff-1.d0
          cof4_fff=cof4_fff-1.d0
          yt1_fff=cof1_fff*dof_fff_klo
          yt2_fff=cof2_fff*dof_fff_khi
          ypt1_fff=cof3_fff*dof_fff_khi
          ypt2_fff=cof4_fff*dof_fff_klo
          fik=fik + (yt1_fff+yt2_fff)*h2sixth_fff
          fikp=fikp + ( ypt1_fff - ypt2_fff )*hsixth_fff
         else
              klo_ggg=tt_ggg
              khi_ggg=klo_ggg+1
           khi_fff=klo_fff+1
              cof_ggg_klo=cof_ggg(klo_ggg)
           cof_fff_klo=cof_fff(klo_fff)
              dof_ggg_klo=dof_ggg(klo_ggg)
           dof_fff_klo=dof_fff(klo_fff)
              b_ggg=tt_ggg-klo_ggg
           b_fff=tt_fff-klo_fff
              a_ggg=1.d0-b_ggg
           a_fff=1.d0-b_fff
              cof_ggg_khi=cof_ggg(khi_ggg)
           cof_fff_khi=cof_fff(khi_fff)
              dof_ggg_khi=dof_ggg(khi_ggg)
           dof_fff_khi=dof_fff(khi_fff)
              b2_ggg=b_ggg*b_ggg
           b2_fff=b_fff*b_fff
              gjik=a_ggg*cof_ggg_klo
           fik=a_fff*cof_fff_klo
              gjikp=cof_ggg_khi-cof_ggg_klo
           fikp=cof_fff_khi-cof_fff_klo
              a2_ggg=a_ggg*a_ggg
           a2_fff=a_fff*a_fff
              cof1_ggg=a2_ggg-1.d0
           cof1_fff=a2_fff-1.d0
              cof2_ggg=b2_ggg-1.d0
           cof2_fff=b2_fff-1.d0
              gjik=gjik+b_ggg*cof_ggg_khi
           fik=fik+b_fff*cof_fff_khi
              gjikp=hi_ggg*gjikp
           fikp=hi_fff*fikp
              cof3_ggg=3.d0*b2_ggg
           cof3_fff=3.d0*b2_fff
              cof4_ggg=3.d0*a2_ggg
           cof4_fff=3.d0*a2_fff
              cof1_ggg=a_ggg*cof1_ggg
           cof1_fff=a_fff*cof1_fff
              cof2_ggg=b_ggg*cof2_ggg
           cof2_fff=b_fff*cof2_fff
              cof3_ggg=cof3_ggg-1.d0
           cof3_fff=cof3_fff-1.d0
              cof4_ggg=cof4_ggg-1.d0
           cof4_fff=cof4_fff-1.d0
              yt1_ggg=cof1_ggg*dof_ggg_klo
           yt1_fff=cof1_fff*dof_fff_klo
              yt2_ggg=cof2_ggg*dof_ggg_khi
           yt2_fff=cof2_fff*dof_fff_khi
              ypt1_ggg=cof3_ggg*dof_ggg_khi
           ypt1_fff=cof3_fff*dof_fff_khi
              ypt2_ggg=cof4_ggg*dof_ggg_klo
           ypt2_fff=cof4_fff*dof_fff_klo
              gjik=gjik + (yt1_ggg+yt2_ggg)*h2sixth_ggg
           fik=fik + (yt1_fff+yt2_fff)*h2sixth_fff
              gjikp=gjikp + ( ypt1_ggg - ypt2_ggg )*hsixth_ggg
           fikp=fikp + ( ypt1_fff - ypt2_fff )*hsixth_fff
         endif
      endif
! end optimized version

        tt=fij*fik
        dens3=dens3+tt*gjik

        t1=fijp*fik*gjik
        t2=sij*(tt*gjikp)
        f3ij(1,jkcnt)=fxij*t1 + (fxik-fxij*costheta)*t2
        f3ij(2,jkcnt)=fyij*t1 + (fyik-fyij*costheta)*t2
        f3ij(3,jkcnt)=fzij*t1 + (fzik-fzij*costheta)*t2

        t3=fikp*fij*gjik
        t4=sik*(tt*gjikp)
        f3ik(1,jkcnt)=fxik*t3 + (fxij-fxik*costheta)*t4
        f3ik(2,jkcnt)=fyik*t3 + (fyij-fyik*costheta)*t4
        f3ik(3,jkcnt)=fzik*t3 + (fzij-fzik*costheta)*t4
        endif
3000        continue
2000        continue
        dens=dens2+dens3
        call splint(cof_uuu,dof_uuu,tmin_uuu,tmax_uuu, &
             hsixth_uuu,h2sixth_uuu,hi_uuu,8,dens,e_uuu,ep_uuu)
        ener_iat=ener_iat+e_uuu

! Only now ep_uu is known and the forces can be calculated, lets loop again
        jcnt=0
        jkcnt=0
        do 2200,jbr=lsta(1,iat),lsta(2,iat)
        jat=lstb(jbr)
        jcnt=jcnt+1
        txyz(1,iat)=txyz(1,iat)-ep_uuu*f2ij(1,jcnt)
        txyz(2,iat)=txyz(2,iat)-ep_uuu*f2ij(2,jcnt)
        txyz(3,iat)=txyz(3,iat)-ep_uuu*f2ij(3,jcnt)
        txyz(1,jat)=txyz(1,jat)+ep_uuu*f2ij(1,jcnt)
        txyz(2,jat)=txyz(2,jat)+ep_uuu*f2ij(2,jcnt)
        txyz(3,jat)=txyz(3,jat)+ep_uuu*f2ij(3,jcnt)

! 3 body embedding term
        do 3300,kbr=lsta(1,iat),lsta(2,iat)
        kat=lstb(kbr)
        if (kat.lt.jat) then
        jkcnt=jkcnt+1

        txyz(1,iat)=txyz(1,iat)-ep_uuu*(f3ij(1,jkcnt)+f3ik(1,jkcnt))
        txyz(2,iat)=txyz(2,iat)-ep_uuu*(f3ij(2,jkcnt)+f3ik(2,jkcnt))
        txyz(3,iat)=txyz(3,iat)-ep_uuu*(f3ij(3,jkcnt)+f3ik(3,jkcnt))
        txyz(1,jat)=txyz(1,jat)+ep_uuu*f3ij(1,jkcnt)
        txyz(2,jat)=txyz(2,jat)+ep_uuu*f3ij(2,jkcnt)
        txyz(3,jat)=txyz(3,jat)+ep_uuu*f3ij(3,jkcnt)
        txyz(1,kat)=txyz(1,kat)+ep_uuu*f3ik(1,jkcnt)
        txyz(2,kat)=txyz(2,kat)+ep_uuu*f3ik(2,jkcnt)
        txyz(3,kat)=txyz(3,kat)+ep_uuu*f3ik(3,jkcnt)

        endif
3300        continue
2200        continue

        tener=tener+ener_iat
        tener2=tener2+ener_iat**2
        tcoord=tcoord+coord_iat
        tcoord2=tcoord2+coord_iat**2

1000    continue

        return
        end

```

In addition to the energy and the forces the program still returns the coordination number as well as the variance of the energy per atom and the coordination number. The coordination number is calculated using a soft cutoff between the first and second nearest neighbor distance. These extra calculations are very cheap and not visible as an increase in the CPU time

## VI Parallel performance results

Table I shows the final overall speedups obtained by the program. The results were obtained for an 8000 atom system, but the CPU time per call and atom is nearly independent of system size.

Obtaining such high speedups was not straightforward. Only the Compaq Fortran90 compiler was able to use in the original version of the program the OpenMP ’parallel do’ directive to obtain a good speedup. Both the IBM compiler and the Intel compiler failed. In order to get the performances of Table I, it was necessary to encapsulate the workload of the different threads into the subroutines ’sublstias’ and ’subfen’, which amounts to doing the parallelization quasi by hand. Using allocatable arrays in connection with OpenMP turned also out to be tricky. Because of these problems, the parallelization was much more painful that one might expect for a shared memory model.

## VII Conclusions

The results show that simulations for very large silicon systems are feasible on relatively cheap serial or parallel computers accessible to a large number of researches.

I thank Tom Lenosky for his help in the implementation of the interatomic potential and Cyrus Umrigar for providing me with the timings on the IBM machine.


format: text

text:
# 1 Introduction

## 1 Introduction

HERA is an $`ep`$ collider which has especially high sensitivity to new particles coupling to lepton-quark pairs. In 1994-97 HERA collided 27.5 GeV positrons on 820 GeV protons. In 1998 the proton energy was raised to 920 GeV increasing the center-of-mass energy $`\sqrt{s}`$ from 300 GeV to 318 GeV. In 1998 and in the first months of 1999, HERA ran with electrons. In May 1999 HERA switched back to $`e^+p`$ collisions.

The three main colliding periods as well as the corresponding luminosities for each experiment are summarized in table 1.

Recent results on searches for new particles based on different models beyond the SM are reported here, as well as searches for rare event topologies with high $`P_T`$ leptons and missing momentum.

Several models have been tested analysing the full HERA-1 data corresponding to an integrated luminosity of 115 and 128 $`\text{pb}^1`$ for H1 and ZEUS respectively.

## 2 Contact Interactions

This indirect search for new processes allows to test the effects of new physics implying the exchange of particles with a mass higher than the center-of-mass energy $`s`$. An effective mass scale $`\mathrm{\Lambda }`$ is introduced: $`\mathrm{\Lambda }\sqrt{s}`$.

Different models beyond the Standard Model (SM) can thus be tested in the frame of contact interactions: for instance, in compositeness models, limits on the mass scale $`\mathrm{\Lambda }`$ can be set. Limits on high mass leptoquarks can also be determined via the ratio $`M/\lambda `$ where $`\lambda `$ is the coupling between leptoquarks and SM particles. Finally results from contact interactions allow also to put limits on the radius of the quark $`R_q`$.

Contact interactions are represented by lagrangian (1) which is added to the SM lagrangian.

$$_{CI}=\underset{a,b=L,R}{}\eta _{ab}^q(\overline{e_a}\gamma _\mu e_a)(\overline{q_b}\gamma _\mu q_b)$$
(1)
$`\eta _{ab}^q=ϵ\frac{g^2}{(\mathrm{\Lambda }_{ab}^q)^2}`$ are model-dependant coefficients depending on the quark flavour $`q`$ and the quark and electron chiralities $`a`$ and $`b`$. Different combinations of these coefficients lead to different contact interaction models. $`ϵ=\pm 1`$ is the sign of the interference with respect to the SM, $`\mathrm{\Lambda }_{ab}^q`$ is the effective mass scale, and $`g`$ is the coupling constant, conventionally fixed so that $`g^2=4\pi `$.

Contact interactions are searched for by analysing high $`P_T`$ samples of deep inelasting scattering candidates (DIS samples). Deviations from the SM in the distributions of usual kinematic variables of the DIS candidates (e.g. the virtuality $`Q^2`$ of the exchanged photon) are looked for.

As an example, figure 1 shows a fit of the measured neutral current (NC) $`Q^2`$ distribution to the VV compositeness model (this model assumes the same couplings $`\eta _{ab}^q`$ for all particle states), in both cases where new physics interferes constructively or destructively with the SM DIS processes. No deviation at high $`Q^2`$ is observed, therefore limits are deduced.

Figure 2 shows the H1 and ZEUS exclusion limits on the ratio $`ϵ/\mathrm{\Lambda }^2`$ for general compositeness models . H1 and ZEUS data are in good agreement with the SM expectation.

The resulting limits on the effective scale $`\mathrm{\Lambda }`$ are summarized in table 2.

Measurements of high $`Q^2`$ DIS can also be used to set limits on the radius of the quark $`R_q`$: the measurement of the spatial charge distribution of the quark is performed using the classical form factor approximation. In this approximation, the structure of the quark leads to the cross section deviation from the SM presented in equation (2).

$$\frac{d\sigma }{dQ^2}=\frac{d\sigma ^{SM}}{dQ^2}\left[1\frac{R_e^2}{6}Q^2\right]^2.\left[1\frac{R_q^2}{6}Q^2\right]^2$$
(2)
$`R_e`$ and $`R_q`$ are the root of the mean square radius of the electroweak charge of the electron and the quark respectively.

Assuming that electrons are point-like particles, an upper limit of $`\mathrm{0.73.10}^{18}`$ m is obtained for the quark radius from a fit to the ZEUS data .

## 3 Extra-Dimensions

The Large Extra-Dimensions (LED) model, proposed by Arkani-Hamed, Dimopoulos and Dvali assumes a space-time containing $`4+n`$ dimensions : the SM particles including the strong and electroweak bosons are supposed to be confined in a 4D world while gravitons can propagate in $`n`$ compactified extra-dimensions.

The effective coupling strength describing the interaction between the gravitons and SM particles is $`1/M_s`$, where $`M_s`$ is the fundamental Planck mass scale in the full space. If extra-dimensions are large enough, then $`M_s`$ may be in the TeV range. The contribution of graviton exchange to $`eq`$ scattering can be described by a contact interaction. Therefore the model is tested using NC samples and fitting the cross section to a formula including the graviton exchange, considering $`M_s`$ as a free parameter.

The corresponding limits on the Planck mass scale $`M_s`$ for both experiments are summarized in table 3 .

## 4 Leptoquarks

In various models (e.g. some models of Grand Unification), triplets of coloured bosons, called leptoquarks (LQs), are supposed to be produced, which could be observed at HERA by an $`eq`$ fusion as illustrated in figure 3.

LQs are scalar or vector bosons characterized by a weak Isospin $`T=0,1/2,1`$ and a fermionic number $`F=0,2`$ (defined as $`F=L+3B`$ where $`L`$ and $`B`$ are respectively the lepton and baryon numbers). The model proposed by Buchmüller, Rückl and Wyler (BRW)  describes 14 LQs (each coupling to only one combination of fermion chiralities):

* 7 LQs with a fermion number $`F=0`$ searched for in $`e^+p`$ collisions.
* 7 LQs with a fermion number $`|F|=2`$ searched for in $`e^{}p`$ data.

In the BRW model, first generation LQs decay into $`eq`$ and $`\nu q`$ with fixed branching ratios.

As the LQs decay topologies are identical to the Deep Inelastic Scattering (DIS) processes, LQs decaying into $`eq`$ are searched for in NC samples whereas LQs decaying into $`\nu q`$ are searched for in Charged Current (CC) samples.

The presence of LQs would manifest itself as an excess of DIS-like events with a large angle for the final state lepton. Therefore an angular cut varying with the LQ mass is applied in order to optimize the sensitivity to LQs. Mass distributions in the NC and CC channels are respectively shown in figures 4 and 5 .

In the NC channel, the previous excess observed at masses around 200 GeV in H1 and ZEUS 1994-97 data has not been confirmed by new 1999-2000 data. Therefore the mass distribution in the full data sample is compatible with the SM expectation. In the CC channel, both experiments show also a good agreement with the SM expectation. Therefore constraints on LQ masses and their coupling $`\lambda `$ to $`eq`$ pairs are derived.

Figure 6 presents the resulting limits on the coupling $`\lambda `$ as a function of the LQ mass for scalar LQs, in the context of the BRW model (fixed branching ratios). HERA limits are seen to be competitive with those obtained at other colliders. Moreover HERA extends the exclusion region at high mass of the LQ.

In the more general case where the branching ratio $`\beta _{eq}`$ ($`\beta _{\nu q}`$) for the LQ decay into $`eq`$ ($`\nu q`$) is not fixed anymore, constraints on $`\beta _{eq}`$ can be set as a function of the LQ mass, for different values of the coupling $`\lambda `$.

H1 results are shown in figure 7 for two example scalar LQs. This shows that HERA limits are compatible with Tevatron for low values of $`\lambda `$ and superior at large $`\lambda `$. Moreover HERA extends the limits to small values of $`\beta _{eq}`$.

## 5 Lepton-Flavour Violation

Experiments searching for neutrino oscillations show an increasing evidence for Lepton Flavour Violation, therefore the lepton number is likely to be violated. Processes like $`eq\mu q^{}`$ or $`eq\tau q^{}`$ via a LQ production or exchange are searched for (these processes are also illustrated by figure 3). Both these transitions have been looked for in various data sets, and constraints have been set .

As an example, constraints on the process $`eq\mu q^{}`$ which involves two couplings $`\lambda _{eq}`$ and $`\lambda _{\mu q}`$ are illustrated here.

At low LQ mass, upper limits on the product $`\lambda _{eq}\times \sqrt{BR_{LQ\mu q}}`$ can be deduced as a function of the LQ mass. Assuming $`\lambda _{eq}=\lambda _{\mu q}`$, the branching ratio $`BR_{LQ\mu q}`$ is fixed to 0.5 and limits on $`\lambda _{eq}`$ versus the LQ mass are shown in figure 8 . This figure shows that HERA limits are better than low-energy experiments for $`ed\mu b`$ transitions via the process $`B\mu e`$ and for a LQ mass lower than 280 GeV.

At high LQ mass, the cross section is proportional to $`\left(\frac{\lambda _{eq}\lambda _{\mu q}}{M^2}\right)^2`$ where $`M`$ is the LQ mass. Bounds on this ratios were also set in  and .

## 6 Squarks produced by R-parity violation

In the Minimal Supersymmetric Standard Model (MSSM) , R-parity is defined as $`R_p=(1)^{3B+L+2S}`$. SM particles have $`R_p=1`$ whereas their supersymmetric partners have $`R_p=1`$.

If $`R_p`$-violation would be allowed, superparticles could be produced at HERA by an $`eq`$ fusion: $`ep\stackrel{~}{q}X`$. This process is illustrated by figure 9. An $`R_p`$-violating coupling $`\lambda ^{}`$ is involved corresponding to an $`R_p`$-violation term $`\lambda _{ijk}^{}L_iQ_j\overline{D}_k`$ added to the lagrangian.

The $`e^+p`$ data can mainly probe the coupling $`\lambda _{1j1}^{}`$ via the process $`e^++d\stackrel{~}{u}_L`$.

Both decay modes of the squark are investigated at HERA: $`\mathit{}_p`$decays and “gauge” decays.

The $`\mathit{}_p`$decay mode is illustrated in figure 10. As the squark can decay into an electron and a quark or a neutrino and a quark, the main backgrounds are respectively NC and CC processes.

Example final states resulting from a squark “gauge” decay (mediated by a $`R_p`$-conserving coupling) are illustrated in figure 11. The table indicates all signatures considered, together with their corresponding backgrounds for $`e^+p`$ scattering. No deviation from the SM has been observed, therefore limits on the coupling $`\lambda ^{}`$ are derived.

.

The results have been interpreted in various SUSY models. The first one is the unconstrained MSSM. In this model, sfermion masses are considered as free parameters. The other free parameters are: the squark mass $`M_{\stackrel{~}{q}}`$, the soft SUSY breaking mass term $`M_2`$, the mixing mass term for Higgs doublets $`\mu `$ and the ratio of the vacuum expectation values of the 2 neutral Higgs bosons $`\mathrm{tan}\beta `$. A scan of the SUSY parameter space $`(M_2,\mu ,\mathrm{tan}\beta )`$ is performed. The resulting limits on $`\lambda ^{}`$ of the H1 and ZEUS experiments are shown in figure 12 .

As they are calculated for a large range in $`\mu `$ and $`M_2`$, these limits are almost model independent. They are seen to be more restrictive than the indirect bounds from Atomic Parity Violation (APV) for the coupling $`\lambda _{131}^{}`$, which would lead to stop production at HERA. For a given value of $`\lambda _{1j1}^{}=0.3`$, squark masses below 260 GeV are excluded.

Another SUSY model is the minimal Super-Gravity model (mSUGRA)  which contains the additional constraint of a universal mass $`m_0`$ introduced for all sfermions. In this model the electroweak symmetry breaking is supposed to be driven by radiative corrections. Therefore the free parameters are mainly $`m_0`$ and the unified gaugino mass $`m_{1/2}`$ at the GUT scale, the sign of $`\mu `$ and $`\mathrm{tan}\beta `$.

The resulting limits in the plane ($`m_0,m_{1/2}`$) are presented in figure 13. H1 limits are competitive with (better than) Tevatron for low (medium) $`\mathrm{tan}\beta `$. For the coupling $`\lambda _{131}^{}`$, H1 bounds are comparable to LEP limits111 LEP limits should not depend too much on $`\mathrm{tan}\beta `$ for small and medium values of $`\mathrm{tan}\beta `$. for intermediate values of $`m_0`$, $`m_{1/2}`$ and $`\mathrm{tan}\beta `$.

## 7 Excited fermions

In various extensions of SM, excited states of fermions are expected ($`e^{}`$, $`\nu ^{}`$ and $`q^{}`$). The excited fermions can then de-excite by emission of electroweak bosons $`\gamma `$, $`Z`$ or $`W`$ as shown in figure 14 producing various final state topologies to be searched for. This process is described by the effective lagrangian (3) .

$$_{int}=\frac{1}{2\mathrm{\Lambda }}\overline{F}_R^{}\sigma ^{\mu \nu }\left[gf\frac{\tau ^a}{2}W_{\mu \nu }^a+g^{}f^{}\frac{Y}{2}B_{\mu \nu }+g_sf_s\frac{\lambda ^a}{2}G_{\mu \nu }^a\right]F_L$$
(3)
$`F_R^{}`$ and $`F_L^{}`$ are the left-handed and right-handed components of the excited form weak isodoublets. $`W_{\mu \nu }^a`$, $`B_{\mu \nu }`$ and $`G_{\mu \nu }^a`$ are the field-strength tensors of the respective gauge groups SU(2), U(1) and SU(3). $`g`$, $`g^{}`$ and $`g_s`$ are the SM couplings. $`f`$, $`f^{}`$ and $`f_s`$ determine the coupling strength between the excited fermions and the bosons corresponding to the gauge groups U(1), SU(2) and SU(3). $`\tau ^a`$ are the Pauli matrices, $`Y`$ the weak hypercharge operator, $`\lambda ^a`$ the Gell-Mann matrices and $`\mathrm{\Lambda }`$ the compositeness scale.

Assuming relations between $`f`$, $`f^{}`$ and $`f_s`$, the branching ratios are fixed and limits on the ratio $`f/\mathrm{\Lambda }`$ versus the excited fermion mass can be deduced from measurements.

Figure 15 shows the mass distribution of the excited electron, neutrino and quark candidates in $`e\gamma `$, $`eW`$ and $`q\gamma `$ channels. . No deviation from the SM is observed, therefore limits are established.

The resulting limits on $`f/\mathrm{\Lambda }`$ for each channel are presented in figure 16 as a function of the excited fermion mass . In the excited electron channel, HERA limits are competitive with LEP limits for high excited electron masses. In the excited neutrino channel, HERA has a unique sensitivity at high excited neutrino masses. Finally, in the excited quark channel, HERA limits are more sensitive than Tevatron limits if small values of $`f_s`$ are considered.

Less model dependent limits are also obtained by allowing an arbitrary ratio $`f^{}/f`$ in the excited neutrino analysis. Figure 17 shows the limits on $`f/\mathrm{\Lambda }`$ as a function of the ratio $`f^{}/f`$ for different values of the excited neutrino mass. These limits are then converted into limits versus the excited neutrino mass for a range of $`f^{}/f`$ between $`5`$ and $`5`$.

## 8 Isolated leptons

High $`P_T`$ isolated lepton events with missing $`P_T`$ were observed in the 1994-97 H1 data . This analysis has now been extended to the full HERA-1 data. The most likely SM interpretation of these events is W production with a lepton decay of the W as illustrated in figure 18.

The LO cross section of the W production has been calculated and evaluated to be at the level of 1.3 pb. NLO calculations reduce the scale uncertainty to about $`10\%`$ at high $`P_T^X`$ .

W production events are characterized by missing $`P_T`$ (due to the presence of the neutrino), a high energy lepton and possibly a jet. The electron and muon channels of the W decay ($`We\nu `$ and $`W\mu \nu `$) have been analysed by H1 and ZEUS.

Table 4 shows the H1 and ZEUS event rates with the full HERA-1 statistics in both channels as a function of the transverse momentum of the hadronic system $`P_T^X`$ . The distribution of the lepton-neutrino transverse mass and the hadronic system transverse momentum $`P_T^X`$ are presented in figure 19. At low $`P_T^X`$ the observed events are consistent with the SM expectation. However more events than expected are observed by H1 at high $`P_T^X`$.

Figure 20 shows the 2-dimensional ($`P_T^X`$,$`M_T`$) distribution. Few data events are observed in regions where the SM expectation is low. A possible interpretation of these events is the anomalous production of single top events.

## 9 Single top production

The top quark production within the SM is negligible at HERA. Therefore if any top signal would be found, this would be a sign of new physics. The top quark could be produced by Flavour Changing Neutral Current (FCNC) through a $`\gamma u`$ fusion . In the semi-leptonic decay channel ($`tbW`$ with $`Wl\nu `$) the final state is composed of an isolated lepton, missing transverse momentum and high hadronic transverse momentum $`P_T^X`$.

The W selection is further tightened in order to select top candidates . Figure 21 presents the resulting mass distribution of the lepton-neutrino-jet system for the electron and muon channels. Several lepton events are in the top mass region.

A complementary top search has been performed via the hadronic channel ($`tbW`$ with $`Wq\overline{q}`$) with three high $`P_T`$ jets expected in the final state.

In addition to the existing analysis, further cuts are applied. A 2-jet mass is reconstructed in a W mass window in order to reduce the signal/noise ratio. A 3-jet mass is also reconstructed as shown in figure 22.

The results of the top searches are shown in table 5 for both semi-leptonic and hadronic decay channels. No visible excess is observed in the hadronic channel compared to the SM, therefore limits on anomalous top couplings are derived.

HERA resulting limits on $`tu\gamma `$ coupling are the following:

* H1: $`\kappa _{tu\gamma }<0.305`$
* ZEUS: $`\kappa _{tu\gamma }<0.19`$

Figure 23 summarizes the limits on the anomalous coulings $`V_Z`$ ($`tuZ`$ vectorial coupling) and $`\kappa _\gamma `$ ($`tu\gamma `$ magnetic coupling) obtained at HERA, LEP and Tevatron. HERA sensitivity to $`\kappa _\gamma `$ is competitive with other colliders.

## 10 Conclusion and future prospects

HERA is the unique collider to test direct $`eq`$ interactions. About $`110\text{pb}^1`$ of $`e^+p`$ data and $`15\text{pb}^1`$ of $`e^{}p`$ data have been collected per experiment at HERA-1.

No evidence of new physics has been observed in various models in inclusive analyses (contact interactions, extra-dimensions, leptoquarks) and exclusive analyses (lepton-flavour violation, $`R_p`$-violating SUSY, excited fermions), therefore new constraints have been set. HERA limits are seen to be competitive with and complementary to the LEP and Tevatron searches. The status of isolated lepton events with missing $`P_T`$ is still intriguing and will become clearer with the new HERA-2 data.

HERA has been shutdown since fall 2000 for a general upgrade: new focussing magnets have been installed in order to increase the luminosity and many improvements have been performed in the detectors in order to increase their sensitivity. Moreover, the lepton beam will be longitudinally polarised in the H1 and ZEUS interaction regions. The first luminosity runs are predicted for beginning 2002 and HERA-2 is expected to accumulate 1 $`\text{fb}^1`$ in the next 5 years. The anticipated factor of ten increase in the integrated luminosity will give an outstanding discovery potential for HERA.

## 11 Acknowledgements

I which to thank my H1 and ZEUS colleagues who contributed to the results presented here as well as people who helped me in preparing this talk.


text:
# Dislocation Emission around Nanoindentations on a (001) fcc Metal Surface Studied by STM and Atomistic Simulations

## Abstract

We present a combined study by Scanning Tunneling Microscopy and atomistic simulations of the emission of dissociated dislocation loops by nanoindentation on a (001) fcc surface. The latter consist of two stacking-fault ribbons bounded by Shockley partials and a stair-rod dislocation. These dissociated loops, which intersect the surface, are shown to originate from loops of interstitial character emitted along the \<110\> directions and are usually located at hundreds of angstroms away from the indentation point. Simulations reproduce the nucleation and glide of these dislocation loops.

Nanoindentation is a well-recognized technique that elucidates mechanical properties of surfaces and thin films at atomic scale . Although it is the small-scale counterpart of the traditional hardness test, it is far from being understood at the same level. Experimental evidence is still rather fragmentary although a number of studies involving atomic force microscopy and related techniques have provided information about such remarkable phenomena as the passivating properties of self-assembled molecular single layers  or the softening role of steps in the onset of plasticity . Although most of our initial knowledge of the processes taking place under a nanoindenter comes from molecular dynamics simulations, in which much progress has been made in recent years , a limiting factor of the comparison of atomistic simulations and experiment has been the different length scales accessible to each. In the present work we have succeeded to bridge the length-scale gap by increasing the experimental resolution around the nanoindentation point down to atomic scale.

The nucleation of discrete dislocations and their interactions are crucial to understand the first steps of crystal plasticity. In the present Letter, we report the first direct observation by Scanning Tunneling Microscopy (STM) of the emission of individual dislocation loops by nanoindentation of a Au(001) surface and show that a combined STM and atomistic simulation investigation explains how these loops are originated and how they glide away hundreds of angstroms from the indentation point. We fully characterize these loops and present evidence of the atomic displacement processes that are responsible for their generation. Furthermore, standard dislocation theory is shown to provide a suitable and convenient framework for encompassing the experimental results and the atomic simulations.

The experiments were performed in ultra high vacuum conditions on a Au(001) crystal cleaned by repeated sputtering and annealing in a system described in more detail elsewhere . Nanoindentations were done pushing the tungsten tip against the sample for distances of a few nanometers from the stable tunnel distance with the control feedback switched off. The last layer of the Au(001) surface is reconstructed with the well known “hex” $`5\times 20`$ reconstruction which appears in STM images as fringes oriented along a \<110\> direction. These fringes arise from a moiré-like effect between the topmost layer with an hexagonal arrangement and the lower bulk-like layer with square symmetry.

A typical image of the surface of the crystal after performing nanoindentations is shown in Fig. 1. Indentations themselves are imaged as multi-storied pits, and some material is shown as a pile-up surrounding the pits. Rows of bump-like features with a height of 0.6 $`\pm 0.1`$ Å (called from now on hillocks) are apparent along the compact \<110\> directions on the (001) surface, at distances of hundreds of angstroms from the indentation points. We have also observed similar hillocks in the course of ion irradiation on both Au(001) and Ag(001) followed by gentle annealing: in these cases the hillock distribution is random, with no alignments in rows as the ones surrounding nanoindentations.

Close-up STM images of hillocks are shown in Fig. 2 together with an interpretation of their sub-surface structure, which consists of two stacking-faults on intersecting {111} planes each bounded by two parallel Shockley partial dislocations. The whole configuration is held up by a stair-rod dislocation parallel to the surface. We argue that the origin of the hillocks can be traced to perfect dislocation loops that are punched into the crystal by the tip displacement into the surface . In our case, a family of such dislocation loops would consist of V-shaped half-loops intersecting the surface with Burgers vector parallel to the latter. They would, on energetic grounds, be split into two pairs of Shockley partials giving rise to the configuration shown in Fig.2b, i.e. resulting in hillocks. We can ascertain that these loops are of interstitial character based on the fact that for Au (Fig. 2a) we observe a missing reconstruction fringe on top of the hillock. As the interatomic spacing of the substrate is larger than the one of the uppermost reconstructed layer, the position of the missing-fringe in the moiré-like pattern corresponds to the position of the extra row of interstitial atoms below. Furthermore, for the Ag(001) atomic resolution image of Fig. 2c, it is clear that each partial dislocation produces a mismatch of one half interatomic unit between the rows on both sides of each stacking-fault, the two mismatches adding up to one extra row of atoms in the inner side of the hillock.

To gain insight into the atomic processes involved in the creation of these hillocks, atomistic simulations were carried out. A repulsive potential was used to model a spherical indenter of radius 4 nm penetrating the surface of a Au(001) crystal modeled with the embedded atom method potential . The reconstructed layer is thought to behave like a floating layer . Thus, it is not expected to affect dislocation generation and behavior during nanoindentation and was not included in the simulations. A top view of the surface after the simulated nanoindentation is shown in Fig. 3. In agreement with experiment, it is observed that two hillocks have been generated around the nanoindentation trace. To make a more quantitative comparison between simulations and dislocation models, a quantity closely related to the Burgers vector can be defined for every atom. This quantity is called the slip-vector, $`𝐬_i`$, being defined for atom $`i`$ as:

$$𝐬_i=\frac{1}{N_s}\underset{ji}{\overset{N_{nn}}{}}(𝐫_{ij}𝐫_{ij}^0)$$
(1)
where $`𝐫_{ij}`$ and $`𝐫_{ij}^0`$ are the vectors linking atom $`i`$ and all its $`N_{nn}`$ nearest neighbors $`j`$ in their current and reference (prior to the indentation) positions, respectively. $`N_s`$ stands for the number of slipped neighbors. The spatial distribution of the slip vectors moduli $`𝐬_i`$ around the nanoindentation trace, with a suitable color scale, is shown in Fig. 4. Green atoms with $`𝐬_i\frac{a_0}{\sqrt{6}}`$ are atoms on a stacking fault ($`a_0`$ is the Au lattice parameter). Blue atoms correspond to values of the slip-vector modulus between the stacking fault value and zero, while yellow ones have a slip-vector ranging between a complete lattice parameter and the stacking fault value. Blue and yellow atoms are, consequently, around the core of a Shockley partial dislocation. The configuration shown in Fig. 4, with a height displacement for the topmost atoms of 0.6 Å, strikingly reproduces the previously proposed subsurface structure of a hillock.

A dynamical picture of the process is obtained by recording successive frames of the simulated atomic events. To obtain the sequence, the indenter is first lowered in increments of 0.01 nm down to a depth of 0.58 nm in a quasi-static way at zero temperature. At this point dislocation loops are observed below the indenter, in agreement with previous results. Then, a constant energy molecular dynamics simulation follows the evolution of the system for 36 ps. A different dislocation configuration (again reproducing the previously proposed subsurface structure of a hillock) is created close to the tip and, then, glides away in a \<110\> direction. We stress that, although the experimentally observed hillocks are usually much larger than the simulated ones (and the indentation itself is also much deeper), we observe also hillocks of the same size in both experiment and simulation.

Hillocks are seen to glide as a whole unit. This behavior can be understood on the grounds that Shockley partial dislocations bounding a stacking fault are expected to glide easily on the {111} gliding planes. In our simulation they indeed glide away dragging with them the stair-rod dislocation (motion of structures formed by stair-rod dislocations and Shockley partial dislocations in thin films has been recently reported). We argue that the rows of hillocks appearing in the STM image of Fig. 1 are the result of successive emission of loops that glide away from the nanoindentation trace. Once started into motion due to the high stress close to the indentation, the hillocks would glide away from the indentation point until they collide or interact with other defects in the crystal. The simulation cell is too small to observe in detail this effect, although the hillocks move with a constant velocity within the unit cell once they are far enough from the indentation point. Hillock-like structures (although interpreted in a different way) have been observed to drift in highly stressed regions. The creation of hillocks following ion irradiation can also be explained in terms of the above model: after long ion irradiation and further annealing, the surface is known to exhibit a multi-storied pit structure, whereas the presence of a large supersaturation of adatoms and, probably, sub-surface interstitials, is suspected. These extra-atoms can cluster on crystallographic planes and, after relaxation, give rise to perfect loops that start the above mechanism.

In the course of the many simulations performed, a variety of hillock structures have been found. They differ in the exact arrangement of the loop below the surface: more complex configurations than the simplest one observed in Fig. 4 (an edge V-shaped loop dissociated along {111} planes) are possible. But the defining characteristics are the same for all of them: the Burgers vector of the undissociated loop is a lattice vector, all the sections of the loop appear dissociated along {111} planes, and the loop glides as a whole unit in a \<110\> direction. It is also worth remarking that, in the simulations, the hillocks remain in place once the tip is retracted from the surface, in agreement with the experiment.

Comparing different experimental images, it is observed that the size distribution of both the span $`s`$ (distance between stacking faults at the surface), and width $`w`$ (separation between partial dislocations again at the surface) of the observed hillocks is rather broad, ranging from a lowest resolvable size of about one reconstruction period ($``$ 1.4 nm) to a size of several ones. However, the parameter $`w`$ (in fact, the width of the extended dislocation) is experimentally seen to level off with increasing span $`s`$ of the hillock. The repulsive interaction force between two Shockley partial dislocation segments, at a given distance, increases initially with segment length (which is proportional to the observed span) but attains a constant value when this length becomes much larger than the distance between the segments. Within a model using elemental dislocation theory and taking into account image dislocations to include the effect of the free surface, the exact form of the $`w(s)`$ curve can be predicted and is found to be in agreement with the experimental data. It is worth emphasizing that the parameters (Burgers vectors, geometry…) of this new configuration can be explained in terms of continuum dislocation theory.

In summary, we have shown that the initial stages of plastic deformation around a nanoindentation result in the emission from near the contact point of dislocation half-loops intersecting the surface; they can, alternatively, be created by accretion of irradiated-in interstitials. These loops are split into pairs of Shockley partial dislocations giving rise to peculiar configurations at the surface (hillocks), involving four Shockley partial dislocations and a stair-rod which can be analyzed in terms of standard dislocation theory. Dissociated loops glide across compact planes and this provides a novel mechanism for matter transport away from nanoindentations. The observation of these defects in other systems, such as Ag(001), suggests that our results can be generalized to other (001) fcc metal surfaces, thus providing a consistent description with an unprecedented resolution of the incipient plastic deformation mechanisms. This might help the interpretation of elementary yielding events observed in previous nanoindentation investigations.

This research was supported by Spanish CICyT through Project No. PB96-0652 and by the Office of Basic Energy Sciences, Division of Materials Sciences, U.S. Department of Energy under contract No. DE-AC04-94AL85000. O.R.d.l.F gratefully acknowledges support from the Spanish MEC.


text:
# Single production of doubly charged Higgs bosons at hadron colliders

## Acknowledgements

Work is supported by the Academy of Finland under the contract no. 42328 and by RFFI grant 01-02-17152 (Russian Fund of Fundamental Investigations).


text:
# 1 The parton density 𝑓 (𝑓={𝑢_𝑣,𝑑_𝑣,𝑢̄,𝑑̄,𝑔}) of our model at a scale 𝜇₀≈0.168 GeV. The 𝑓(𝑥) in our model is simulated using Monte Carlo without any parameter. All we need in our model are principle of detailed balance and assumption of equal probability for every energy configuration of an 𝑛-parton Fock state.

## Abstract

Taking proton as an ensemble of quark-gluon Fock states and using the principle of detailed balance, we construct a simple statistical model for parton distribution of proton. The recent observed Bjorken-$`x`$ dependent light flavor sea quark asymmetry $`\overline{d}(x)\overline{u}(x)`$ can be well reproduced by Monte Carlo simulation as a pure statistical effect.

Parton Distribution of Proton in a Simple Statistical Model

Yong-Jun Zhang111zyj@pubms.pku.edu.cn

Department of Physics, Peking University, Beijing 100871, China and

Institute of High Energy Physics, CAS, Beijing 100039, China

Bing-Song Zou222zoubs@mail.ihep.ac.cn

CCAST (World Laboratory), P.O. Box 8730, Beijing 100080 and

Institute of High Energy Physics, CAS, P. O. Box 918(4), Beijing 100039, China

Li-Ming Yang

Department of Physics, Peking University, Beijing 100871, China

PACS: 12.40.Ee, 12.38.Lg, 14.20.Dh, 14.65.Bt

Proton is the simplest system in which the three colors of QCD neutralize into a colorless object, but its internal quark-gluon structure is still not well understood. The complication comes from the presence of sea quarks in the proton. In all global analyses of parton distribution in nucleons before 1990, a symmetric light-quark ($`\overline{u}`$, $`\overline{d}`$) sea was assumed, based on the usual assumption that the sea of quark-antiquark pairs is produced perturbatively from gluon splitting. However, a surprisingly large asymmetry between the $`\overline{u}`$ and $`\overline{d}`$ sea quark distributions in the proton has been observed in recent deep inelastic scattering and Drell-Yan experiments.

There have been many theoretical attempts trying to find the origins for this asymmetry. It is believed that the asymmetry cannot be produced from perturbative QCD and mesonic degrees of freedom play an important role for the effect.

In this paper we follow a new idea for the origin of the light flavor sea quark asymmetry to reproduce the recent observed $`\overline{d}(x)\overline{u}(x)`$ distribution with a simple statistical model. The basic idea in Ref. is rather simple: while sea quark-antiquark pairs are produced flavor blindly by gluon splitting, $`\overline{u}`$ quarks have larger probability to annihilate than $`\overline{d}`$ quarks due to the fact that there are more $`u`$ quarks than $`d`$ quarks in the proton, which hence causes the asymmetry. Taking proton as an ensemble of quark-gluon Fock states and using the principle of detailed balance for transitions between various Fock states through creation or annihilation of partons, the probabilities $`\rho _{i,j,k}`$ of finding the quark-gluon Fock states $`|\{uud\}\{i,j,k\}`$ have been obtained and given in Table 2 of Ref., with $`i,j,k`$ the number of $`\overline{u}u`$ pairs, the number of $`\overline{d}d`$ pairs, the number of gluons, respectively. With the density matrix $`\rho _{i,j,k}`$ for the quark-gluon Fock states $`|\{uud\}\{i,j,k\}`$, the sea-quark flavor asymmetry was calculated as $`\overline{d}\overline{u}0.124`$, which is in surprisingly agreement with the experimental data $`\overline{d}\overline{u}=0.118\pm 0.012`$. Encouraged by this success, here we want to extend the model to calculate the Bjorken-$`x`$ distribution of partons to study the $`x`$-dependence of the flavor asymmetry in the nucleon sea, $`\overline{d}(x)\overline{u}(x)`$, which has recently been well measured by the FNAL E866/NuSea Collaboration.

For a quark-gluon Fock state $`|\{uud\}\{i,j,k\}`$, the total number of partons is $`n=3+2i+2j+k`$. If the $`n`$ partons were free particles without mutual interactions, their momentum distribution $`d\rho _n^F(p_1,\mathrm{},p_n)`$ would simply follow the $`n`$-body phase space, i.e.,

$$d\rho _n^F(p_1,\mathrm{},p_n)=d\mathrm{\Phi }_n(P;p_1,\mathrm{},p_n)=\delta ^4(P\underset{i=1}{\overset{n}{}}p_i)\underset{i=1}{\overset{n}{}}\frac{d^3p_i}{(2\pi )^32E_i}$$
(1)
with $`P`$ the 4-momentum of the proton, $`p_1,\mathrm{},p_n`$ the 4-momenta of $`n`$ partons. If we ignore the mass of partons, then $`E_i\sqrt{\stackrel{}{p}^2+m_i^2}=|\stackrel{}{p_i}|P_i`$ and we have

$$d\rho _n^F(p_1,\mathrm{},p_n)=\delta ^4(P\underset{i=1}{\overset{n}{}}p_i)\underset{i=1}{\overset{n}{}}\frac{P_idP_id\mathrm{\Omega }_i}{2(2\pi )^3}=\delta ^4(P\underset{i=1}{\overset{n}{}}p_i)\underset{i=1}{\overset{n}{}}\frac{E_idE_id\mathrm{\Omega }_i}{2(2\pi )^3}.$$
(2)
However, we know that partons are not free particles and are confined in the proton. In potential picture, partons are almost free only near the center of the proton and their momenta decrease when moving away from the center. Partons with smaller momenta at the center stay longer close to the center while partons with larger momenta at the center stay shorter around the center. So compared with total free particles, partons in the proton should have larger probability for smaller momenta. In this paper we assume the $`n`$-parton momentum distribution $`d\rho _n(p_1,\mathrm{},p_n)`$ to be

$`d\rho _n(p_1,\mathrm{},p_n)`$ $`=`$ $`{\displaystyle \frac{1}{_{i=1}^nP_i}}d\mathrm{\Phi }_n(P;p_1,\mathrm{},p_n)`$ (3)
$`=`$ $`\delta ^4(P{\displaystyle \underset{i=1}{\overset{n}{}}}p_i){\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{dP_id\mathrm{\Omega }_i}{2(2\pi )^3}}=\delta ^4(P{\displaystyle \underset{i=1}{\overset{n}{}}}p_i){\displaystyle \underset{i=1}{\overset{n}{}}}{\displaystyle \frac{dE_id\mathrm{\Omega }_i}{2(2\pi )^3}}.`$
This is equivalent to assuming equal probability for any energy configuration $`(E_1,\mathrm{},E_n)`$ of $`n`$ partons in the proton.

With the $`n`$-parton momentum distribution $`d\rho _n(p_1,\mathrm{},p_n)`$ and taking the light front formula for the Bjorken-$`x`$

$$x=\frac{E_{\mathrm{parton}}p_{z\mathrm{parton}}}{M_{\mathrm{proton}}},$$
(4)
we can easily get the $`x`$-distribution of partons, $`\rho _n(x)`$, for an $`n`$-parton Fock state of the proton by using a simple Monte Carlo simulation with a Monte Carlo event generator program called GENEV from the CERN computer program library (CERNLIB). The $`\rho _n(x)`$ is normalized as $`_0^1𝑑x\rho _n(x)=1`$.

For a $`n`$-parton Fock state $`|\{uud\}\{i,j,k\}`$, the total number of partons is $`n`$ including $`u,d,\overline{u},\overline{d}`$ and $`g`$ with parton number of $`2+i`$, $`1+j`$, $`i`$, $`j`$ and $`k`$, respectively. Then the $`x`$-distribution for each kind of partons, $`u_{i,j,k}(x)`$, $`d_{i,j,k}(x)`$, $`\overline{u}_{i,j,k}(x)`$, $`\overline{d}_{i,j,k}(x)`$ and $`g_{i,j,k}(x)`$, are

$`u_{i,j,k}(x)`$ $`=`$ $`\rho _n(x)(2+i),`$ (5)
$`d_{i,j,k}(x)`$ $`=`$ $`\rho _n(x)(1+j),`$ (6)
$`\overline{u}_{i,j,k}(x)`$ $`=`$ $`\rho _n(x)i,`$ (7)
$`\overline{d}_{i,j,k}(x)`$ $`=`$ $`\rho _n(x)j,`$ (8)
$`g_{i,j,k}(x)`$ $`=`$ $`\rho _n(x)k.`$ (9)
Summing up parton $`x`$-distribution of all possible Fock states $`|\{uud\}\{i,j,k\}`$ according to their weights ($`\rho _{i,j,k}`$) in Table 2 of Ref. to get the parton $`x`$-distribution of proton, we have

$`u(x)={\displaystyle \underset{i,j,k}{}}\rho _{i,j,k}u_{i,j,k}(x),`$ (10)
$`d(x)={\displaystyle \underset{i,j,k}{}}\rho _{i,j,k}d_{i,j,k}(x),`$ (11)
$`\overline{u}(x)={\displaystyle \underset{i,j,k}{}}\rho _{i,j,k}\overline{u}_{i,j,k}(x),`$ (12)
$`\overline{d}(x)={\displaystyle \underset{i,j,k}{}}\rho _{i,j,k}\overline{d}_{i,j,k}(x),`$ (13)
$`g(x)={\displaystyle \underset{i,j,k}{}}\rho _{i,j,k}g_{i,j,k}(x),`$ (14)
which satisfy normalization condition

$`{\displaystyle _0^1}x\left[u(x)+d(x)+\overline{u}(x)+\overline{d}(x)+g(x)\right]𝑑x=1.`$ (15)
Our $`u(x)`$ and $`d(x)`$ include both valence and intrinsic sea quarks which are identical and not distinguishable in our approach. The $`x`$-distribution of valance quarks ($`u_v`$, $`d_v`$) can be easily obtained as

$`u_v(x)`$ $`=`$ $`u(x)\overline{u}(x),`$ (16)
$`d_v(x)`$ $`=`$ $`d(x)\overline{d}(x).`$ (17)
In addition, we have

$`\overline{n}`$ $`=`$ $`{\displaystyle _0^1}\left[u(x)+d(x)+\overline{u}(x)+\overline{d}(x)+g(x)\right]𝑑x=5.57,`$ (18)
$`\overline{E}`$ $`=`$ $`{\displaystyle \frac{M_{\mathrm{Proton}}}{\overline{n}}}={\displaystyle \frac{0.938\mathrm{GeV}}{5.57}}=0.168\mathrm{GeV},`$ (19)
where $`\overline{n}`$ is the average number of partons in proton and $`\overline{E}`$ is the average energy of partons in proton.

From above equations, we get the parton $`x`$-distribution densities $`f(x,\mu _0^2)`$ $`(f=u,d,\overline{u},\overline{d},g)`$ as shown in Fig.1 with the scale $`\mu _0\overline{E}=0.168`$ GeV.

The corresponding momentum distribution, $`xf(x,\mu _0^2)`$, is shown in Fig.2. As a comparison, the $`xf(x,Q^2)`$ distribution of GRV at some higher $`Q^2`$ is shown in Fig.3. It will be interesting to check whether with some QCD evolution equation our $`xf(x,\mu _0^2)`$ distribution could evolve into a momentum distribution $`xf(x,Q^2)`$ at higher $`Q^2`$ to be similar to the GRV’s.

The quarks and gluons in the Fock states are the “intrinsic” partons of the proton, multi-connected non-perturbatively to the valence quarks. Such partons are different from “extrinsic” partons generated from the QCD hard bremsstrahlung and gluon-splitting as part of the lepton scattering interaction. Partons measured at certain $`Q^2`$ by experiments include both “intrinsic” and “extrinsic” ones. Since “extrinsic” partons are generated flavor blindly, the light flavor sea quark asymmetry is mainly due to “intrinsic” partons and is practically $`Q^2`$ independent, although correlations between “intrinsic” and “extrinsic” partons can cause some small $`Q^2`$ dependence for the asymmetry. Experimental data at various $`Q^2`$ values also show little $`Q^2`$ dependence. Hence we can compare our model prediction of $`\overline{d}(x)\overline{u}(x)`$ at the low scale $`\mu _0`$ directly with experimental data at higher scales, as shown in Fig.4. Our prediction is in good agreement with recent experiment data of FNAL E866/NuSea at $`Q^2=54`$ GeV2 and HERMES at $`<Q^2>=2.3`$ GeV2.

In summary, following Ref. we take proton as an ensemble of quark-gluon Fock state, i.e., $`|p=_{i,j,k}c_{i,j,k}|\{uud\}\{i,j,k\}`$, with $`\rho _{i,j,k}|c_{i,j,k}|^2`$ determined by the principle of detailed balance. By further assuming equal probability for any energy configuration ($`E_1,\mathrm{},E_n`$) of $`n`$-parton Fock state in the proton, we get parton $`x`$-distribution functions at a scale of $`\mu _0\overline{E}=0.168`$ GeV with a Monte Carlo simulation of a simple statistical model. The corresponding light flavor sea quark asymmetry $`\overline{d}(x)\overline{u}(x)`$ reproduces the recent experiment data quite well. This is a further support of the new origin of the light flavor sea quark asymmetry as a pure statistical effect due to the fact that there are more $`u`$-quarks than $`d`$-quarks in the proton.

Acknowledgment: The work is partially supported by the CAS Knowledge Innovation Project (KJCX2-N11) and by National Natural Science Foundation of China.


text:
# Automating the Synthetic Field Method: Application to Sextans A

## 1. Introduction

To observe the opacity and distribution of dust independent of theoretical models for the light or dust distribution, the dust needs to be backlighted by a more distant source. The group led by Keel and White used an occulted galaxy for this purpose and their studies (Domingue et al. (1999, 2000), Keel et al. (2001) and White et al. (1992, 2000)) showed insight into the fine structure of the dust in a few galaxies.

Gonzalez et al. (1998) developed a method to use the distribution of field galaxies which they called the “Synthetic Field Method” (SFM); it is not restricted to the case of nearby galaxy pairs and it can supply values for the dust opacities of galaxies in the disks of spirals and irregulars in general.

## 2. The Synthetic Field Method

The method consists of comparing the numbers of the background galaxies found in Wide Field Camera images of nearby spiral and irregular galaxies with those of simulations. These simulations are a Hubble Deep Field image, extincted by a certain opacity added to the original data. (see figure 1.) By adding the extincted HDF frame to the original data, all the effects of clustering and confusion by the foreground galaxy are calibrated. All the three HDF frames are used in seperate simulations to account as much as possible for clustering biases; however, clustering of the field galaxies remains the main uncertainty. Visual identification of the background galaxies in both the data and the simulations is very time consuming. Therefore we have searched for ways to automate the procedure.

### 2.1. Automating the SFM

Sextans A was the test case for automating the SFM together with NGC4536. In both cases results from manual galaxy counts were available. The automation works briefly as follows: a modified Source Extractor v2.2.2 (Bertin & Arnouts 1996) provides catalogues of objects with a variety of structural parameters. (Source Extractor’s neural network classification, concentration, contrast, flux radius, area, ellipticity, color and asymmetry). 111Concentration, contrast and asymmetry were defined as in Abrahams (1993).

Candidate objects are selected based on a score which is determined by a series of criteria for the structural parameters; values typical for background galaxies earn an object a higher score while those typical for stars and globular clusters result in a lower score. The advantages of the system are that objects does not need to meet all the criteria to be selected, and criteria can involve combinations of parameters.

The selected objects in the original data were checked visually to exclude diffraction spikes and blended stars. The selected objects in the simulations were automatically checked for the data-objects, thus leaving only simulated galaxies.

## 3. Sextans A data

Sextans A provided an excellent opportunity to test the automated Synthetic Field Method with data on most of the galaxy with varying exposure times but with the two filters (V (F555W) and I (F814W)) used by Gonzalez et al. (1998) (tabel 1).

HI maps from Skillmann et al. (1988; figure 3) showed that the HST fields covered a range of HI column densities; figure 3) This allows us to compare dust to gas ratios similar to the approach of Cuillandre et al. (1999) in the outskirts of M31.

## 4. Results

The results for the WF3 and 4 are consistent with 0 magnitude opacity and WF2 with 0.5 magnitude. (see figure 2) However the error in the number of background galaxies due to their clustering is such that all can be considered consistent with no dust opacity. As can be seen from figure 4 the determined opacities for Sextans A all correspond to the relation Cuillandre et al. (2001) found experimentally in Sextans A.

## 5. Discussion

We conclude first that the SFM can be applied with almost the same accuracy on fields with widely different exposure times. The limitations come from background galaxy clustering and confusion, not sensitivity. Secondly, the SFM can be automated to a satisfactory degree. As the clustering is the main uncertainty, applying the SFM to larger areas with similar properties (HI column density, inside or outside the spiral arm or at a certain radius) will lower these errors and give a general picture of dust distribution in irregular and spiral galaxies.

We intend now to carry out a survey using the SFM of suitable WFPC data currently in the HST archive.

## References

Bertin, E. and Arnouts, S., 1996, AAPS, 117, 393

Cuillandre, J., Lequeux, J., Allen, R.J., Mellier, Y. and Bertin, E., 2001,
ApJ, 554, 190

Domingue, D.L., Keel, W.C., Ryder, S.D. and White, R.E., 1999, AJ, 118, 1542

Domingue, D.L., Keel, W.C. and White, R.E., 2000, ApJ, 545, 171

Fruchter, A. and Hook, R.N., Proc. SPIE, 3164, 120

González, R.A., Allen, R.J., Dirsch, B., Ferguson, H.C., Calzetti, D.
and Panagia, N.,1998, ApJ, 506, 152

Keel, W.C. and White, R.E., 2001, AJ, 121, 1442

Skillman, E.D., Terlevich, R., Teuben, P.J. and van Woerden, H., 1988,
AAP, 198,33

White, R.E. and Keel, W. C., 1992, Nat, 359, 129

White, R.E., Keel, W.C. and Conselice, C.J., 2000, ApJ, 542, 761


text:
# Large scale extragalactic jets powered by very-high energy gamma rays

## Abstract

The radiative cooling of electrons responsible for the nonthermal synchrotron emission of large scale jets of radiogalaxies and quasars requires quasi continuous (in time and space) production of relativistic electrons throughout the jets over the scales exceeding 100 kpc. While in the standard paradigm of large scale jets this implies in situ acceleration of electrons, in this letter we propose a principally different “non-acceleration” origin of these electrons, assuming that they are implemented all over the length of the jet through effective development of electromagnetic cascades initiated by extremely high energy $`\gamma `$-rays injected into the jet from the central object. This scenario provides a natural and very economic way to power the jets up to distances of 100 kpc and beyond.

Although there is little doubt in nonthermal origin of radiation of large scale extragalactic jets, it remains a theoretical challenge to explain how the relativistic electrons responsible for the radio-to-X-ray synchrotron emission could be distributed quite uniformly all over the huge length of the jet $`L_{jet}100`$ kpc. Although the recent exciting discoveries by the Chandra observatory added much to our knowledge of structures of these fascinating components of powerful radiogalaxies and quasars, they didn’t solve the old problems and, if fact, brought new puzzles. In particular, there remain yet substantial problems concerning the identification of the X-ray emission mechanism. The most natural mechanism, the synchrotron radiation of multi-TeV electrons, applied in a standard manner to large scale jets, faces certain difficulties to interpret the observed X-ray features of a significant fraction of the Chandra jets (see e.g. Ref ). Therefore the recent “trend” which invokes relativistic bulk motion with Lorentz factor as large as $`\mathrm{\Gamma }10`$ (on $``$ 100 kpc scales !), and thus increases significantly the efficiency of the inverse Compton scattering on 2.7 K cosmic background radiation (CMBR) as a source of X-rays , was readily accepted by the extragalactic jet community . The very idea of relativistic bulk motion could be very productive also for other possible jet models . In particular, it would allow explanation of the broad-band data of the jet of 3C 273 by a single population of electrons . Even so, this effect does not solve the problem of continuous (in space and time) acceleration of multi-TeV electrons. The association of the observed X-rays to very high-energy (VHE) protons could be a possible solution .

In this letter we suggest a “non-acceleration” version of the electron synchrotron model, namely assuming that electrons are not accelerated in the jet, but are result of pair production of extremely high energy $`\gamma `$-rays interacting with the CMBR . The pair production in the field of CMBR and diffuse infrared background photons as a source of ultra-relativistic electrons created far from the central object (AGN) has been suggested  in the context of formation of giant extragalactic pair halos. The same process can provide ultrarelativistic electrons for nonthermal synchrotron X- and $`\gamma `$-rays in clusters of galaxies and beyond . Recently, highly collimated neutral beams (neutrons and $`\gamma `$-rays) have been noticed as possible supplier of synchrotron-emitting electrons in the extended jets of FRII radio galaxies .

The mean free path of $`\gamma `$-rays in the field of CMBR has a minimum $`L_\gamma 8`$ kpc at $`E_\gamma 10^{15}\mathrm{eV}`$. At both lower and higher energies $`L_\gamma `$ increases – sharply (exponentially) at $`E_\gamma 10^{15}\mathrm{eV}`$ (due to the threshold effect $`\gamma `$-rays interact with the Wien tail of CMBR), and slowly, $`L_\gamma 14.6E_{\gamma ,16}T_{2.7}^2[1+0.7\mathrm{ln}(E_{\gamma ,16}T_{2.7})]^1`$ kpc at $`E_\gamma 10^{15}`$ eV due to the decrease of the cross-section with the parameter $`E_\gamma E_{\mathrm{CMBR}}m_\mathrm{e}^2c^4`$. Hereafter, $`E_{\gamma ,16}=E_\gamma /10^{16}`$ eV, $`T_{2.7}=T/2.7`$ K. Thus, a $`\gamma `$-ray beam with a broad spectrum extending to $`10^{18}\mathrm{eV}`$ can supply the jet with “desirable” VHE electrons along the jet, and thus power the jet up to distances $`100`$ kpc or even more.

In what follows we analyze a simple model in which a VHE photons are injected into a cylindrical jet and, interacting with CMBR as well as with the low-frequency radiation of the jet itself, initiate an electromagnetic cascade, the relevant processes being pair-production (PP), inverse Compton scattering (ICS) and synchrotron radiation (see Fig. 1). For numerical calculations we use the code described in Ref. which solves kinetic equations for VHE particles cascading on soft photon backgrounds.

The first generation of pair-produced electrons upscatter the soft background photons. Propagation length of electron with respect to ICS in Klein-Nishina regime is $`L_{e,\mathrm{KN}}1.5T_{2.7}^2E_{e,16}`$ kpc. The electrons suffer also synchrotron losses.

High resolution radio observations show that the magnetic field in large scale jets may dominate by the regular component parallel to the jet axis, although there could be anomalous regions near the knots with oblique field components (see e.g. Ref.). Here we adopt a simplified picture assuming that the field consists of two - random and regular

(aligned with the jet axis) - components, and that the regular component significantly exceeds the random component. We normalize these fields as $`B_0=B_{0,5}\times 10^5`$ G and $`B_r0.01B_0=B_{\mathrm{r},7}\times 10^7`$ G, respectively. Generally, if the first generation electrons appear at small angle to the regular magnetic field (i.e. the $`\gamma `$-ray beam from the central engine is directed along the jet), they will spiral with a small pitch angle $`\theta `$ in the magnetic field $`B_0`$ and at the same time get gradually deflected by the random magnetic field $`B_\mathrm{r}`$. As long as the deflection of electrons does not exceed a few degree, the synchrotron losses of electrons are dominated by the random magnetic field. At larger angles the synchrotron radiation in the regular field becomes more important. Moreover, for $`\theta 30^{}`$ the latter dominates also over the Compton losses. At the same time, because of the Klein-Nishina effect, the electrons of extremely high energies, $`E10^{17}\mathrm{eV}`$ are basically cooled through the synchrotron radiation, even when they move at small pitch angles.

Deflection of electrons in the jet can be approximately described as diffusion in pitch angle . The diffusion length $`L_{\mathrm{diff}}`$ can be estimated as  $`L_{\mathrm{diff}}5B_{\mathrm{r},7}^2E_{e,14}^2l_{\mathrm{pc}}^1(\theta _0/3^{})^2\text{ kpc},`$ where $`l_{\mathrm{pc}}`$ is the correlation length of the random B-field in parsecs and $`\theta _0`$ is the initial opening angle of the primary photon beam. When electrons cool down to $`10^{14}`$ eV, their propagation length is determined by ICS $`L_e=0.3(T_{2.7})^4(E_{e,14})^1\text{ kpc}.`$ When $`L_{\mathrm{diff}}`$ becomes of order of $`L_e`$ the cascade electrons are essentially deflected by the random B-field. Comparing the diffusion length to the propagation elngth we find that the electron trajectories below the energy $`E_{\mathrm{crit}}=4\times 10^{13}B_{\mathrm{r},7}^{2/3}l_{\mathrm{pc}}^{1/3}T_{2.7}^{4/3}(\theta _0/3^{})^{2/3}\text{ eV}`$ are randomized. Such electrons are cooled through synchrotron and ICS losses, and form a bright ”shell” around the ”main stream” of the cascade as it is demonstrated in Fig. 1.

The energy losses of ”shell” electrons are dominated

by the regular magnetic field which establish a standard $`dn_e/dE_eE_e^2`$ type spectrum of electrons. Correspondingly, the spectrum of synchrotron radiation behaves as $`\nu F_\nu \nu ^{0.5}`$ with a maximum at $`h\nu _{\mathrm{max}}=1.6B_{0,5}(E_{\mathrm{crit}}/10^{14}\mathrm{eV})^2\text{ keV}.`$ Note that at low energies a possible energy dependent-escape might make the electron spectrum, and subsequently the synchrotron radio spectrum somewhat steeper.

Because the spectra of Chandra jets typically extend beyond 1 keV, for a reasonable value of $`E_{\mathrm{crit}}10^{14}\mathrm{eV}`$, which is a free parameter in our model, the regular B-field should exceed $`10^5\mathrm{G}`$. This explains the need in the regular field in our model. We need sufficiently strong magnetic filed in order to explain the X-ray data. The random field cannot play this role because such a strong random field would destroy the cascade from very beginning. On the other hand we need a cascade in order to have a more or less homogeneous distribution of electrons, both in density and energy spectrum, across the jet. Assuming that the maximum energy of synchrotron radiation from the jet lies in the X-ray energy domain, which minimizes the energy requirement to the source, we find the following relation between the random and regular magnetic fields: $`B_{\mathrm{r},7}^2B_{0,5}^{3/2}l_{\mathrm{pc}}7T_{2.7}^4(\theta _0/3^o)^2\left(h\nu _{\mathrm{max}}/1\mathrm{k}\mathrm{e}\mathrm{V}\right)^{3/2}`$. An example of numerical calculation of evolution of photon and electron spectra in the ”main stream” and in the ”shell” are shown in Fig. 2.

So far we have assumed that the synchrotron background in the jet is negligible compared to the external CMBR background. This seems to be the case for the jet in PKS 1127-145 (see Fig. 3) because the radio power of the jet is small . However, as one can see from Fig. 2, the densities of CMBR and synchrotron photons can be of the same order, as, for example in the knot A of the jet in 3C 273 . Therefore, the interactions of the highest energy $`\gamma `$-rays with the radio synchrotron radiation may well dominate, at least in the bright knots. Besides, the random B-field in the knots can be stronger than in the rest of the jet, which would result in increase of $`E_{\mathrm{crit}}`$ at which the electron trajectories are randomized. To demonstrate this effect we performed numerical calculations assuming very hard $`\gamma `$-ray spectrum with energy concentrated around $`10^{17}\mathrm{eV}`$. In Fig. 4 we show the evolution of rate of production of electrons in the knot, which determines the synchrotron luminosity of the ”shell”, together with the profile of the synchrotron X-radiation radiation measured by Chandra.

One can see that in the initial part of the jet, where synchrotron background is quite low, the rate of injection of the photoproduced electrons into the ”shell” is small. But in the knot A the rate increases rapidly. An observer who detects the jet from a side and does not see the ”main stream” of the cascade, may conclude that electrons are effectively accelerated in the knot.

Apparently, the formation of knots in the suggested model is a nonlinear process. The increase of synchrotron luminosity at some point of the jet would lead to an increase of the rate of ejection of electrons from the ”main stream”. This, in its turn leads to the further increase of synchrotron luminosity and formation of bright knots. Note that although the radio luminosity of the jet is directly related to the injection rate of electrons in the ”shell”, and therefore to the VHE $`\gamma `$-ray luminosity of the central source, it depends on some other factors as well, like the time of “operation” of the central source, the strength of the magnetic field in the jet, the escape of electrons from the jet, etc. While the X-ray data tell us about the VHE $`\gamma `$-ray luminosity of the central source at the present epoch, the radio data rather reflect the history of evolution of the the central source. An important question is whether photons with energies larger than $`10^{16}`$ eV can be effectively produced and at the same time freely escape the dense photon environments in the central engines of AGN which are believed to be powerful particle accelerators up to $`E_{\mathrm{max}}10^{19}`$ eV .

The accelerated protons can produce VHE $`\gamma `$-rays interacting with the ambient photon fields (supplied, for example by accretion disk around the massive black hole) through photo-meson process. Since we need a beamed $`\gamma `$-radiation emitted at a small angle to the MHD jet, the protons should cross the photon field almost rectilinearly. Therefore, the condition of a high proton-to-gamma conversion efficiency in the production region of a linear size $`R`$ implies $`\tau _{p\gamma }=\sigma _{p\gamma }Rn_{\mathrm{ph}}1`$. On the other hand, the produced $`\gamma `$-rays can effectively escape the production region if $`\tau _{\gamma \gamma }=\sigma _{\gamma \gamma }Rn_{\mathrm{ph}}<1`$. Thus, the production region can freely escape only those $`\gamma `$-rays for which $`\sigma _{\gamma \gamma }<\sigma _{p\gamma }`$. This is possible for very energetic $`\gamma `$-rays in a “hot” ambient photon gas. As an example, we show in Fig. 5 the spectra of protons and $`\gamma `$-rays emerging from a source filled with thermal radiation with $`T=10^4`$ K which serves a target for protons within a narrow energy interval between $`10^{18}`$ and $`10^{19}`$ eV. Since the center of mass energy of $`\gamma \gamma `$ interaction significantly exeeds the threshold of pair production, $`(E\mathrm{kT})/m_\mathrm{e}^2c^410^6E_{17}`$, at $`E10^{17}`$ eV the pair-production cross-section becomes less than the (almost energy-independent) photo-meson production cross-section, $`\sigma _{p\gamma }10^{28}\mathrm{cm}^2`$. Therefore the $`E10^{17}`$ eV $`\gamma `$-rays are not only effectively produced, but also are able to escape the source without significant losses. The numerical calculations shown in Fig. 5 illustrate this possibility. Note that a broader, e.g. power-law type proton spectrum would result in effective production of less energetic, $`E10^{16}\mathrm{eV}`$ $`\gamma `$-rays as well. However, due to the increase of the $`\gamma \gamma `$ cross-section, only a small fraction of these photons can escape the source.

To conclude, the electromagnetic cascade initiated by VHE photons interacting with ambient radiation fields in the large scale extragalactic jets is an attractive mechanism for production of ultra-relativistic electrons (with almost 100 percent efficiency) which can be responsible for the observed radio-to-X-ray spectra of jets. The trajectories of electrons with energies below $`E_{\mathrm{crit}}`$ are isotropized by the random magnetic field $`B_r`$. Such electrons form a ”shell” around the cascade. Observer, who looks at the jet from a side, sees synchrotron and inverse Compton radiation only form the ”shell” electrons. The cascade can be developed effectively in the jet provided that the strength of the random B-field does not exceed $`1\mu \mathrm{G}`$. On the other hand, when at a very large distances from the central source the random field is reduced to a very low level comparable with the intergalactic field, $`B10^9`$ eV or less, the cascade continues almost rectilinearly until the 10-100 TeV $`\gamma `$-rays start to interact effectively with the diffuse infrared background photons. This interactions would lead to the formation of observable giant pair halos with specific angular and energy distributions depending on the intensity of diffuse infrared background at the cosmological epoch corresponding to the redshift of the central source $`z`$ . And finally, if the central source is a blazar, i.e. the jet is pointed to the observer, we may expect beamed $`\gamma `$-ray emission with a characteristic for cascades $`E^{1.5}`$ type spectrum extending to 100 TeV. However, due to significant intergalactic absorption, the $`\gamma `$-rays will arrive with significantly distorted spectra. The possible implications of this mechanism for the TeV blazars like Mkn 421 and Mkn 501 are discussed elsewhere.

A.N. would like to thank SFB 375 der DFG for financial support.