text:
# Fluctuating Topological Defects in 2D Liquids: Heterogeneous Motion and Noise

\[

## Abstract

We measure the defect density as a function of time at different temperatures in simulations of a two dimensional system of interacting particles. Just above the solid to liquid transition temperature, the power spectrum of the defect fluctuations shows a $`1/f`$ signature, which crosses over to a white noise signature at higher temperatures. When $`1/f`$ noise is present, the 5-7 defects predominately form string like structures, and the particle trajectories show a 1D correlated motion that follows the defect strings. At higher temperatures this heterogeneous motion is lost. We demonstrate this heterogeneity both in systems interacting with a short ranged screened Coulomb interaction, as well as in systems with a long range logarithmic interaction between the particles.

\]

In liquids and glassy systems there has been considerable interest in dynamical heterogeneities which occur when certain regions of the sample have a higher mobility than the rest of the sample. Particle motion in these systems often occurs by means of correlated motion of a group of particles along a string like structure. . Most of the studies of heterogeneities in 2D and 3D have focused on systems that are supercooled near the glass transition in simulations  and experiments . A recent study considered a relatively simple system of a 2D dense monodisperse colloidal assembly which shows an ordered to disordered transition as a function of density . In this system the colloids form a defect free triangular lattice at high densities and disorder for lower densities. Near the disordering transition at the close packing density, collective heterogeneous motion of the particles appears, and the system consists of domains of sixfold coordinated particles, $`n_c=6`$, surrounded by grain boundaries composed of strings of $`n_c=7`$ and $`n_c=5`$ dislocations. The formation of grain boundaries or defect condensation in 2D monodisperse liquids has also been observed in colloidal assemblies  and in dusty plasmas .

Since the dislocations have a complex interaction with one another, consisting of a long range logarithmic repulsion and a short range attractive core interaction, one could expect complicated dislocation dynamics for systems with even a low density of dislocations. In the liquid state, creation and annihilation of dislocation pairs also occurs due to the thermally induced motions of the underlying particles. Although there have been numerous studies of the fluctuations of the particles in 2D liquids, the fluctuations in the defect density, defined as the number of particles with $`n_c6`$, have not been investigated. This measure would be easy to access in experiments on systems such as colloids and dusty plasmas where the individual particles can be directly imaged. It is not known how the density of $`n_c6`$ particles would fluctuate as the ordering transition is approached. If the defects are concentrated in clumps or grain boundaries, then it is likely that the creation or annihilation of defects will be highly correlated, which can give rise to $`1/f^\alpha `$ fluctuation spectra. Conversely, if the defects are appearing and disappearing in an uncorrelated manner, a white noise spectra would arise. We also wish to connect the formation of strings of defects with the appearance of correlated particle motion along 1D strings.

In this work we show that for 2D systems which form triangular lattices at low temperatures, for increasing temperature there is a transition from a non-defected regime to a defected regime where there is a proliferation of defects. Here we do not attempt to examine the nature of the disordering transition, such as whether there is an intermediate hexatic phase . Instead, we concentrate on the motions of the particles and the defect fluctuations in the liquid phase. Close to the disordering transition, in the defected regime, the system consists of regions of particles with sixfold order surrounded by strings or grain boundaries of 5-7 fold defects, and there are few free dislocations. We also find that in this region, collective particle motion occurs in 1D strings along the grain boundaries. When temperature is fixed, the density and structure of the defect strings fluctuates with time. Near the disordering transition, the defect density shows large fluctuations with $`1/f`$ spectra. For increasing temperatures, the defect density increases and the string like structures of the defects and the 1D string like motion of the particles are lost. At the same time, the low frequency noise power increases and the spectrum crosses over to a white noise signature, indicating the loss of correlations in the creation or annihilation of defects.

We numerically study 2D systems at finite temperature using Langevin dynamics. We consider monodisperse particles in a sample with periodic boundary conditions and have investigated two types of interactions. Most of the results presented here are for particles interacting via a Yukawa or screened Coulomb interaction potential, $`V(r_{ij})=(Q^2/|𝐫_i𝐫_j|)\mathrm{exp}(\kappa |𝐫_i𝐫_j|)`$. Here $`𝐫_{i(j)}`$ is the position of particle $`i(j)`$, $`Q`$ is the charge of the particles, $`1/\kappa `$ is the screening length, and we use

$`\kappa =2/a_0`$, where $`a_0`$ is the lattice constant. We choose to study this interaction since the quantities we study in this work can be accessed in systems such as 2D colloidal assemblies and dusty plasmas, in which the particles interact via a screened Coulomb interaction. We also consider particles interacting via a logarithmic potential, $`V(r)=\mathrm{ln}(r)`$. Unlike the Yukawa potential, this interaction is long-range so that an interaction cutoff cannot be implemented. For increased computational efficiency we use the summation method of Ref. . A physical example of logarithmically interacting particles with overdamped dynamics is vortices in thin-film type-II superconductors.

The equation of motion for particle $`i`$ is:

$$\frac{d𝐫_i}{dt}=\underset{ji}{\overset{N_c}{}}_iV(r_{ij})+𝐟_T$$
(1)
where $`𝐟_T`$ is a randomly fluctuating force due to thermal kicks with $`<𝐟^T(t)>=0`$ and $`<𝐟^T(t)𝐟^T(t^{})>=2k_BT\delta (tt^{})`$. In all the simulations we initially start from an ordered triangular configuration. We fix the temperature $`𝐟_T`$ for $`10^6`$ time steps before we begin to take data. We note that previous simulations have shown that extremely long time transients can arise near the order to disorder transition of up to $`2\times 10^6`$ time steps for system sizes of $`36864`$ particles . In this work we limit ourselves to system sizes of $`N_c=2000`$ or less, so that $`10^6`$ time steps is adequate for equilibration. In addition, histograms of the time series of the defect density are Gaussian, which is further evidence that our systems are equilibrated. We measure temperature in units of $`T/T_d`$ where $`T_d`$ is the temperature at which the first free dislocations appear.

In Fig. 1(a) we show the Voronoi construction of a system of particles with a screened Coulomb interaction for $`T/T_d=1.04`$. The Voronoi construction is similar to the Wigner-Seitz cell construction and an individual particle is located in the center of each cell. If a particle has six nearest neighbors, $`n_c=6`$, then the polygon has six sides. In Fig. 1 particles with $`n_c=6`$ are white, $`n_c=5`$ are dark gray, and $`n_c=5`$ are black. We analyze a series of such images at fixed $`T/T_d=1.04`$. We find that 29% of the particles are defected and that most (94.5%) of the defects are part of a cluster or condensate rather than free. In some regions, the defects form strings or grain boundaries. There are also a small number of free disclinations present, indicating that we are in the liquid phase rather than a hexatic phase, according to 2D continuous melting theories. Clustering of defects has been observed in experiments on dusty plasmas  which find comparable numbers of free dislocations. In Fig. 1(b) we show the same system at a higher temperature, $`T/T_d=9`$. Here, the number of defects is higher with most of the defects again predominately appearing in clusters; however, there is no string like characteristic to the clusters. Instead, the defects form clump like objects.

In Fig. 2(a) we show the positions of the particles (black dots) and trajectories (lines) for the same system in Fig. 1(a). The trajectories are taken for the particle motion over 10000 time steps. Fig. 2(a) shows that the motion is heterogeneous with certain regions in which the particles do not move. These areas are also the defect free regions. In the areas where there is motion, the particles move by $`a`$, and there is some tendency for the motion to occur along string like paths. We also often find places where the motion occurs in a circular path with an immobile particle in the center. A similar kind of motion in 1D like strings and rotations was seen in experiments on dense colloidal liquids . If we take trajectories for longer times, eventually all the particles take part in some motion. In Fig. 2(b) we show the positions and the trajectories for the high temperature phase in system Fig. 1(b) for $`2000`$ time steps. Here the motion occurs everywhere with no evidence for correlated stationary regions. We find that even for very short times,

there are no correlations in the motion.

We next consider a way to characterize the behaviors of the system at the different temperatures by measuring the fluctuations of the defect density as a function of time. We compute the defect configuration and density every 20 time steps for 20000 frames and obtain a time series of the defect density vs time for several temperatures from $`T/T_d=1.04`$ to $`10`$. We have considered a variety of sampling rates and find consistent results. In Fig. 3(a) we show a portion of the time series of the fraction of $`n_c=6`$ particles, $`P_6(t)`$, for $`T/T_d=1.04`$ (upper curve) and $`T/T_d=7.0`$ (lower curve). For $`T/T_d=1.04`$ the fluctuations show long time variations, while for the higher temperature the fluctuations are very rapid. In Fig. 3(b) we plot the power spectrum $`S(f)`$ of $`P_6(t)`$ for $`T/T_d=1.04`$, which fits well to a $`1/f^\alpha `$ scaling over more than three decades with the best fit $`\alpha =1.04`$, close to $`1/f`$ noise. As the temperature increases, the spectrum changes from $`1/f`$ to white noise ($`\alpha =0`$) for low frequencies. The small frequencies correspond to the long

time correlations in the system indicating that correlation times of the defect creation or annihilation decrease at higher temperatures. In Fig. 3(c), the power spectrum for $`T/T_d=7.0`$ is white with $`\alpha =0`$ for a large portion of the curve. We have also considered $`P_6(t)`$ for $`T/T_d<1.0`$ where there are a small number of bound dislocation pairs present. Here we find white noise with small amplitude fluctuations.

The magnitude of the noise power, obtained by integrating $`S(f)`$ over the first octave of frequencies, is related to the defect density and the disordering transition. In Fig. 4(a) we show $`P_6`$ as a function of $`T`$ and in Fig. 4(b) we show the corresponding noise power $`S_0`$. The peak in $`S_0`$ coincides with the onset of the defect proliferation. For increasing $`T`$, the noise power falls and saturates when the spectrum becomes white. There is a finite $`S_0`$ for $`T/T_d<1.0`$ since, as noted, pairs of bound dislocations can still be thermally created. We have repeated the same set of simulations for particles interacting with a long range logarithmic potential and again find a noise power peak coinciding with the defect proliferation

(insets of Fig. 4) and a $`1/f`$ power spectrum just above $`T_d`$.

We have also considered a sudden quench from $`T/T_d>1.0`$ to $`T/T_d<1.0`$. Here the system is out of equilibrium, and approaches equilibrium by means of the annihilation of the defects. In Fig. 5 we show the particle motions after a quench from $`T/T_d=1.04`$ to $`T/T_d=0.2`$. We find that the number of defects decreases rapidly to a saturation point where a few defects or grain boundaries remain. The motion of the particles corresponds to the annihilation of the defects. As seen in Fig. 5, the motion during the defect annihilation process has the same string like nature as the equilibrium $`T>T_d`$ motion near the defect proliferation transition. This result suggests that the heterogeneous particle motion is produced by the motion of defects, particularly the creation or annihilation of these defects in a correlated manner.

In summary, we have proposed a new way to examine dynamical heterogeneities in a liquid by measuring the fluctuations in the topological defect density. We find that near the onset of defect proliferation, the particle motion is heterogeneous and the defects cluster together into grain boundaries or strings. The defect fluctuations in this phase have a $`1/f`$ character and a large noise power. For higher temperatures the heterogeneities are lost and the fluctuation spectrum is white. We have also examined the defect annihilation after a quench from the disordered phase to the ordered phase and find that the particle motion also shows heterogeneities as the defects are annihilated. We have considered systems with screened Coulomb interactions and logarithmic interactions and find similar behavior. Our predictions can be easily tested in 2D dense colloidal liquids. It would also

be interesting to study the fluctuations of defects in 3D. This would be experimentally possible in 3D colloidal assemblies using confocal microscopy.

We acknowledge helpful discussions with S. Bhattacharya, D. Grier, V. Vinokur, E. Weeks, and M. Weissman. This work was supported by the US Department of Energy under Contract No. W-7405-ENG-36.


format: text

text:
# Interaction between dust grains near a conducting wall

## 1 Introduction

A good deal of experiments in dusty plasma physics are performed with aerosol grains suspended in the plasma sheath. A negatively charged dust grain in the sheath levitates above a horizontal electrode due to the balance of two forces: the gravity force directed downwards and the sheath electric field that pushes the grain upwards. The near-sonic or the supersonic ion flow in the sheath creates the wake field downstream of the grain. Since the latter is confined by a certain Mach cone, it is commonly accepted that the wake field affects the motion of grains, which are situated downstream only. The usual assumption is that the intergrain potential is smooth in the horizontal direction, i.e., two grains levitating at the same height repel one another via the screened Coulomb potential. The structure of the wake field and the grain interaction in an ion flow were studied in details \[1-5\]. The asymmetric interaction of vertically aligned grains was also observed experimentally .

Analytical theory and computer simulations cited above assumed that the plasma density is constant and the influence of the conducting electrode is negligible. Although both assumptions evidently fail under conditions of a real plasma sheath, taking into account the plasma non-uniformity seems to be a very difficult problem. In order to estimate the influence of the electrode upon grain interaction, here I use the zeroth approximation, which however, seems the only one treatable analytically.

In the present paper, the following simple model is accepted. Let there be the monoenergetic ion stream entering a conducting electrode (or a grid) located at the horizontal plane, $`z=0`$. The stream velocity, $`u`$, exceeds the ion thermal velocity but it is much less than the electron thermal velocity. Electrons are Boltzmann distributed. Two problems are addressed: first, how does the electrode modify the interaction between two grains levitating at the same height and, second, how this affects the spectrum of dust acoustic waves propagating along a single dust layer.

## 2 Intergrain interaction

The electrostatic potential produced by a point charge, $`Q`$, located at $`𝐫=𝐫_0`$ is given by the solution of the Poisson equation

$$\mathrm{\Delta }\widehat{\epsilon }\phi (𝐫)=4\pi Q\delta (𝐫𝐫_0),$$
(1)
where $`\widehat{\epsilon }`$ is the operator of the static dielectric permittivity of an ambient plasma. Within the accepted model, the spatial Fourier transform of $`\widehat{\epsilon }`$ is given by

$$\epsilon (𝐤)=1\frac{\omega _i^2}{(k_zu+i0)^2}+\frac{k_D^2}{k^2},$$
(2)
where $`u`$ is the velocity of the ion flow, $`\omega _i`$ is the ion plasma frequency and $`k_D`$ is the inverse electron Debye length. The ion flow is parallel to the $`z`$ axis and directed downwards, $`u<0`$.

In unbounded plasma, the natural boundary condition for Eq. (1) is $`\phi |_r\mathrm{}0`$. It is convenient to express the solution to Eq. (1) in terms of its Green function, i.e., $`\phi (𝐫)=QG^0(𝐫𝐫_0)`$, where the Fourier transform of $`G^0(𝐫)`$ with respect to the transverse coordinate $`𝐫_{}=(x,y)`$ is

$$G_k_{}^0(z)=\frac{dk_z}{2\pi }e^{ik_zz}\frac{4\pi }{k^2\epsilon (𝐤)}.$$
(3)
Here $`𝐤_{}=(k_x,k_y)`$ stands for the transverse components of a wave vector. With the dielectric function given by Eq. (2), the zeros of denominator in Eq. (3) are

$`k_z`$ $`=`$ $`\pm k_Dqi0,`$
$`k_z`$ $`=`$ $`\pm ik_D\kappa ,`$ (4)
where

$`q^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left\{\sqrt{(q_{}^2+1\mu ^2)^2+4\mu ^2q_{}^2}q_{}^2+\mu ^21\right\},`$ (5)
$`\kappa ^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left\{\sqrt{(q_{}^2+1\mu ^2)^2+4\mu ^2q_{}^2}+q_{}^2\mu ^2+1\right\},`$ (6)
$`q_{}=k_{}/k_D`$ is the normalized transverse wave vector, and $`\mu =\omega _i/k_Du`$ is the inverse Mach number.

The integral in Eq. (3) is readily evaluated resulting in

$$G_k_{}^0(z)=\frac{2\pi }{k_D}\frac{1}{q^2+\kappa ^2}\left[\kappa e^{\kappa k_D|z|}+2\theta (z)q\mathrm{sin}\left(qk_Dz\right)\right].$$
(7)
The spatial structure of the potential is recovered by the Fourier transform of the Green function (7) with respect to $`k_{}`$. The first term in parenthesis in expression (7) gives rise to the Debye-Hückel potential distorted by the ion flow, while the second term represents the wake field situated downstream of the charge \[1-3\]. The opening of the Mach cone confining the wake and the field structure inside it depend on the stream velocity.

Now we turn to the evaluation of the electric potential of a charge located near a conducting wall. Let the wall be situated at $`z=0`$ plane, while the charge is placed above it, $`z_0>0`$. Then the potential is given by the solution of Eq.(1) supplemented with the boundary condition $`\phi |_{z=0}=0`$. As it is well-known from electrostatics, we can make allowance for this boundary condition by introducing the surface charge density, $`\sigma _k_{}^{ind}`$, induced at the conducting surface. Then, the potential of a unit charge is written as

$$G_k_{}(z,z_0)=G_k_{}^0(zz_0)+G_k_{}^0(z)\sigma _k_{}^{ind}.$$
(8)
Taking into account the boundary condition, $`G_k_{}(0,z_0)=0`$, we find the surface charge density $`\sigma _k_{}^{ind}`$ and, finally,

$$G_k_{}(z,z_0)=G_k_{}^0(zz_0)\frac{G_k_{}^0(z)G_k_{}^0(z_0)}{G_k_{}^0(0)}.$$
(9)
One may doubt whether description a bounded dispersive medium in terms of the response function of an unbounded medium is justifiable. However, more accurate and lengthy calculations give the same result. The physical reason is that there are no ions reflected by the conducting wall within the present model. The mathematical reason is that Eq. (1) actually masks a set of partial differential equations with two real characteristics directed downwards.

Of particular interest for the following is the interaction potential of two charges placed at the same height, $`z_0`$. Since

$$G_k_{}(z_0,z_0)=\frac{4\pi }{k_D}\frac{e^{\kappa a}}{\kappa ^2+q^2}\left[\kappa \mathrm{sinh}(\kappa a)+q\mathrm{sin}(qa)\right],$$
(10)
where $`a=k_Dz_0`$, the normalized potential is

$$w(\rho )=2\underset{0}{\overset{\mathrm{}}{}}q_{}𝑑q_{}J_0(q_{}\rho )\frac{e^{\kappa a}}{\kappa ^2+q^2}\left[\kappa \mathrm{sinh}(\kappa a)+q\mathrm{sin}(qa)\right],$$
(11)
where $`G(r_{},z_0,z_0)=k_dw(\rho )`$ and $`\rho =k_dr`$. The asymptotic behaviour of $`w(\rho )`$ is determined by Eq. (10) at $`q_{}0`$. The latter depends essentially on whether the the ion flow is supersonic or subsonic. In the case of the supersonic flow, $`\mu <1`$, the roots (5,6) at $`q_{}0`$ are approximated as

$$qq_{}\frac{\mu }{\sqrt{1\mu ^2}},\kappa \sqrt{1\mu ^2},$$
(12)
and the leading term of the asymptotic expansion of $`w(\rho )`$ (11) is

$$w(\rho )|_\rho \mathrm{}\frac{2}{\mu }(1\mu ^2)^{3/2}e^{a\sqrt{1\mu ^2}}\frac{\theta (a\mu \rho \sqrt{1\mu ^2})}{\sqrt{a^2\mu ^2\rho ^2(1\mu ^2)}}+O\left(e^r\right).$$
(13)
This expression may be interpreted as a mirror reflection of the wake field produced by the grain at $`0<z<z_0`$. More detailed numerical investigation of the potential (13) shows that $`w(\rho )`$ is always positive for $`\mu <1`$.

Quite another behaviour is observed for the case of the subsonic flow, $`\mu >1`$. The roots (5,6) at $`q_{}0`$ are now

$$q\sqrt{\mu ^21},\kappa q_{}\frac{\mu }{\sqrt{\mu ^21}},$$
(14)
while the potential at infinity behaves like

$$w(\rho )|_\rho \mathrm{}\frac{2a\mu }{\mu ^21}\mathrm{sin}\left(a\sqrt{\mu ^21}\right)\frac{1}{\rho ^3}.$$
(15)
The most important distinction between the expressions (13) and (15) is that in the subsonic regime the potential is attractive if

$$\mathrm{sin}\left(a\sqrt{\mu ^21}\right)<0.$$
(16)
The numerically evaluated example of the potential (15) demonstrating the attraction is depicted in Fig. 1. It should be pointed out that the inequality (16) guarantees the long-scale attraction between grains. With the opposite inequality imposed on $`a`$ and $`\mu `$, the potential is repulsive at large distance but the attractive branch may appear at smaller scales, as shown in Fig. (2).

Evidently, the existence of the attractive branch of the interparticle interaction may result in formation of various patterns and clusters even in the absence of the external potential well confining grains in the horizontal direction. Also, the even distribution of grains in a dust layer may become unstable. The latter possibility is discussed in the next section.

## 3 Dust layer

Now consider a two-dimensional gas consisting of dust grains hovering over the conducting electrode. Ignoring intergrain correlations, the linearized equations of motion are written as

$`{\displaystyle \frac{\sigma }{t}}`$ $`+`$ $`_{}(\sigma _0𝐯)=0,`$ (17)
$`{\displaystyle \frac{𝐯}{t}}`$ $`=`$ $`\nu 𝐯{\displaystyle \frac{Q^2}{M}}k_d_{}{\displaystyle 𝑑𝐫_{}^{}w(k_D(𝐫_{}𝐫_{}^{}))},`$ (18)
where $`\sigma _0`$ is the unperturbed value of the surface density and $`\sigma `$ is the density perturbation. Here I consider horizontal motions only, i.e., $`v_z=0`$. The term $`\nu 𝐯`$ in Eq. (18) corresponds to the grain friction on an ambient neutral gas. The intergrain interaction is described by the potential given by Eqs. (11). Although the grain charge, $`Q`$, generally depends on the ambient plasma parameters, for simplicity we ignore its variability.

Assuming that all quantities are proportional to $`\mathrm{exp}(i\omega t+i𝐤_{}𝐫_{})`$ we easily get the dispersion relation for the gas oscillations:

$$\omega (\omega +i\nu )=g_d\frac{k_{}^2}{2\pi }G_k_{}(z_0,z_0),$$
(19)
where $`g_d=2\pi Q^2\sigma _0/M`$.

This expression describes dust sound waves in the continuous medium approximation. In the long-wave limit, $`q_{}0`$, the dispersion relation is of the form

$$\omega (\omega +i\nu )=\frac{g_dk_dq_{}^2}{\sqrt{|1\mu ^2|}}\{\begin{array}{cc}e^{a\sqrt{1\mu ^2}}\mathrm{sinh}(a\sqrt{1\mu ^2}),\hfill & \mu <1\hfill \\ \mathrm{sin}(a\sqrt{\mu ^21}),\hfill & \mu >1\hfill \end{array}$$
(20)
Evidently, the layer is unstable, i.e., $`\text{Im }\omega >0`$, if $`G_k_{}(z_0,z_0)<0`$. In the long-wave limit this is possible in the subsonic flow only and the corresponding constraint coincides with Eq. (16). More detailed investigation shows that the potential, $`G_{}(z_0,z_0)`$, is always positive if $`\mu <1`$. However, in the subsonic regime, $`\mu >1`$, there are regions of instability, that is, $`G_k_{}(z_0,z_0)<0`$ if $`0<q_{}<q_{max}(a,\mu )`$. The latter are shown as shadowed areas in Fig. 3. With the increasing distance ($`a\mathrm{}`$) to the wall or decreasing stream velocity ($`\mu \mathrm{}`$) the instability regions shrink to zero, $`q_{max}(a,\mu )a^{1/2},\mu ^{1/2}`$.

## 4 Conclusion

To summarize, we have shown that the presence of the conducting wall may drastically change the electrostatic interaction of the dust grains in an ion flow. In particular, the electrostatic image of the grain wake field may result in attraction between grains levitating at the same height that, in its turn, yields the Jeans-type instability of the dust layer.

It would be unduly naive to draw quantitative conclusions from the present calculations. However, it seems reasonable that even in a real plasma sheath, which is essentially non-uniform, the electrostatic image of the grain wake field may also affect the motion of another grains outside the Mach cone. Although in most experiments the screened Coulomb interaction is observed, there are indications that the intergrain potential may be more complicated. It was recently reported that, under certain conditions, a void, i.e., a dust-free region, appears in a central part of a single dust layer . The emergence of a two-dimensional void in a layer consisting of some hundreds of grains can hardly be explained in the same manner as a three-dimensional void; the latter requires strong influence of the dust component upon the discharge structure . Although currently one cannot exclude that some additional external forces appeared in the experiment , we can conjecture that the 2D void formation is provided by complicated intergrain interaction, for example, the one described in this paper.

## Acknowledgements

This study was supported in part by the Russian Foundation for Basic Research (project no. 02-02-16439) and the Netherlands Organization for Scientific Reseach (grant no. NWO 047.008.013).


text:
# Electrons in Dry DNA from Density Functional Calculations

## I Introduction

The idea of using DNA chains as molecular wires in nano-devices has been proposed and is being explored in the nano-science community Mirkin96 ; Alivisatos96 ; Braun98 , given the extremely high selectivity of a DNA strand to pair with its complementary partner, which facilitates the design of complex circuitry at the nanometer scale. Even though some researchers have chosen to coat Braun98  or modify Rakitin01  the DNA wires in order to improve its conduction properties, the original idea was to use pure DNA, assuming that electric current could flow along the wires. This assumption has proven to be less trivial than expected, however, and has originated an important controversy.

Direct measurements show from low resistance DC conduction Fink99  (even superconduction Kasumov01 ), to no conduction at all Braun98 ; dePablo00 ; Gomez02 , including results with non-ohmic characteristics Porath00 ; Yoo01  or appreciable ac conductance Gruner00 . The transport has been proposed to be band-structure-like Porath00 , by tunneling Beratan97 , solitonic BenJacob98 , polaronic Schuster99 ; Conwell00 ; Yoo01 , or facilitated by fluctuations Ladik99 ; Bruinsma00 . This rich scenario is not surprising when considering the comparably small amount of clean-cut information available (mostly theoretical) about basic aspects of the electronic structure of dry DNA Ladik71 ; Ladik83 ; Pablo97 ; Yang98 ; Ye00 , (in all the experiments mentioned above the samples were dried in one way or another). In addition to the fundamental complexity of the system (large system size, weak interactions, structure floppiness), there is a substantial lack of detailed knowledge about the environment the DNA molecules are left with after the drying process. Most likely, residual water molecules will stay adhered to the chain, together with the cations counterbalancing the negative charge of the molecule, all quite uncontrolled.

A theoretical simulation benefits from absolute control at the atomic scale. In this particular case, the system can be prepared in the best possible way for conduction. Previous simulations Ladik71 ; Ladik83 ; Pablo97 ; Yang98 ; Ye00  have offered very insightful information for related problems, in spite of being based on relatively drastic approximations in the electronic structure, due to the system size. They were mostly aimed at the description of DNA properties in live conditions, as were the recent first-principles studies in wet DNA Landman01 ; Zhang02 ; Parrinello02 .

Recent developments in linear-scaling first-principles Soler02 ; Bowler02 ; Scuseria00  methods based on density-functional theory Payne92  (DFT) allow the study of an infinite DNA double helix from first principles. DFT functionals based on the Generalized Gradient Approximation (GGA) have shown to provide very satisfactory accuracy in molecules akin to DNA Sponer96 ; Hutter96 . In this work, an extensive study is presented on a set of 30 nucleic base pairs that validates the theoretical method. Then, an infinite DNA chain of trivial sequence is studied: all guanines (G) on one strand, all cytosines (C) on the other, in the Watson-Crick pairing arrangement. We simulate it completely dry and in vacuo, and in an A-helix structure, the most relevant one for these conditions. The electronic structure is analyzed after relaxing the structure with the first-principles forces. Genetic sequence change effects are explored on a similar system where one guanine-cytosine pair is swapped every eleven base pairs, then relaxed and analyzed. Preliminary results of this study were published in a joint experiment-theory paper dePablo00 .

## II Method

The first-principles DFT calculations reported here were performed using the Siesta method Soler02 ; Pablo96 ; Daniel97 , which has already been applied to a large variety of systems websiesta , including biomolecules dePablo00 ; Marivi00 . It is a fully linear-scaling method, both in the building of the DFT Hamiltonian and in its solution. GGA PBE  was used for electron exchange and correlation. Core electrons were replaced by norm-conserving pseudopotentials TroullierMartins91  in their fully non-local formulation KB . A uniform mesh with a planewave cutoff of 150 Ry was used for integrations in real space Soler02 . The basis sets used are made of numerical atomic orbitals of finite range and are described below.

The unit cell of the A-helix studied here includes eleven poly(G)-poly(C) (pGpC) base-pairs with a total of 715 atoms. Periodic boundary conditions are used to simulate infinitely long DNA chains well separated from each other, with negligible mutual interactions. Only the $`\mathrm{\Gamma }`$ $`k`$-point of the 715-atom cell was included in the calculations, which corresponds to 11 $`k`$-points, if the helical symmetry is considered. The initial geometry was that of an experimental low-resolution A-helix structure for the pGpC sequence Landoldt-Bornstein , The geometry was relaxed by means of ab initio linear-scaling DFT, as described elsewhere ADNA-PRB . The final structure had still well defined A-helix characteristics ADNA-PRB ; Calladine97 . A second DNA realization (swapped hereafter) was built by swapping one single G-C pair into C-G out of the eleven in the cell, followed by another full relaxation.

The electronic structure is analyzed in terms of the Kohn-Sham eigenvalues and eigenfunctions, based on Janak’s theorem Janak78 . They are obtained by a single cube-scaling diagonalization, using the final geometry and electron density obtained with the linear-scaling method. The known band-gap problem of DFT is not important here, since the particular value of the band gap is of no relevance for the conclusions of this study as long as there is a gap, which is the case. The band widths have been shown to be well described by this method, at least up to the kind of accuracy needed in this work. A much more costly calculation including explicit electronic correlation would not be justified. Mott-like correlations, that could seem to be relevant for the narrow band-widths characteristic of the problem (see below), are not important given the small density of carriers (electrons or holes, see discussion below) expected in the system. The effect of disorder and vibrations is explicitly considered below. The next section is devoted to a study of nucleic base pairs where the focus is mainly on the hydrogen bonds, which might in principle be worst described by our methodology. It will be followed by the results for A-DNA.

## III Study on base pairs

A systematic study of 30 nucleic base pairs has been performed with the methodology presented above, in order to assess the reliability of the approximations, but also obtaining new results. For this study we have used a standard doble-$`\zeta `$ polarized (DZP) basis set, namely, a double (split valence) basis for each valence orbital plus polarization functions in all the atoms. The cutoff radii for the atomic orbitals of each element were obtained for an energy shift Emilio99  of 50 meV.

A significant range of configurations of the four bases guanine, cytosine, adenine and thymine (G, C, A, T) are considered, the same as those studied by $`\stackrel{ˇ}{\mathrm{S}}\mathrm{poner}`$ et al. Sponer96  in their MP2 study. The Watson-Crick configurations are designated WC, and the Hoogsteen, reversed Hoogsteen and reversed Watson-Crick appear as H, RH, RWC respectively. Other configurations are distinguished simply with numbers, eg. AA1, AA2, etc. following the nomenclature of Hobza and Sandorfy Hobza . The structures of the bases and base-pairs studied in this work can be found in Figures 1 and 2 of Ref. Sponer96, . For the numbering of the atoms we followed Ref. Sponer2, . In the following our results are compared with those of other methods. The direct comparison with experimental geometries has been shown to be misleading because of the important charge rearrangements that happen in the crystal phases used in experiments Celia00 .

The interaction energy $`E_{int}`$ in the following is defined as the energy of the base-pair minus the energy of each base with the same geometry it has in the pair. The total stabilization energy $`E_t`$ is defined as the difference between the energy of the pair and that of each base in its isolated optimal geometry. The difference between both is thus the deformation energy, i.e., the increase in intramolecular energy due to the geometry change when the base-pair is formed.

In this study all the energies have been corrected for basis set superposition error (BSSE) using the standard Boys-Bernardi counterpoise correction Boys . The correction found for $`E_{int}`$ is used for $`E_t`$ as well, given the fact that the relaxation of the isolated base molecules has been found to change the BSSE correction by less than 10%. No BSSE correction was included in the forces. The effect of this approximation was gauged for the AA1 base-pair, for which BSSE corrected and uncorrected forces were calculated for independent relaxations, showing a deviation of 0.02 Å in the final acceptor-hydrogen distances.

Our results for $`E_{int}`$ as compared with the MP2 results of $`\stackrel{ˇ}{\mathrm{S}}\mathrm{poner}`$ and coworkersSponer96  are presented in table 1, as calculated for the same geometries obtained at the Hartree-Fock (HF) 6-31G\** level Sponer96 . For all the base-pairs except GG4, the agreement is quite satisfactory, with differences smaller than 8% and much less in most cases. GG4 is an exception to the general trend with a relative difference of 26% that is reduced to 16% by extending the basis range to the longer radii given by an energy shift of 10 meV. A similar correction for the other base pairs gives much smaller variations. The standard deviation of our results compared to the MP2 values is of 0.73 kcal/mol, the same value found by $`\stackrel{ˇ}{\mathrm{S}}\mathrm{poner}`$ et al.Sponer96 using the B3LYP functional b3lyp  on the same geometry. The standard deviation between both sets of DFT results (our PBE and their B3LYP) is of 0.85 kcal/mol.

Concerning geometries, the only pair relaxed with MP2 is cytosine-cytosine Sponer96 . Table 2 compares the length of the hydrogen bond and both $`E_{int}`$ and $`E_t`$ for various methods, showing again a very satisfactory agreement.

The results shown in this section validate the approximations used in our method for the system under study, including the sensitive hydrogen bonds. Even if the piling up of bases was not investigated, other studies show Hutter96  that GGA functionals can describe them with the required accuracy when connected by the sugar-phosphate chains.

The whole set of base pairs as well as the isolated bases have been relaxed with our method. The results for geometry and interaction energies can be found in Maider00 . We extract in table 3 the results for the guanine-cytosine pair, of relevance for the DNA calculations in the next section. The results of several authors Sponer96 ; Celia00 ; Santamaria  are compared, which include different GGA functionals, namely B3LYP b3lyp , BP86 p86 ; b86 , and VWN-BP p86 ; b88 ; vwn . The results of the DZP basis of this study are also compared with those of the basis set used for DNA below, DZ(P), a double-$`\zeta `$ basis Emilio99 ; radii , which was polarized for phosphorous and for the atoms involved in the hydrogen bonds. The Kohn-Sham eigenvalues around the Fermi level for an isolated GC pair were checked for DZ(P), which compare well with the corresponding values for the DZP basis. The convergence of the spatial range of the basis orbitals, important for band widths in the DNA molecule, was also tested radii .

## IV Results for A-DNA

Figure 1 shows the bands of pGpC close to the Fermi level. As expected for a well saturated molecule, the Fermi level is in a clear band gap, which is of 2.0 eV in this case (this number is to be taken with caution since usual DFT functionals tend to underestimate the band gaps of insulators by a substantial amount). Electronic transport will thus be mediated by carriers in either the top of the valence band (holes, or radical cations in biochemical language) or in the bottom of the conduction band (electrons or radical anions), introduced by doping.

In the absence of disorder, the electronic structure close to the Fermi level shows well defined bands with eleven states per unit cell, one per base pair in the unit cell. The top-most valence band has a very small bandwidth of 40 meV (see figure 1), the energy separation with the next band below being ten times larger. This first band is associated to the $`\pi `$-like highest occupied molecular orbital (HOMO) of the guanines. The continuous line in the figure is the tight-binding band of a one-dimensional system with one orbital per unit cell. The excellent agreement of this extremely simple model with the ab initio points \[see figure 1 (b)\] demonstrates that each base pair contributes to this band with one single orbital that interacts negligibly with other orbitals in the pair. This allows the modelling of the basic physics of transport in this molecule with one single orbital per base pair with an electron hopping interaction of 10 meV. Other terms for a model Hamiltonian are quantitatively justified below.

Figure 2 displays the charge density of the states associated to this band, which appears almost exclusively on the guanines, with weight neither in the backbones nor in the cytosines. The lowest conduction band has a width of 270 meV, similar to the width of the spectroscopic feature observed by Porath et al. Porath00 . It is separated from the next band by 0.7 eV. Similarly to the HOMO situation, this band is made of the lowest unoccupied molecular orbital (LUMO) of the cytosines, as shown in figure 2, although the overlap in this case is more pronounced, as corresponds to an eight times larger bandwidth.

Note that any matrix element between HOMO and LUMO states is very much depressed by the very small spatial overlap between them. Optical absorption will thus be very weak for photon energies corresponding to the HOMO-LUMO gap (adequately corrected including correlation and excitonic effects). The unoccupied band above the LUMO band remains of a marked cytosine character. The third unoccupied band 3.9 eV above the HOMO band (see figure 1) is the first one with substantial weight on the guanines, and thus with an appreciable oscillator strength for absorption. The comparison of these results with those obtained for wet conditions Parrinello02  shows that the LUMO is quite sensitive to the environment, moving from the cytosines to the cations when in presence of Na<sup>+</sup> and water.

The effect of the swap on the electronic structure of the chain is dramatic, as can be seen in figure 3. The HOMO of the swapped guanine (indicated by an asterisk in figure 3) sinks 0.55 eV (14 times the HOMO bandwith) into lower valence band levels. This stabilization is of electrostatic origin ADNA-PRB . Figure 4 (left) shows the electronic density associated to the eleven highest occupied states for the swapped structure, showing the cut in the HOMO-state channel produced by the swapped pair. The situation is similar for the unoccupied band, albeit less dramatic, with the LUMO of the swapped cytosine raising 0.5 eV above the LUMO band.

All the results presented so far refer to fixed, relaxed geometries. It has been proposed, however, that the electron transfer could be mediated by phonons Schuster99 ; Conwell00 ; Ladik99 ; Bruinsma00 . Indeed, there are many modes in DNA that are very soft and related to the nucleic bases, since these flat, rigid molecules hang from the backbone by a single bond. It is beyond the scope of this study to calculate phonon modes and their interactions with electrons, but interesting qualitative information on the matter can be inferred from the conjugate-gradient (CG) minimizations. If close to the minimum, any snapshot during the minimization gives the atomic positions that correspond to the freezing of vibration modes with various amplitudes. The closer to the minimum, the more important the low-frequency modes versus the high-frequency ones valley . We study thus the electronic structure of geometries chosen from the CG path with an energy of around 0.9 eV above the minimum, that would correspond to room temperature or lower temperature . The highest occupied energy levels of such a geometry (called ragged) is presented in figure 3. If not as dramatic as the swapping, the effect of these “vibrations” on the electrons is also important, since the guanine HOMO one-particle states spread over an energy width eight times the original band width.

## V Discussion and Conclusions

Electronic transport in DNA requires the presence of carriers. In controlled biochemical experiments (in solution), doping is induced by intercalating or attaching to the DNA specific molecules of known tendency to donate or accept electrons Barton99 . In the dry experiments introduced earlier there is uncontrolled doping of DNA. Since the negative charge of each phosphate group along both chains has to be compensated by a cation from solution, in the drying process cations precipitate along the DNA helix decorating it. In addition to a very probable disorder in the cation positions and variability in the amount of water molecules hydrating them, there will also be deficiency or excess of counter-cations. Missing cations introduce holes in the valence band and cations in excess (or cations with higher charge) introduce electrons in the conduction band. The drying procedures of the dry experiments do not allow estimating the concentration of carriers, not even their sign, even if (in view of the ionization potential of G and the electron affinity of C) the presence of holes has been argued to be more likely. It is reasonable to expect in any case that carrier concentrations will remain low, and that for the more drastic drying procedures the likelihood is higher for the presence of holes, since cations will be pulled away with the water. pGpC has been reported to be $`p`$-type in at least one experimental setup Yoo01 .

It has been shown above that the HOMO band width is one order of magnitude smaller than the energy scale associated to sequence disorder. In terms of a one-dimensional Anderson model, this leads to electronic localization over very few base pairs. The localization is observed when plotting the separated eigenstates. Figure 4 (right) shows one electronic state within the HOMO energy range. The isosurface plotted is for a density 50 times smaller that that on the left, showing clearly the localization of the state within a few base pairs. Sequence disorder will therefore dramatically depress the possible band-structure-like DC conduction in DNA of arbitrary sequence ($`\lambda `$-DNA). The possibility Ratner00  of electron hopping or tunneling, like in proteins Beratan97  is open, but it is not expected to show the conductances close to one quantum reported for $`\lambda `$-DNA by some groups Fink99 ; Kasumov01 . For electron doping, disorder fluctuations are still comparable or larger than the LUMO band width, and the localization of the electronic states is still considerable, albeit not as extreme as for the hole-doped case.

It is known Phillips91  that certain types of disorder do not localize electrons as in the Anderson model. They represent, however, very particular realizations that cannot be expected from $`\lambda `$-DNA in general. Very interesting recent results Carpena02  show possible electron delocalization in 1D systems with more general long-range correlated disorder, as the disorder of natural DNA. More research is needed to see how robust is this result when including other effects like off-diagonal disorder or vibrations. In any case, the very large disorder-to-bandwidth ratio presented above makes extended states for holes very unlikely, the typical state being localized as displayed in figure 4 (right).

The particular experiment of Porath et al. Porath00  used a controlled DNA sequence, namely pGpC, for which the sequence disorder arguments do not apply. They explain their results in terms of band-like conduction. The results and analysis of Berlin and coworkers Ratner00  propose band-like conduction in this particular system as well. It is tempting to support this band-like condution assuming it proceeds via electrons in the LUMO, given the agreement of the widths of our LUMO band and the feature they observe spectroscopically (the sensitivity of the LUMO to the environment mentioned above has to be kept in mind, however). Such an explanation would be by no means conclusive in the light of our “vibration” results, since a deterioration of conduction should be expected with temperature, which does not seem to happen. The spread of HOMO states with frozen-in vibrations indicates the energy scale of electron-phonon interactions, again very substantial as compared with the band widths. Even if a more rigorous study is needed to quantify these interactions and their effect, our rough estimations above allow us to expect a substantial reduction of conduction with temperature even below room temperature.

Considering the floppy modes in DNA and their effect on the electron states, and the fact that a hole on a guanine represents a positive charge amidst negatively charged guanines, the idea of self trapped (small) polarons seems not at all unlikely. Polarons have already been suggested Schuster99 ; Conwell00  and supported Yoo01  in the literature. Some of the models used have assumed the polarons to be of an origin similar to that found in polyacetylene, i.e. of off-diagonal (hopping) origin. Polarons of electrostatic origin are described by an electronic Hamiltonian in which the diagonal matrix elements are the ones that change with inter-base separation and stabilize the polaron, as in Holstein’s model Holstein59 . A purely classical electrostatic model based on the ab initio results for the geometry and ground-state electrostatic characteristics of the system ADNA-PRB  displays polaronic charge excitations with a spread of a few base pairs. DFT calculations on the matter will be presented elsewhere Simone .

Acknowledgments. We are grateful to J. $`\stackrel{ˇ}{\mathrm{S}}\mathrm{poner}`$ and R. Santamaría for making their coordinates of bases and base-pairs available to us. We are also indebted to R. Weht for many useful discussions and his help during the first stages of this work. This work has been supported by the Spanish Ministerio de Ciencia y Tecnología under grant BFM2000-1312, and by the Fundación Ramón Areces. DSP acknowledges support from Spain’s MCyT under the Ramón y Cajal program. Calculations were performed at the CCCFC of the Universidad Autónoma de Madrid and at the PSMN of the École Normale de Lyon (SESC).


text:
# Project SANC (former CalcPHEP): Support of Analytic and Numeric calculations for experiments at Colliders

## 1 SANC PROJECT AND ITS ROOTS

The main goal of this project is the creation of a software product, accessible via Internet, for the automatic computation of pseudo- and realistic observables with a one-loop precision for various processes of elementary particle interactions, such as: $`12`$, $`13`$, $`22`$, $`14`$, $`23`$, etc. .

It has two roots: 1) Codes aimed at a theoretical support of HEP experiments, such as TOPAZ0 , ZFITTER  and similar ones.
2) Numerous FORM2-codes, written by the authors of ref.  while they were working on it.

It is supposed that the main software products of the project should be the Monte-Carlo event generators that are being created in collaboration with S. Jadach and Z. Was from INF (Krakow, Poland) and B.F.L. Ward from University of Tennessee (Knoxville, USA). On top of these generators, for some processes the semi-analytic codes are also provided, the latter being used both for the cross-checks of MC generators and for fits of inclusive observables to the theory predictions.

## 2 BASIC NOTIONS

### 2.1 The SM Lagrangian in the $`R_\xi `$ gauge

All the calculations start from the SM Lagrangian in the $`R_\xi `$ gauge. It depends on 25 input parameters (IPS), fields, and on three gauge parameters:

$$=(\text{IPS: 25 parameters, fields,}\xi __A,\xi __Z,\xi ),$$
(1)
where the fields are: fermions, vector bosons, physical Higgs field, $`H`$, and unphysical fields, $`\varphi ^0,\varphi ^\pm ,Y^A,Y^Z,X^\pm `$ . We give the example of Feynman Rules for vector boson propagators:

$$\begin{array}{cc}A:\hfill & \frac{1}{p^2}\left\{\delta _{\mu \nu }+\left(\xi __A^21\right)\frac{p_\mu p_\nu }{p^2}\right\},\hfill \\ Z:\hfill & \frac{1}{p^2+M__Z^2}\left\{\delta _{\mu \nu }+\left(\xi __Z^21\right)\frac{p_\mu p_\nu }{p^2+\xi __Z^2M__Z^2}\right\},\hfill \\ W^\pm :\hfill & \frac{1}{p^2+M__W^2}\left\{\delta _{\mu \nu }+\left(\xi ^21\right)\frac{p_\mu p_\nu }{p^2+\xi ^2M__W^2}\right\}.\hfill \end{array}$$

### 2.2 Scalar functions and reduction

At present, SANC knows how to deal with up to third-rank tensorial reduction to the usual scalar functions: $`A_0,B_0,C_0`$ and $`D_0`$; and to the auxiliary scalar functions: $`a_0,b_0,c_0`$ and $`d_0`$, which are due to particular form of photonic propagator in $`R_\xi `$ gauge. A new fortran library for numerical calculations of all these functions but one ($`D_0`$) is created and thoroughly tested by means of comparison with the other codes. For the $`D_0`$ function, we use an old coding of TOPAZ0 .

### 2.3 Amplitude basis, scalar form factors, helicity amplitudes

SANC computes the one-loop covariant amplitude ($`𝒜`$) of a process parametrized in a certain basis by a certain number of scalar form factors (SFF). Next, it computes helicity amplitudes (HA) in terms of SFF. As an example, we present three covariant amplitudes of the decays $`B(Q)f(p_1)\overline{f}(p_2)`$, where $`B=H,Z,W`$:

$`Hf\overline{f}:𝒜I_S,`$ (2)
$`\text{1 structure (S-basis), 1 SFF, 1 HA};`$
$`Zf\overline{f}:𝒜i\gamma _\mu \gamma _6_𝐋+i\gamma _\mu _𝐐+m_fD_\mu _𝐃,`$
$`\text{3 structures (LQD-basis), 3 SFF, 3 HA};`$
$`Wu\overline{d}:𝒜i\gamma _\mu \gamma _6_𝐋+i\gamma _\mu \gamma _7_𝐑`$
$`+m_uD_\mu \gamma _6_{\mathrm{𝐋𝐃}}+m_dD_\mu \gamma _7_{\mathrm{𝐑𝐃}},`$
$`\text{4 structures (LRD-basis), 4 SFF, 4 HA},`$
where $`D_\mu =\left(p_1p_2\right)_\mu `$. Besides SFF, the HA depend on kinematical factors and coupling constants. For instance, for $`Z`$ decay three HA are:

$`𝐀_{\mathrm{𝟎}\pm \pm }^𝐙`$ $`=`$ $`{\displaystyle \frac{gm_f}{c__W}}\left[a_f_𝐋+\delta _f_𝐐+{\displaystyle \frac{1}{2}}a_f\beta _f^2M__Z^2_𝐃\right],`$
$`𝐀_{\pm \pm }^𝐙`$ $`=`$ $`{\displaystyle \frac{gM__Z}{\sqrt{2}c__W}}\left[a_f(1\beta _f)_𝐋+\delta _f_𝐐\right],`$ (3)
where $`\beta _f^2=14{\displaystyle \frac{m_f^2}{M__Z^2}},`$
and $`\delta _f=v_fa_f=2Q_fs__W^2,a_f=I_f^{(3)}.`$
We note that the number of SFF and the number of independent non-zero HA always coincide. In general, a $`2f2f`$ process with four different external fermion masses is described by 16 SFF and 16 independent HA. In case of $`e^+e^{}t\overline{t}`$ process, one has six SFF and six independent HA if the electron mass is ignored .

## 3 STATUS OF THE PROJECT

### 3.1 Four levels of computer system

Here we explain what the SANC system does at its four levels in the case of calculation of pseudo-observables (decay rates, event distributions) for the simplest decays: $`H(Z,W)f_1\overline{f}_2`$.

Level 1: the chain of calculations from $`_{SM}`$ to the ultraviolet-free helicity amplitudes is realized. (All codes are written in FORM3 ). It contains four sublevels, which calculate:

1. the scalar form factors;
2. the soft and hard photonic contributions to the decay rates (analytic chain);
3. the helicity amplitudes for basic processes;
4. the helicity amplitudes for accompanying bremsstrahlung processes;

Level 2: an s2n.f software generates automatically the fortran codes for $`\mathrm{\Gamma }^{(1)}=\mathrm{\Gamma }^{\mathrm{Born}}+\mathrm{\Gamma }^{\mathrm{Virt}}+\mathrm{\Gamma }^{\mathrm{Soft}}+\mathrm{\Gamma }^{\mathrm{Hard}}`$;

Level 3: an infrared rearrangement (or exponentiation) procedure should work here (it is still at the stage of development);

Level 4: a Monte Carlo event generator works out. For the time being, we have here a ‘manually written’ fortran code.

### 3.2 Some keywords, characterizing SANC

SANC is an Internet-based system (address http://brg.jinr.ru/) the use of which is supposed to be as simple as the usual surf on Internet. It is a database-based system. This means that there is a storage of source codes written in several languages, which talk to one another, being placed into a homogeneous environment written in JAVA (linker). It uses the precomputation or intermediate access principle which means that:
1) All the relevant one-loop diagrams are precomputed and stored. This greatly saves the CPU.
2) It has several intermediate ‘entries’ to bypass CPU consuming computations. The user may also access its final product, a MC generator. From the other side, the full chain of calculation ‘from the Lagrangian to realistic distribution’ can be worked out in real time upon the user request. This option is intended to demonstrate the self-reproducibility of the full chain of calculation at any time.

## 4 First versions of SANC

* v0.01, 3/01 realizes a part of analytic calculations of level 1 for the decays $`H`$ $`(Z,W)f_1\overline{f}_2`$ (demonstration of viability);
* several versions, v0.02c/d, were tested on the road toward realization of levels 2,4 for decays $`H(Z,W)f_1\overline{f}_2`$;
* the current version, v0.03 (summer’02), realizes the full chain of calculations, returns numbers and distributions for these decays at the one-loop level of precision (demonstration of workability);
* the next version, v0.10 (at work), should contain all FORM3 codes of level 1 for the processes $`2f2f`$ and the decays $`F3f`$.

### 4.1 Comparison with other codes

A lot of comparisons with the other codes were made. For instance, we compared numbers with those of the ZFITTER code for all channels of $`e^+e^{}f\overline{f}`$ with light fermion masses and found the following agreement: within 8–9 digits for the SFF; within 7–8 digits for the one-loop differential cross-sections $`d\sigma ^{(1)}/d\mathrm{cos}\vartheta `$; within 6–7 digits for the total cross section and $`\sigma ^{FB}`$.

Several comparisons were made for the $`e^+e^{}t\overline{t}`$ process. First is an internal one, between two fortran codes: s2n.f and eeffLib. The latter is a ‘manually written’ code aimed at benchmarking of s2n.f software. For the SFF and for the complete one-loop differential cross-section $`d\sigma ^{(1)}/d\mathrm{cos}\vartheta `$ the numbers agree within 12–13 digits. A comparison between s2n.f and FeynArts  for the one-loop cross-section without soft photons revealed an agreement within 11 digits. Finally, we compared s2n.f numbers with those of the topfit code  for the differential one-loop cross-section with soft photons and found an 8-digit agreement.

## 5 CONCLUDING REMARKS

SANC clearly is a long-term project. Its first phase is realized in demonstration version 0.10.

Its second phase assumes the creation of a complete software product, accessible via an Internet-based environment, and realizing the full chain of calculations at the one-loop level of precision including processes 2 $``$ 3 and decays 1 $``$ 4.