text:
# Diffusion-controlled annihilation 𝐴+𝐵→0 with initially separated reactants: The death of an 𝐴 particle island in the 𝐵 particle sea

\[

## Abstract

We consider the diffusion-controlled annihilation dynamics $`A+B0`$ with equal species diffusivities in the system where an island of particles $`A`$ is surrounded by the uniform sea of particles $`B`$. We show that once the initial number of particles in the island is large enough, then at any system’s dimensionality $`d`$ the death of the majority of particles occurs in the universal scaling regime within which $`4/5`$ of the particles die at the island expansion stage and the remaining $`1/5`$ at the stage of its subsequent contraction. In the quasistatic approximation the scaling of the reaction zone has been obtained for the cases of mean-field ($`dd_c`$) and fluctuation ($`d<d_c`$) dynamics of the front.

82.20.-w, 05.40.-a

\]

The reaction front formed in the $`A+B0`$ reaction-diffusion systems with initially separated reactants is of great interest since it represents a pattern for a wide spectrum of processes in physics, chemistry, and biology . Since the seminal article of Galfi and Racz  much work has been devoted to studying this problem by different approaches\[3-20\]. A standard way to treat the problem analytically is to solve the system of equations

$`a/t=D_A^2aR,b/t=D_B^2bR`$ (1)
with the initial state given by

$`a(x,0)=a_0\theta (x),b(x,0)=b_0\theta (x),`$ (2)
where $`a(x,t)`$ and $`b(x,t)`$ are the mean local concentrations of $`A^{}s`$ and $`B^{}s`$, $`R(x,t)`$ is the macroscopic reaction rate and $`\theta (x)`$ denotes the Heaviside step function, so that the $`A^{}s`$ are initially uniformly distributed on the left side ($`x<0`$), and the $`B^{}s`$ on the right side ($`x>0`$) of the initial boundary.

Dimensional  and renormalization group  analyses show that at $`dd_c=2`$ one can adopt the mean-field approximation $`R(x,t)=ka(x,t)b(x,t)`$($`k`$ being the reaction rate constant) with logarithmic corrections in the $`2d`$ case, whereas in $`1d`$ systems fluctuations play the dominant role and the explicit form of $`R`$ remains unknown. There are, however, several techniques which enable one to derive a lot of information from (1) even for $`d<d_c`$. They are focused on the long-time limit $`kt\mathrm{}`$ and include, as two basic concepts, the scaling ansatz  and the quasistatic approximation . According to the scaling ansatz (SA) the long-time behavior of the system inside the reaction zone may be represented in the form

$`R=R_f𝒬\left({\displaystyle \frac{xx_f}{w}}\right),`$ (3)
where $`x_ft^{1/2}`$ denotes the position of the reaction zone center, $`R_ft^\beta `$ is the height, and $`wt^\alpha `$ is the width of the reaction zone. At $`dd_c`$ the scaling exponents take the values $`\alpha =1/6,\beta =2/3`$ whereas at $`d=1`$ they appear equal to $`\alpha =1/4,\beta =3/4`$\[9-14\], so that at any $`d`$ the width of the reaction zone on the diffusion length scale $`\mathrm{\Lambda }_Dt^{1/2}`$ asymptotically unlimitedly contracts:

$$w/\mathrm{\Lambda }_D0\mathrm{as}t\mathrm{}.$$

The quasistatic approximation (QSA) consists in the assumption that for sufficiently long times the kinetics of the front is governed by two characteristic time scales. One time scale, $`t_J=(d\mathrm{ln}J/dt)^1`$, controls the rate of change in the diffusive current $`J=J_A=J_B`$ of particles arriving at the reaction zone. The second time scale, $`t_fw^2/D`$ is the equilibration time of the reaction front. Assuming $`t_f/t_J1`$ from the QSA in the mean-field case with $`D_{A,B}=D`$ it follows

$`R_fJ/w,w(D^2/Jk)^{1/3},`$ (4)
whereas in the $`1d`$ case $`w`$ acquires the $`k`$-independent form $`w(D/J)^{1/2}`$ . The most important feature of the QSA is that $`w`$ and $`R_f`$ depend on $`t`$ only through the time dependent boundary current $`J(t)`$, which can be calculated analytically without knowing the concrete form of $`𝒬`$, i.e., in fact, representing the reaction zone on the scale $`\mathrm{\Lambda }_D`$ in the form $`R(x,t)=J\delta (xx_f)`$. On the basis of the QSA the general description of spatiotemporal behavior of the system $`A+B0`$ has been obtained for arbitrary nonzero diffusion coefficients and initial species concentrations . These results are in full agreement with extended numerical calculations and experiments and were generalized recently to the cases of reversible reaction $`A+BC`$ \[21-24\] and to several more complex reactions \[25-29\].

Until now, however, the main attention has been focused on the systems with $`A`$ and $`B`$ domains unlimited extension, i.e., with unlimited number of $`A^{}s`$ and $`B^{}s`$ particles. The aim of this work is to develop a new line in the study of the $`A+B0`$ dynamics under the assumption that the particle number of one of the species, say $`A`$, per unit of the initial boundary is finite. More precisely, we will consider the problem on the dynamics of death of an $`A`$ particle island surrounded by the uniform sea of particles $`B`$ and will reveal the defining features of this process.

Let particles $`A`$ with concentration $`a_0`$ be uniformly distributed in the island $`x(L,L)`$ surrounded by the unlimited sea of particles $`B`$ with concentration $`b_0`$ on the left $`x(\mathrm{},L)`$ and on the right $`x(L,\mathrm{})`$ of the island. By symmetry our problem is reduced to the solution of the system (1) in the interval $`x[0,\mathrm{})`$ at the initial conditions

$$a(x,0)=a_0\theta (Lx),b(x,0)=b_0\theta (xL)$$

with the boundary conditions

$$(a,b)_{x=0}=0,b(\mathrm{},t)=b_0.$$

To simplify the problem essentially we will assume, as usually, $`D_A=D_B=D`$. Then by measuring the length, time, and concentration in units of $`L`$, $`L^2/D`$, and $`b_0`$, respectively, i.e. assuming $`L=D=b_0=1`$, and defining the difference concentration $`s(x,t)=a(x,t)b(x,t)`$ and the ratio of initial concentrations $`a_0=r`$ we come from (1) to the simple diffusion equation

$`s/t=^2s`$ (5)
with the initial conditions

$$s_0(x[0,1))=r,s_0(x(1,\mathrm{}))=1,$$

and the boundary conditions

$$s_{x=0}=0,s(\mathrm{},t)=1.$$

The solution to Eq.(5) has the form

$`s(x,t)={\displaystyle \frac{r+1}{2}}\left[\mathrm{erf}\left({\displaystyle \frac{1+x}{2\sqrt{t}}}\right)+\mathrm{erf}\left({\displaystyle \frac{1x}{2\sqrt{t}}}\right)\right]1,`$ (6)
whence, according to condition $`s(x_f,t)=0`$  there immediately follows the equation defining the law of motion of the reaction front center, $`x_f(t)`$

$`\mathrm{erf}\left({\displaystyle \frac{1+x_f}{2\sqrt{t}}}\right)+\mathrm{erf}\left({\displaystyle \frac{1x_f}{2\sqrt{t}}}\right)={\displaystyle \frac{2}{r+1}}.`$ (7)
Let us assume that $`dd_c`$ and the reaction rate constant $`k`$ is sufficiently large  so that at times $`tt_{GR}k^11`$ the annihilation goes to the scaling Galfi-Racz regime, i.e. in the vicinity of $`x_f`$ there forms a narrow reaction zone $`w/\mathrm{\Lambda }_D1`$. At $`t1`$ and $`|1x_f|1`$ from (7) in accordance with  we find

$$x_f=1+c_f\sqrt{t}+\mathrm{},$$

where $`\mathrm{erf}(c_f/2)=(r1)/(r+1)`$ and, hence, at $`r1`$ the island contracts, whereas at $`r>1`$ the island expands. By virtue of the fact that the number of particles in the island (per unit of the initial boundary) is finite, the stage of its expansion always goes to the stage of its subsequent contraction to end in an instant of time $`t_c`$, when the reaction front center approaches the origin of coordinates, $`x_f(t_c)=0`$:

$`\mathrm{erf}(1/2\sqrt{t_c})=1/(r+1).`$ (8)
According to (8) at $`r1`$ the ”lifetime” of the island $`t_c1`$ so the majority of the particles die at times $`t1`$, when the diffusive length exceeds appreciably the initial island size, $`\mathrm{\Lambda }_D1`$. The evolution of the island in such a large-$`t`$ regime is of principal interest to us here, and its analysis is the main goal of the present report.

In the limit $`r,t1`$ from (6) we find

$`s(x,t)={\displaystyle \frac{(r+1)}{\sqrt{\pi t}}}e^{x^2/4t}(1\chi )1,`$ (9)
where $`\chi =(1x^2/2t)/12t+\mathrm{}`$, and, hence, the law of the front motion is

$`x_f=2\sqrt{t}(1+ϵ)\mathrm{ln}^{1/2}[\left({\displaystyle \frac{r+1}{\sqrt{\pi t}}}\right)(1ϵ)],`$ (10)
where $`ϵ=1/12t+\mathrm{}`$. By assuming $`t1>t_{GR}`$, at the stage of the island expansion and at a considerable time period of its subsequent contraction, the duration of which will be estimated below, the reaction zone width is $`wx_f<\mathrm{\Lambda }_D`$, therefore in terms of the QSA we can take $`w/x_f0`$, and, hence, $`R(x,t)=J\delta (xx_f)`$ so that $`a=s,b=0`$ when $`x<x_f`$, and $`a=0,b=|s|`$ when $`x>x_f`$. Thus, by neglecting the width of the reaction zone, for the number of particles in the island we have $`N(t)=_0^{x_f}a(x,t)𝑑x=_0^{x_f}s(x,t)𝑑x`$ whence, with account taken of (9) and neglecting the contribution of $`\chi `$, we find

$`N(t)=(r+1)\mathrm{erf}(x_f/2\sqrt{t})x_f.`$ (11)
Let us now calculate the maximal amplitude $`x_f^M`$ of the island expansion and the time $`t_M`$ of its achievement. According to (7) for $`r>1`$ from condition $`\dot{x}_f=0`$ we have exactly

$$x_f^M=2t_M\mathrm{Arcoth}(x_f^M)$$

whence for $`r1`$ it follows

$`x_f^M=\sqrt{2t_M}(1+1/12t_M+\mathrm{}).`$ (12)
Neglecting the terms $`O(1/r^2)`$ from (10),(12) we find

$$t_M=(r+1)^2/\pi e,x_f^M=(r+1)\sqrt{2/\pi e}.$$

Comparing $`t_M`$ with the lifetime of the island (8),

$$t_c=(r+1)^2/\pi ,$$

we conclude that independently of $`r`$ the ratio $`t_c/t_M=e=\mathrm{const}`$. Going with allowance for this to the reduced coordinate $`\zeta =x/x_f^M`$ and time $`\tau =t/t_c`$ we, finally, come to scaling relations for the distribution of particles $`s(\zeta ,\tau )`$

$`s(\zeta ,\tau )=e^{\zeta ^2/2e\tau }/\sqrt{\tau }1,`$ (13)
for the center of the reaction front

$`\zeta _f=\sqrt{e\tau |\mathrm{ln}\tau |},`$ (14)
and for the number of particles in the island

$`N/N_0=\gamma _r𝒢(\tau ),`$ (15)
where $`\gamma _r=(r+1)/r1`$ as $`r\mathrm{}`$ and scaling function

$$𝒢(\tau )=\mathrm{erf}(\sqrt{|\mathrm{ln}\tau |/2})\sqrt{2\tau |\mathrm{ln}\tau |/\pi }.$$

From (13) and (15) we find that in the turning point $`\tau _M=1/e`$

$$a(0,\tau _M)=\sqrt{e}1=0.64872\mathrm{},$$

$$\gamma _r^1N_M/N_0=0.19886\mathrm{},$$

and, hence, independently of the initial number of particles, $`N_0=r`$, $`4/5`$ of the particles die at the stage of the island expansion and the remaining $`1/5`$ at the stage of its subsequent contraction. In Figs.1 and 2 are shown the calculated according to Eqs.(6),(7) plots of the behavior of the particle distribution $`|s(\zeta ,\tau )|`$, of the center of the reaction front $`\zeta _f(\tau )`$, and the number of particles $`\gamma _r^1N(\zeta _f)/N_0`$, which give the full insight into the scaling regime of the death of the island (13)-(15).

Let us come now to the question of the scaling behavior of the reaction rate $`R(x,t)`$ in the vicinity of the reaction front center $`x_f`$. Assuming according to  that within the applicability of the inequalities $`t_f/t_J1`$ and $`wx_f`$ the reaction rate $`R(x,t)`$ can be described in terms of the QSA by Eqs.(3),(4) to calculate $`R_f`$ and $`w`$ it only remains for us to calculate the diffusive current $`J(t)`$ of particles arriving at the reaction zone. According to (9),(11) we have

$$J=\dot{N}=s/x_{x=x_f}=x_f/2t$$

whence it follows

$`J/J_M=\sqrt{|\mathrm{ln}\tau |/e\tau },`$ (16)
where $`J_M=1/x_f^M=\sqrt{\pi e/2}/(r+1)`$. Substituting (16) into Eqs.(4) we come to the scaling of the reaction zone

$`R_f(J^4k)^{1/3}=R_f^M\left({\displaystyle \frac{|\mathrm{ln}\tau |}{e\tau }}\right)^{2/3}`$ (17)
$`w(Jk)^{1/3}=w_M\left({\displaystyle \frac{e\tau }{|\mathrm{ln}\tau |}}\right)^{1/6},`$ (18)
where $`R_f^M(k/r^4)^{1/3}`$ and $`w_M(r/k)^{1/3}`$ so that $`wR_fR_{global}=_0^{\mathrm{}}R𝑑x=J`$. From (18) it is seen that at the stage of the island contraction the logarithmic term becomes dominant, and as $`\tau 1`$, i.e., $`\delta \tau =1\tau 0`$, the width of the reaction front $`w`$ diverges as $`w(\delta \tau )^{1/6}`$. Comparing the characteristic times $`t_J=(d\mathrm{ln}J/dt)^1r^2\delta \tau `$ and $`t_fw^2(r/k)^{2/3}(\delta \tau )^{1/3}`$ we have $`t_f/t_J(r^2k)^{2/3}(\delta \tau )^{4/3}`$ whence it follows that the characteristic time at which the QSA is violated is

$$\delta \tau _Q1/r\sqrt{k}.$$

Comparing then $`w`$ and $`x_f`$ we obtain $`w/x_f\sqrt{t_f/t_J}(r^2k)^{1/3}(\delta \tau )^{2/3}`$ whence it follows that the width of the reaction zone becomes comparable with the island size at the times

$`\delta \tau _w\delta \tau _Q1/r\sqrt{k}.`$ (19)
From (19) it is seen that at sufficiently large values of $`k1`$ the characteristic times $`\delta \tau _w\delta \tau _Q1/r1`$ and, hence, practically all the particles die in the scaling regime (15).

The analogous estimations in the case of the $`k`$-independent fluctuation $`1d`$ regime with $`w1/\sqrt{J}`$ yield the expressions (17),(18) with $`R_f^{MF}r^{3/2},w_M^F\sqrt{r}`$, and exponents $`3/4`$ instead of $`2/3`$, and $`1/4`$ instead of $`1/6`$, respectively. In the end for the $`1d`$ case we have

$`\delta \tau _w^F\delta \tau _Q^F1/r^{2/3}.`$ (20)
According to (20), despite the much more ”blurred” structure of the $`1d`$ island, in this case, too, in the limit of sufficiently large $`r1`$ the majority of the particles die in the scaling regime (15).

In summary, the problem of death of an $`A`$ particle island in the $`B`$ particle sea at equal diffusivities of $`A^{}s`$ and $`B^{}s`$ particles has been first considered. It has been found that at sufficiently large initial number of particles in the island, $`r1`$, and sufficiently large reaction rate constant, $`k1`$, the death of the majority of the particles at any $`d`$ occurs in the universal scaling regime, and the most essential features of this regime have been revealed. The obtained results can have many applications, especially in surface science, and, in particular, they can provide fresh insight into the dynamics of the Ovchinnikov-Zeldovich hierarchic $`AB`$ structures . By analogy with  it may be expected that the analysis presented can be easily extended over the general case $`D_AD_B0`$ (note that in the static case $`D_B=0`$, which belongs to the separate universality class, with initially $`N_0`$ particles $`A`$ at the origin of coordinates, in  the formally coincident with (15) scaling, $`N=N_0(t/N_0^2)`$, was found). An investigation into the general case of nonzero diffusivities and a most interesting generalization for the anisotropic diffusion of reactants is expected to be presented in a future report.

I would like to thank Alexey Nekrasov for help in computer work. This research was financially supported by the RFBR through Grant No. 02-03-33122.


format: text

text:
# Rotational dynamics of CO solvated in small He clusters: a quantum Monte Carlo study

## I Introduction

Thanks to recent progresses in Helium nano-droplet isolation (HENDI) spectroscopy,Scoles2001  the infrared and microwave spectra of small molecules solvated in HeN clusters are now becoming accessible in the small- and intermediate-size regimes.Jaeger2002 ; Jaeger2003 ; McKellar2003-OCS ; McKellar2003-CO  In the particular case of carbon monoxide, the roto-vibrational spectrum of the molecule solvated in small He clusters (CO@HeN) has been recently studied in the size range $`N=220`$.McKellar2003-CO  The infrared spectrum in the 2145 cm-1 region of the C–O stretch consists of two $`R(0)`$ transitions which smoothly correlate with the $`a`$-type ($`K=00`$) and $`b`$-type ($`K=10`$) $`R(0)`$ lines of the binary complex, CO@He1. The series of $`b`$-type transitions—which starts off about 7 times stronger for $`N=1`$—progressively looses intensity as $`N`$ increases, until it disappears around $`N=78`$. Around this size, just before it disappears, the $`b`$-type line seems to split in two. Analogously, around $`N=15`$ the $`a`$-type line also seems to split, and the assignment of experimental lines becomes uncertain for larger clusters. Elucidating the relation existing between the position, number, and intensity of the rotational lines and the size and structure of the cluster is the principal goal of the present paper.

Computer simulations of quantum many-body systems have also considerably progressed in recent years, allowing in some cases to determine the low-lying spectrum of excited states. The rotational dynamics of small molecules solvated in He clusters and nanodroplets is one of these favorable instances. The scarcity of low-lying excited states typical of superfluid systems makes it possible in this case to extract information on the location and intensity of the spectral lines from an analysis of the time series generated by quantum Monte Carlo random walks.NoAntri-Cornell ; NoAntri-PRL ; POITSEH  The rotational spectrum of OCS@HeN has been studied along these lines in Refs. \[Whaley2003, ; NoAntri-2003, \]. Among the many different flavors of quantum Monte Carlo available in the literature, we adopt Reptation Quantum Monte CarloNoAntri-Cornell ; NoAntri-PRL  which we believe presents distinctive advantages in the present case and which will be briefly introduced in Sec. II. In Sec. III we present and discuss our results, whereas Sec. IV contains our conclusions.

## II Theory, algorithms, and technical details

Virtually all the ground-state quantum simulation methods are based on the prior knowledge of some approximate wave-function, $`\mathrm{\Phi }_0`$, for the system under study. In the Variational Monte Carlo Method (VMC) one contents oneself of this knowledge and the simulation simply aims at calculating the complicated multi-dimensional integrals which are needed to estimate ground-state approximate expectation values, $`\mathrm{\Phi }_0|\widehat{A}|\mathrm{\Phi }_0`$ (here and in the following quantum-mechanical operators are indicated with a hat, $`\widehat{}`$). To this end, a random walk in configuration space is generated according to the Langevin equation:

$$d\text{x}=ϵ\text{f}_0(\text{x})+d𝝃,$$
(1)
where $`\text{x}\{x_\alpha \}`$ indicates the coordinates of the system, $`ϵ`$ is the step of time discretization,

$$\text{f}_0(\text{x})=2\frac{\mathrm{log}\left(\mathrm{\Phi }_0(\text{x})\right)}{\text{x}},$$
(2)
and $`d𝝃`$ is a Gaussian random variable of variance $`2ϵ`$: $`d\xi _\alpha d\xi _\beta =2ϵ\delta _{\alpha \beta }`$. Approximate ground-state quantum expectation values are then estimated as time averages over the random walk, Eq. (1).

Within Reptation Quantum Monte Carlo (RQMC), exact ground-state expectation values and imaginary-time correlation functions are calculated as appropriate derivatives of the pseudo partition function, in the low temperature (large $`T`$) limit:

$$𝒵_0=\mathrm{\Phi }_0|e^{T\widehat{H}}|\mathrm{\Phi }_0,$$
(3)
where $`\widehat{H}`$ is the Hamiltonian of the system. By breaking the time $`T`$ into $`P`$ intervals of length $`ϵ=T/P`$, Eq. (3) can be given a path-integral representation:

$$𝒵_0\mathrm{\Phi }_0(\text{x}_0)\mathrm{\Pi }_{i=0}^{P1}\text{x}_i|e^{ϵ\widehat{H}}|\text{x}_{i+1}\mathrm{\Phi }_0(\text{x}_P)d^{P+1}\text{x}.$$
(4)
For the relatively small systems considered here, it is sufficient to use the primitive approximation to the imaginary-time propagator:

$$\text{x}|e^{ϵ\widehat{H}}|\text{y}\mathrm{e}^{(\text{x}\text{y})^2/2ϵϵ[V(\text{x})+V(\text{y})]/2}+𝒪(ϵ^3),$$
(5)
where $`V(\text{x})`$ is the potential energy at point x. The dynamical variables of the statistical-mechanical system whose partition function is given by Eq. (4) are segments of the VMC random walk generated from Eq. (1), $`\text{x}(\tau )`$ of length $`T`$, which we call reptiles. As the random walk proceeds, the reptile is allowed to creep back and forth: new configurations of the reptile are accepted or rejected according to a Metropolis test made on the integrand of the path-integral representation of $`𝒵_0`$, Eq. (4). It can be demonstratedNoAntri-Cornell ; NoAntri-PRL  that in the large-$`T`$ limit the sample of reptile configurations thus generated is such that the sample average of quantities like

$$𝒜[\text{x}(\tau )]=\frac{1}{T}_0^TA(\text{x}(\tau ))𝑑\tau $$
(6)
converges without any systematic bias (but those due to the finite values of $`ϵ`$ and $`T`$) to $`\widehat{A}=\mathrm{\Psi }_0|\widehat{A}|\mathrm{\Psi }_0`$, $`\mathrm{\Psi }_0`$ being the exact ground-state wavefunction of the system. Even more interesting is the fact that sample averages of reptile time correlations,

$$𝒞_A(\tau )=\frac{1}{T\tau }_0^{T\tau }A(\text{x}(\tau ^{}))A(\text{x}(\tau ^{}+\tau ))𝑑\tau ^{},$$
(7)
provide equally unbiased estimates of the corresponding quantum correlation functions in imaginary time:

$`C_{\widehat{A}}(it)`$ $`=`$ $`\mathrm{\Psi }_0|\widehat{A}(it)\widehat{A}(0)|\mathrm{\Psi }_0`$ (8)
$``$ $`\mathrm{\Psi }_0|\mathrm{e}^{\widehat{H}t}\widehat{A}\mathrm{e}^{\widehat{H}t}\widehat{A}|\mathrm{\Psi }_0,`$
$`𝒞_A(\tau )C_{\widehat{A}}(i\tau )`$. The absorption spectrum of a molecule solvated in a non polar environment is given by the Fourier transform of the autocorrelation function of its electric dipole, d:

$`I(\omega )`$ $``$ $`2\pi {\displaystyle \underset{n}{}}|\mathrm{\Psi }_0|\widehat{\text{d}}|\mathrm{\Psi }_n|^2\delta (E_nE_0\omega )`$ (9)
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\mathrm{e}^{i\omega t}\widehat{\text{d}}(t)\widehat{\text{d}}(0)𝑑t,`$
where $`\mathrm{\Psi }_0`$ and $`\mathrm{\Psi }_n`$ are ground- and excited-state wavefunctions of the system respectively, and $`E_0`$ and $`E_n`$ the corresponding energies. The dipole of a linear molecule—such as CO—is oriented along its axis, so that the optical activity is essentially determined by the autocorrelation function of the molecular orientation versor: $`c(t)C_{\widehat{\text{n}}}(t)=\mathrm{\Psi }_0|\mathrm{e}^{i\widehat{H}t}\widehat{\text{n}}\mathrm{e}^{i\widehat{H}t}\widehat{\text{n}}|\mathrm{\Psi }_0`$. We have seen that RQMC gives easy access to the analytic continuation to imaginary time of correlation functions of this kind. From now on, when referring to time correlation functions, we will mean reptile time correlations, i.e. quantum correlation functions in imaginary time. Continuation to imaginary time transforms the oscillatory behavior of the real-time correlation function—which is responsible for the $`\delta `$-like peaks in its Fourier transform—into a sum of decaying exponentials whose decay constants are the excitation energies, $`E_nE_0`$, and whose spectral weights are proportional to the absorption oscillator strengths, $`|\mathrm{\Psi }_0|\text{d}|\mathrm{\Psi }_n|^2`$. Dipole selection rules imply that only states with $`J=1`$ can be optically excited from the ground state which has $`J=0`$. Information on excited states with different angular momenta, $`J`$, can be easily obtained from the multipole correlation functions, $`c_J(\tau )`$, defined as the reptile time correlations of the Legendre polynomials:

$`c_J(\tau )`$ $`=`$ $`P_J(\text{n}(\tau )\text{n}(0))`$ (10)
$``$ $`{\displaystyle \frac{4\pi }{2J+1}}{\displaystyle \underset{M=J}{\overset{J}{}}}Y_{JM}^{}\left(\text{n}(\tau )\right)Y_{JM}(\text{n}(0)).`$
Both the He-He and the He-CO interactions used here are derived from accurate quantum-chemical calculations.Korona ; Heijmen  The CO molecule is allowed to perform translational and rotational motions, but it is assumed to be rigid. The trial wavefunction is chosen to be of the Jastrow form:

$$\mathrm{\Phi }_0=\mathrm{exp}\left[\underset{i=1}{\overset{N}{}}𝒰_1(r_i,\theta _i)\underset{i<j}{\overset{N}{}}𝒰_2(r_{ij})\right],$$
(11)
where $`\text{r}_i`$ is the position of the i-th atom with respect to the center of mass of the molecule, $`r_i=|\text{r}_i|`$, $`\theta _i`$ is the angle between the molecular axis and $`\text{r}_i`$, and $`r_{ij}`$ is the distance between the i-th and the j-th helium atoms. $`𝒰_1`$ is expressed as a sum of five products of radial functions times Legendre polynomials. All radial functions (including $`𝒰_2`$) are optimized independently for each cluster size with respect to a total of 27 variational parameters. The propagation time is set to $`T=1\mathrm{K}^1`$, with a time step of $`ϵ=10^3\mathrm{K}^1`$. The effects of the length of the time step and of the projection time have been estimated by test simulations performed by halving the former or doubling the latter. These effects were barely detectable on the total energy, and very small on the excitation energies discussed below (we estimate that more converged simulations would actually improve the already excellent agreement with experimentally observed spectra).

The estimate of excitation energies and spectral weights from imaginary-time correlations amounts to performing an inverse Laplace transform, a notoriously ill-conditioned operation which is severely hindered by statistical noise.gubernatis  For each value of $`J`$, we extract the value of the two lowest-lying excitation energies, $`ϵ_{a,b}^J`$i.e. the two smallest decay constants in $`c_J(\tau )`$—as well as the corresponding spectral weights, $`A_{a,b}^J`$, from a fit of $`c_J(\tau )`$ to a linear combination of three decaying exponentials. This fitting procedure does not solve in general the problem of obtaining the spectrum from a noisy imaginary-time correlation function. However, if we know in advance that very few strong peaks, well separated in energy, nearly exhaust the entire spectral weight, their position and strength can be reliably estimated from this multi-exponential fit. In the present study these favorable conditions are usually met, although the limitations of the procedure will show in some cases, as discussed below.

## III Results and discussion

RQMC simulations have been performed for CO@HeN clusters in the size range $`N=130`$. In Fig. 1 we report the values of the He atomic binding energy, $`\mathrm{\Delta }E_N=E_{N1}E_N`$, as a function of the cluster size. $`\mathrm{\Delta }E_N`$ first increases up to $`N=45`$, and it stays roughly constant in the range $`N=58`$; from this size on $`\mathrm{\Delta }E_N`$ starts decreasing, first slowly, then, from $`N=1011`$, rapidly down to a minimum at $`N=19`$. For $`N>19`$ $`\mathrm{\Delta }E_N`$ increases again and slowly tends to the nanodroplet regime (where it coincides with the bulk chemical potential, $`\mu =7.4\mathrm{K}`$ NoAntri-PRL ) which is however attained for much larger cluster sizes than explored here.paesani  This behavior can be understood by comparing the shape of the CO–He potential energy function, $`v(\text{r})`$, with the incremental atomic density distributions, $`\mathrm{\Delta }\rho _N(\text{r})=\rho _N(\text{r})\rho _{N1}(\text{r})`$, where $`\rho _N`$ is the expectation value of the He density operator:

$$\widehat{\rho }(\text{r})=\underset{i=1}{\overset{N}{}}\delta (\text{r}\text{r}_i)$$
(12)
(see Fig. 2). For very small $`N`$ the atomic binding energy is dominated by the He–CO attraction which is strongest in a well located atop the oxygen atom. As He atoms fill this well, $`\mathrm{\Delta }E_N`$ first slightly increases, as a consequence of the attractive He–He interaction, then, for larger $`N`$, the increased He–He interaction is counter-balanced by the spill-out of He atoms off the main attractive well, until for $`N9`$ the reduction of the He–CO interaction overcomes the increased attraction and the binding energy starts decreasing steeply. For $`N`$ in the range 10–14 He density accumulates towards the C pole, while, around $`N=15`$, the first solvation shell is completed and the differential atomic density, $`\mathrm{\Delta }\rho _N`$, is considerably more diffuse starting from $`N=16`$. $`\mathrm{\Delta }E_N`$ reaches a minimum at $`N=19`$. For larger sizes, the trend in the atomic binding energy is dominated by the increase of the He–He attraction related to the increase of the cluster size, until it will converge to the bulk chemical potential.

In Fig. 3 we report the positions and spectral weights of the rotational lines, as functions of the cluster size, $`N`$. In the size range $`N=19`$, analysis of the dipole time correlations clearly reveals the presence of two peaks, with the weight of the higher–energy ($`b`$-type) rapidly decreasing by almost a factor 2. Note that the sum of the spectral weights of these two lines nicely sums to one, indicating that they exhaust all the oscillator strength available for optical transitions originating from the ground state. For $`N`$ between 10 and 12 (shaded area in Fig. 3) the situation is less clear. As the weight of the $`b`$-type line drops to zero, the statistical noise on its position grows enormously. Furthermore the multi-exponential fit introduces some ambiguity, as the results are somewhat sensitive to the number of terms in the sum. However the important information that one line disappears between $`N=10`$ and 12 is clear. For larger $`N`$ only one relevant line remains, and the robustness of the fitting procedure is recovered, with the minor exception of the sizes around $`N=16`$, where the minimum of the $`\chi ^2`$ appears to be less sharp, possibly correlating with the splitting of the line observed in the infrared spectra for $`N=15`$ (see below).

In the upper panel of Fig. 4 we compare the rotational structure of the observed infrared (vibrational) spectrumMcKellar2003-CO  with the rotational excitation energies calculated in this work. Experimental data are referred to the center, $`\nu _0`$, of the vibrational band for $`N=0`$ (CO monomer). In order to better compare our predictions with experiments, we have corrected the former with an estimate of the vibrational shift, $`\mathrm{\Delta }\nu _0`$, i.e. the displacement of the vibrational band origin as a function of the number of He atoms. The vibrational shift can be calculated as the difference in the total energy of the cluster obtained with two slightly different potentials,Heijmen  $`v_{00}`$ and $`v_{11}`$, representing the interaction of a He atom with the CO molecule in its vibrational ground state and first excited states, respectively. Our estimate of the vibrational shift as a function of the cluster size is reported in the lower panel of Fig. 4. Since the evaluation of a small difference between two large energies is computationally demanding for large clusters, $`\mathrm{\Delta }\nu _0`$ has been evaluated perturbatively with respect to the difference $`v_{00}v_{11}.`$paesani  We have used the vibrational shift calculated in Ref. paesani,  after verifying on small clusters that the perturbative treatment is reliable. The agreement between our results and experiments is remarkable. Some of the features of the observed spectrum, however, call for a deeper understanding and theoretical investigation. Two questions, in particular, naturally arise. Why two peaks are observed in the small-size regime, and what determines the disappearance of one of them at $`N=8`$? What determines the split of the higher-frequency ($`b`$-type) line at $`N=7`$ and of the lower-frequency ($`a`$-type) one at $`N=15`$?

The existence of two lines for small $`N`$ is likely due to a larger asymmetry of the cluster in this regime. If the CO@HeN complex is described as a rigid rotor, in fact, one would have one rotational line originating from a $`J=0`$ ground state if the complex has cylindrical symmetry, while this line would double if some of the atomic density accumulates in a longitudinal protrusion. The inertia of the complex would in this case be larger for a rotation about an axis perpendicular to a plane containing the protrusion (end-over-end rotation) than about an axis lying on such a plane. Given that the He density in the ground state of CO@HeN is cylindrically symmetric, any departure from this symmetry can only show up in higher correlation functions. The situation is conceptually similar to that of a fluid whose density is homogeneous and whose structure at the atomic scale is reflected in the pair correlation function. Analogously, we define an atomic angular correlation function, $`C(\varphi )`$, as the probability of finding two He atoms which form a dihedral angle $`\varphi `$ with respect to the molecular axis:

$$C(\varphi )=\frac{1}{N(N1)}\underset{ij}{}\delta (\varphi _i\varphi _j\varphi ).$$
(13)
In the upper panel of Fig. 5 we show $`C(\varphi )`$, for different cluster sizes. The depletion of $`C`$ for $`\varphi `$ larger than $`\pi /2`$, clearly visible for $`N=3`$ (green circles), indicates a tendency of the He atoms to cluster on a same side of the molecular axis. For larger clusters, however, this effect weakens to the extent that it becomes difficult to disentangle from the structural information related to the He–He interaction (the dimple at small $`\varphi `$ and the subsequent maximum around $`\varphi =0.30.4\pi `$). A more sensitive measure of the propensity of He atoms to cluster on a side of the molecule is given by the integral of $`C(\varphi )`$ from 0 to $`\frac{\pi }{2}`$,

$$M=_0^{\frac{\pi }{2}}C(\varphi )𝑑\varphi \frac{1}{2}.$$
(14)
In the lower panel of Fig. 5 we display $`M`$ as a function of the cluster size, $`N`$: one sees that $`M`$ decreases with $`N`$ and reaches a minimum at $`N=14`$. This is the size at which the first solvation shell is completed, and the cluster asymmetry increases again when the second shell starts to build. The rotational spectrum of the solvated molecule, however, is insensitive to this asymmetry for clusters of this and larger sizes because the motion of He atoms in the second and outer solvation shells is decoupled from that in the first and from molecular rotation. The existence of a longitudinal asymmetry is a necessary condition for the doubling of the rotational line. Whether or not this condition is also sufficient depends on the dynamics: if quantum fluctuations make the motion of the protrusion around the molecular axis fast with respect to the molecular rotation, then the asymmetry is effectively washed out. The existence of two lines in the rotational spectrum of the molecule implies therefore that an asymmetry in the classical distribution of He atoms around the molecular axis exists; that the molecular inertia is sensitive to this asymmetry (the protrusion can be ‘dragged’ along the molecular rotation); and that the motion of this protrusion around the molecular axis is not adiabatically decoupled from the molecular rotation.

In order to better characterize the motion of He atoms around the molecule and the coupling of this motion to molecular rotation, we examine the imaginary-time correlations of the versor, $`𝐮`$, of the He center of mass, $`𝐫_{CM}`$, relative to the molecular center of mass:

$$𝒞_𝐮(\tau )=𝐮(\tau )𝐮(0).$$
(15)
For the binary complex, He–CO, $`𝐫_{CM}`$ coincides with the position of the helium atom, and we expect its angular dynamics to be strongly correlated to the molecular rotation, at least in the end-over-end mode. In Fig. 6 we report the frequency of the slowest mode appearing in the spectral analysis of $`𝒞_𝐮(\tau )`$, $`ϵ_𝐮`$, as a function of $`N`$, and compare it with the corresponding frequencies of the molecular rotation. We see that for cluster sizes up to $`N=910`$, $`ϵ_𝐮`$ is degenerate with the $`a`$-type frequency in the molecular rotational spectrum, with a spectral weight which passes from $`A_𝐮1`$ for $`N=1`$ to $`A_𝐮0.7`$ for $`N=10`$. These findings are a manifestation of the fact that He atoms are dragged along the slowest, end-over-end, rotation of the solvated molecule, and that the effect of this dragging decreases when more He states with $`J=1`$ become available and subtract spectral weight to the slowest mode. For $`N>10`$, $`ϵ_𝐮`$ further increases and departs from $`ϵ_a`$, indicating an effective decoupling of the two kinds of motion. In this regime, the effective rotational constant $`B`$ of the solvated molecule is almost independent of the cluster size. Free molecular rotation with an increased moment of inertia with respect to the gas phase is the typical signature of superfluid behavior in He nanodroplets. Extrapolating the result obtained for $`N`$ up to 30 to the nanodroplet limit, we predict a renormalization factor of the $`B`$ value of 0.78. The lowest atomic mode, $`ϵ_𝐮`$, slows down again for $`N=15`$. This is due, however to the slow He motion in the second solvation shell which hardly affects the rotation of the solvated molecule. Although the resolution that can be achieved with our simulations is not sufficient to detect the doubling of the $`a`$ and $`b`$ lines which is experimentally observed for $`N=15`$ and $`N=7`$ respectively, it is interesting to notice that the former occurs in correspondence with the crossing between $`ϵ_𝐮`$ and $`ϵ_a`$, possibly due to the resonant interaction between the two modes. It is tempting to assume that a similar mechanism may be responsible for the doubling of the $`b`$ line at $`N=7`$, involving however higher-energy He states. A deeper study of the He dynamics would clarify this point.

## IV Conclusions

Computer simulations of quantum many-body systems have reached such a degree of sophistication and reliability that in some cases they can be used to provide information, complementary to that which can be obtained in the laboratory, on the dynamical processes probed spectroscopically.

In the case of small polar molecules solvated in He clusters, for instance, the calculation of the time autocorrelation of the molecular dipole (which is the quantity directly coupled to the experimental probe) allows to reproduce rather accurately the roto-vibrational excitation energies which are now becoming experimentally accessible for small clusters ($`N=120`$). Even more importantly, computer simulations give direct access to quantitities and features (such as, e.g., static and dynamic properties of the He matrix) which are not accessible to the experiment, and whose knowledge provides the basis for understanding the relation between structure and dynamics in these confined boson systems.

In the specific case of CO@HeN, which is the subject of the present study and of a recent infrared spectroscopy experiment,McKellar2003-CO  the presence of two spectral lines—$`a`$-type and $`b`$-type, evolving respectively from the end-over-end and from the free-molecule rotations of the binary complex—is related to the propensity of the He atoms to cluster on a same side of the molecular axis, which we measure by an angular pair distribution function: as more He atoms progressively fill the first solvation shell, their clustering propensity weakens; the CO impurity gets more isotropically coated, looses a preferred axis for the free-molecule mode, and the $`b`$–type line disappears.

The time autocorrelation of the versor of the He center of mass provides dynamical information on the He atoms in excited states with $`J=1`$. We find a substantial spectral weight on a He mode whose energy, $`ϵ_𝐮`$, is degenerate with the $`a`$-type line for $`N`$ up to about ten. This indicates that some of the He density is dragged along by the molecular rotation—in other words, part of the angular momentum in the cluster mode involving molecular rotation is carried by the He atoms. We also find that for larger clusters the molecular rotation effectively decouples from this He mode, and its energy $`ϵ_a`$ becomes essentially independent of the number of He atoms. Based on the nearly constant value of $`ϵ_a`$ in the range of $`N`$ between 15 and 30, well beyond completion of the first solvation shell, we predict the effective rotational constant in the nanodroplet limit to be smaller by a factor 0.78 than its gas phase value.

###### Acknowledgements.

We would like to thank A.R.W. McKellar for providing us with a preprint of Ref. \[McKellar2003-CO, \] prior to publication. We are grateful to G. Scoles for bringing that work to our attention, for his continuous interest in our work, and for many useful discussions. Last but not least, we wish to thank S. Fantoni for his encouragement and for a critical reading of our manuscript. This work has been partial supported by the Italian MIUR through PRIN.


text:
# Cluster structures within Fermionic Molecular Dynamics

## 1 FERMIONIC MOLECULAR DYNAMICS

In the Fermionic Molecular Dynamics (FMD) model  the A-body state is given as a Slater determinant $`|Q>`$ of single-particle states $`|q_i>`$

$$|Q>=𝒜\left\{|q_1>\mathrm{}|q_A>\right\}.$$
(1)
The single-particle wave functions are described by Gaussian wave packets localized in phase-space

$$<𝒙|q>=\underset{i}{}c_i\mathrm{exp}\left\{\frac{(𝒙𝒃_i)^2}{2a_i}\right\}|\chi _i>|\xi >.$$
(2)
In contrast to the AMD approach  in the FMD model the width parameter $`a`$ is treated as a complex variational parameter that can be different for each single-particle state. The spin orientation given by the spinor $`|\chi >`$ is also treated as a variational parameter. For light nuclei the description can be improved by using a superposition of two Gaussian wave packets.

The FMD many-particle state is determined by minimizing the intrinsic energy

$$E\left[|Q>\right]=\frac{<Q\left|H_{\mathrm{𝑒𝑓𝑓}}T_{\mathrm{𝑐𝑚}}\right|Q>}{<Q|Q>}$$
(3)
with respect to the parameters of all single-particle states. This corresponds to a Hartree-Fock calculation using the FMD single-particle basis. To improve upon this mean-field result we restore the symmetries of the Hamiltonian and perform a parity and angular momentum projection. This defines the projection after variation (PAV) result.

The effects of the projection can be quite large and a variation after projection calculation (VAP) should be done . But as a full VAP calculation is very expensive we perform VAP calculations in the spirit of the Generator Coordinate Method. The intrinsic energy is minimized under constraints on collective coordinates like radius, dipole, quadrupole or octupole moments obtaining a set of many-body configurations. The VAP minimum is then obtained by minimizing the projected energies of these configurations as a function of the constraints. These configurations can then also be used in a multiconfiguration calculation (MC) where the projected Hamiltonian is diagonalized within the set of these configurations. This provides us with a well defined procedure to study the properties of the ground and the excited states.

## 2 EFFECTIVE INTERACTION

We use an effective interaction that is derived from the realistic Bonn or Argonne interactions by means of the Unitary Correlation Operator Method (UCOM) . The correlated interaction includes the short-range radial and tensor correlations induced by the repulsive core and the tensor force. The correlated interaction no longer connects to high momenta and can thus be used directly within product-state model spaces. The effects of three-body correlations and genuine three-body forces are at this stage simulated by an additional two-body force that contains momentum-dependence and spin-orbit forces. This correction term is adjusted to reproduce the binding energies and radii of $`{}_{}{}^{4}\mathrm{He}`$, $`{}_{}{}^{16}\mathrm{O}`$, $`{}_{}{}^{40}\mathrm{Ca}`$, $`{}_{}{}^{24}\mathrm{O}`$ and $`{}_{}{}^{48}\mathrm{Ca}`$. Here we find for $`{}_{}{}^{16}\mathrm{O}`$ and $`{}_{}{}^{40}\mathrm{Ca}`$ that after angular momentum projection the tetrahedral configurations (see Figure 1) with $`\alpha `$-clusters are favored energetically versus the spherical shell model configurations by about 5 MeV. Altogether we need a 15% correction to the *ab-initio* correlated two-body interaction.

## 3 GROUND AND EXCITED STATES OF 12C

The question of the clustering in $`{}_{}{}^{12}\mathrm{C}`$ has a long history and has been studied in the AMD framework . Our effective interaction is fitted to binding energies and radii of $`{}_{}{}^{4}\mathrm{He}`$ and $`{}_{}{}^{16}\mathrm{O}`$ but has not been adjusted to $`{}_{}{}^{12}\mathrm{C}`$. We perform the FMD variation and obtain the $`(1s_{1/2})^4(1p_{3/2})^{12}`$ shell model configuration (V), displayed in the left part of Figure 2. As it is spherical the angular momentum projection (PAV) does not affect the result. The binding energy obtained with this configuration is 7.5 MeV too low compared to the experimental value (see Table 1). The obtained radius is too small. We can also restrict the intrinsic state to pure $`\alpha `$-cluster configurations. In these configurations the spin-orbit interaction does not contribute and already in the mean-field calculation a triangular configuration is preferred. In a VAP calculation where the positions and the widths of the $`\alpha `$’s are the generator coordinates we obtain 4.3 MeV less binding than with the shell model configuration and the radius is too large for the pure $`\alpha `$-cluster configuration. If we perform a variation after projection calculation with the radius and the octupole moment as generator coordinates we obtain the triangular configuration with polarized $`\alpha `$-clusters denoted by VAP on the l.h.s. of Figure 2. It represents an intermediate between the VAP $`\alpha `$-cluster and the shell model (V/PAV) shape. The description of this polarization is strongly affected by using two Gaussians per single-particle state. We obtain 4.2 MeV more binding compared to a VAP calculation with only one Gaussian per single-particle state. In case of the shell-model configuration the energy gain is only 1.5 MeV. The radius obtained in the VAP calculation is a little bit too small when compared to experiment. For the intrinsically deformed VAP configurations we can project out not only a $`0^+`$ but also $`2^+`$, $`4^+`$ and $`3^{}`$ states (as shown in the VAP column of the spectrum in Figure 2) and can calculate the $`B(E2)`$ value for the $`0^+2^+`$ transition using bare charges. The $`B(E2)`$ value obtained for the VAP configuration is too small compared to the experimental value.

To improve the description we perform multiconfiguration calculations. For the additional configurations quadrupole and octupole moments are used as generator coordinates. In Figure 2 the configurations giving the biggest contributions for the first $`3^{}`$ and the second and third $`0^+`$ states are shown. In a multiconfiguration calculation containing the VAP and these three intrinsic states the binding energy of the ground state is increased by 1.5 MeV. For the electromagnetic properties the effect of admixing these less compact configurations is more pronounced. The radius is increased and the $`B(E2)`$ value is now in agreement with the experimental value. The spectrum obtained with these four intrinsic states is shown in the Multiconfig(4) column of Figure 2. Using 25 configurations the energies for the ground state and the lowest excited states change only slightly. For the higher excited states the effect is much bigger. Comparing to the experimental spectrum our results for the excited $`0^+`$ and $`2^+`$ states are already quite reasonable. Whereas the second $`0^+`$ state seems to be mainly built of three $`\alpha `$ configurations, the third $`0^+`$ appears to be mainly $`{}_{}{}^{8}\mathrm{Be}`$ plus $`\alpha `$. A detailed study of these $`0^+`$ and also $`2^+`$ states will be done in the future. We expect that even more extended configurations will have to be considered, taking into account the proposed Bose condensed nature of these states around the three $`\alpha `$ threshold .

## 4 HELIUM ISOTOPES

To calculate the properties of the He isotopes we performed FMD calculations with two Gaussians per single-particle state using the electric dipole and the mass quadrupole moment as constraints for variation after projection (VAP) and multiconfiguration (Multiconf) calculations. In Table 2 the calculated binding energies and radii are shown together with the experimental values. In Figure 3 the intrinsic configurations corresponding to the VAP minima are displayed. With the two Gaussians the FMD wave function is able to describe the neutron halo around the $`\alpha `$-core. Using only a single Gaussian an imperfect halo together with a distorted $`\alpha `$-core is observed. Compared to the calculation using a single Gaussian the neutron-rich He isotopes gain about 3-5 MeV in the PAV case. For the $`{}_{}{}^{6}\mathrm{He}`$ nucleus we find a quadrupole deformed VAP minimum where the two halo neutrons are sitting on opposite sides of the core and a more bound (by about 1.1 MeV) dipole deformed VAP minimum where the neutrons are located on the same side of the core (see middle row of Figure 3). For the other isotopes we only find minima with respect to the dipole deformation. The multiconfiguration calculations that include shapes of different electric dipole moments reproduce the experimental staggering of the binding energies. The large measured matter radii of $`{}_{}{}^{6}\mathrm{He}`$ and $`{}_{}{}^{8}\mathrm{He}`$ are also well described by the multiconfiguration calculations. We also observe a significant increase in the charge radii of the neutron rich He isotopes compared to $`{}_{}{}^{4}\mathrm{He}`$. This is explained by the zero-point motion of the more or less unchanged $`\alpha `$-core with respect to the center of mass of the nucleus.

We conclude that the zero-point motion of paired neutrons in a soft-dipole mode seen in VAP and multiconfiguration mixing calculations is an essential ingredient for the understanding of neutron-rich He isotopes and probably also for other neutron halo nuclei.


text:
# Single-atom density of states of an optical lattice

## Abstract

We consider a single atom in an optical lattice, subject to a harmonic trapping potential. The problem is treated in the tight-binding approximation, with an extra parameter $`\kappa `$ denoting the strength of the harmonic trap. It is shown that the $`\kappa 0`$ limit of this problem is singular, in the sense that the density of states for a very shallow trap ($`\kappa 0`$) is qualitatively different from that of a translationally invariant lattice ($`\kappa =0`$). The physics of this difference is discussed, and densities of states and wave functions are exhibited and explained.

The achievement of degeneracy in atomic gases has led to a surge of activity. One reason is that these systems provide a test-bed for theories developed in the context of solid-state physics. Indeed there have been experimental observations of Bose-Einstein condensation (BEC) Anderson et al. (1995), vortex lattices Abo-Shaeer et al. (2001); Haljan et al. (2001), superfluid-insulator transitions Greiner et al. (2002), fermion degeneracy DeMarco and Jin (1999) and, very recently, advances have been made towards Bardeen-Cooper-Schrieffer (BCS) superfluidity Greiner et al. . A great advantage of these systems is one’s unprecedented control of the experimental situation. By varying the magnetic field strength, particle-particle interactions can be altered. Moreover, one can use lasers to create an optical lattice whose parameters can be changed at will Greiner et al. (2001). The latter scenario has been the subject of recent theoretical attention Scott et al. (2002); Büchler et al. ; Ho et al. ; Scott et al. (2003).

One key difference between these systems and the usual models of the solid state is that they are not translationally invariant, due to the trapping potential necessary to confine the gas. However, in the usual set-up, when the trap is the only external potential, the local density approximation can be employed, and thermodynamic quantities can be simply related to their values in the absence of the trap Damle et al. (1996). Trapped atoms are, therefore, a good model of the continuum many-particle system — analogous, in the fermion case, to ‘jellium’, the ideal metal of many-electron theory. One thus expects that, when the periodic optical lattice potential is added, an experimental model of electrons in real crystals can be obtained.

In this letter we take a first step towards testing this expectation by obtaining the single-particle density of states (DOS). This is the starting point for understanding many key features of degenerate quantum systems. For example, in an ideal Bose gas it determines whether there is BEC Nozières (1995), while in a weakly attractive Fermi gas it gives the critical temperature for BCS superfluidity de Gennes (1966). The DOS becomes a smooth curve only in the limit in which the trap is very weak, so that the spectrum becomes a continuum. As we shall see, this limit is singular, in the sense that the DOS is qualitatively different from that of an infinite, translationally invariant lattice.

*Model.* We consider a single particle, with no internal degrees of freedom, moving in an infinite tight-binding lattice with one orbital per site and hopping between nearest neighbours. The particle is confined to a finite region by a harmonic trapping potential. For simplicity, we assume a $`d`$-dimensional hypercubic lattice. For $`d>1`$, the problem is separable, so we shall focus our attention on the 1D case, indicating how the results for $`d2`$ follow as the need arises. The Hamiltonian for $`d=1`$ is

$$H=t\underset{j}{}\left(|jj+1|+\text{H.c.}\right)+\frac{1}{2}\kappa \underset{j}{}\left(aj\right)^2|jj|,$$
(1)
where $`j`$ is a site label, $`a`$ the lattice constant, $`\kappa `$ the strength of the trap and $`t`$ the nearest-neighbour hopping integral. (This model also occurs in the related subject of the dynamical diffraction of atoms by static light fields; see, for example, ref. O’Dell (2001) and references therein.) Note that the trap is centred at a lattice site; see refs. Chalbaud et al. (1986); Gallinar and Chalbaud (1991) for a discussion of the non-trivial effect of incommensuration. The model can be generalised also by considering anisotropic lattices and trapping potentials.

Before describing the properties of the model, let us discuss its validity. In experiments Greiner et al. (2001), the lattice is generated by counter-propagating laser fields, giving rise to a static, periodic potential Dahan et al. (1996) of the form $`U(x)=U_0\left[\mathrm{cos}\frac{2\pi x}{a}+1\right]`$. Eq. (1) is valid when the wells of this potential are sufficiently deep: $`U_0\mathrm{}^2/ma^2`$, where $`m`$ is the mass of one atom. In this intense-laser regime, one can approximate the orbitals at the bottom of each well by those of a harmonic oscillator, and the usual method Ashcroft and Mermin (1976) yields the tight-binding parameter as Büchler et al.  $`t=U_0\mathrm{exp}\left[\frac{\pi }{2}\sqrt{U_0\left(\mathrm{}^2/ma^2\right)^1}\right],`$ which is much smaller than the first excitation energy of the orbital. The effect of the overall trapping potential on $`t`$ can be neglected provided that the potential energy difference between the minima of neighbouring wells, $`\mathrm{\Delta }E_j`$, is much less than the barrier height, $`U_0`$. Now, $`\mathrm{\Delta }E_j\kappa a^2j`$ for $`j1`$; hence our model is accurate for states whose wave functions do not extend beyond site $`j_cU_0/\kappa a^2`$. Hence the state’s energy $`ϵ`$ must satisfy $`ϵ\kappa a^2j_c^2`$, i.e. $`ϵϵ_{\mathrm{max}}U_0^2/\kappa a^2`$. We would therefore like to be sure that $`ϵ_{\mathrm{max}}t`$, which is true provided that $`U_0\kappa a^2`$.

We are interested in the density of states (DOS), $`\rho (ϵ)_\nu \delta \left(ϵϵ_\nu \right)`$, which is a smooth curve only in the limit in which the spectrum $`\left\{ϵ_\nu \right\}`$ is a continuum. Let us discuss some limiting behaviours of $`\rho (ϵ)`$. The solution is well known for $`\kappa =0`$, when the lattice forms a large ‘box’ of length $`L`$. In this case one obtains Bloch waves with the band dispersion relation $`ϵ_k=2t\mathrm{cos}ka,`$ where $`k`$ is the wave number. When $`La`$, $`ka`$ runs continuously from $`\pi `$ to $`\pi `$. In this limit, $`\rho (ϵ)`$ is non-zero only for $`|ϵ|2t`$ (this region we call the ‘band’), and has square-root singularities at $`ϵ=\pm 2t`$ (the ‘top’ and ‘bottom’ of the band). For $`d=2`$ the DOS extends from $`ϵ=4t`$ to $`4t`$ and is finite at these points, but exhibits a ‘van Hove’ logarithmic singularity at $`ϵ=0`$ Ashcroft and Mermin (1976).

In our model, the lattice itself is taken to be infinite; the finite extent of the wave functions is determined instead by the strength, $`\kappa `$, of the harmonic trap. The corresponding length scale is Damle et al. (1996) $`l\sqrt{t/\kappa }`$. This suggests that, in this model, ‘continuum limit’ should be taken to mean $`la`$ or, equivalently, $`\kappa a^2t`$. We shall work in this limit henceforth (unless otherwise stated).

Let us consider the low- and high-energy states of (1). For energies near the bottom of the band, the dispersion relation has the free-particle form $`ϵ_k2t+\mathrm{}^2k^2/2m^{}`$, where $`m^{}=\mathrm{}^2/2ta^2`$. Thus we would expect the low-lying eigenstates of (1) to resemble those of the usual continuum harmonic oscillator. Indeed, one can show explicitly that $`|\varphi _0=N_0_je^{aj^2/2\sqrt{2}l}|j`$ is the ground state to second order in $`a/l`$. Its energy is $`ϵ_0=2t+\mathrm{}\omega ^{}/2+o\left([a/l]^2\right),`$ where $`\omega ^{}\sqrt{\kappa /m^{}}`$. To first order in $`a/l`$, this equals the usual result for the harmonic oscillator, measured from $`ϵ=2t`$. The low-lying states are obtained similarly. Since the DOS of the 1D harmonic oscillator is constant and equal to $`1/\mathrm{}\omega ^{}`$, we conclude that the low-energy limit of the DOS of (1) is given by $`\rho \left(ϵ\right)1/\mathrm{}\omega ^{}`$ as $`ϵ2t`$ for $`d=1`$. For $`d=2`$, the DOS of the harmonic oscillator is $`\rho \left(ϵ\right)ϵ`$, so we expect a vanishing density of states as $`ϵ4t`$ in this case.

For high energies, the physical nature of the states is quite different. When $`ϵt`$, the first term in (1) can be neglected, and thus the eigenstates of the Hamiltonian are the position eigenstates $`|j`$, with energies $`ϵ_j=\kappa a^2j^2/2.`$ The high-energy 1D DOS is then easily shown to be $`\rho (ϵt)=\sqrt{2/\kappa a^2}ϵ^{1/2}`$. In 2D, the analogous calculation gives $`\rho (ϵ)2\pi /\kappa a^2`$ as $`ϵ/t\mathrm{}`$.

*WKB analysis.* To extend our treatment of (1) to all energies, we employ the Wentzel-Kramers-Brillouin (WKB) approximation Bender and Orszag (1999); it has already been applied to the momentum-space version of the problem in ref. Polkovnikov et al. (2002). By contrast, our WKB-type analysis will be performed in real space O’Dell (2001). It differs from traditional WKB in that $`\mathrm{}`$ does not appear anywhere in (1); instead we choose $`a`$ as our small parameter. As in the ordinary WKB method, we write the wave function in the form $`\psi (x)\mathrm{exp}\left(i\stackrel{x}{}k(x^{})𝑑x^{}\right)`$. The energy is determined by the usual quantisation condition, $`k(x)𝑑x=2\pi \left(n+\gamma \right)`$, where $`n`$ is an integer and $`\gamma `$ a constant. Differentiating this, we obtain the DOS: $`\rho (ϵ)dn/dϵ=\frac{1}{2\pi }\frac{}{ϵ}k(x)𝑑x`$ (note that $`\gamma `$ is not required).

On the other hand, the orbit equation is modified by replacing the free-particle dispersion by the tight-binding form Lifshits et al. (1973): $`2t\mathrm{cos}\left(ka\right)+\frac{1}{2}\kappa x^2=ϵ`$. The local wavenumber becomes $`k(x)=\pm \frac{1}{a}\mathrm{arccos}\left(\frac{\kappa x^22ϵ}{4t}\right)`$. This makes a significant difference, because it introduces two new turning points. The classical turning points (where $`k=0`$) still exist, and are found at $`x=\pm x_c`$, where

$$x_c=\sqrt{\frac{2ϵ+4t}{\kappa }}=2l\sqrt{\frac{ϵ}{2t}+1};$$
(2)
now, however, two new points appear for energies $`ϵ>2t`$. They are associated with Bragg reflection, and they occur at $`x=\pm x_b`$, where

$$x_b=\sqrt{\frac{2ϵ4t}{\kappa }}=2l\sqrt{\frac{ϵ}{2t}1}.$$
(3)
We shall refer to them as ‘Bragg turning points’, and the region between them as ‘Bragg-forbidden’.

The quantisation condition now becomes more subtle. To see why, note that in the Bragg-forbidden region $`k(x)a=\pi +i\mathrm{arccosh}\left(\frac{2ϵ\kappa x^2}{4t}\right)`$: the wave function continues to oscillate for $`|x|<x_b`$, despite the damping of its amplitude. The quantisation condition must therefore include these oscillations; it becomes:

$$2\underset{x_b}{\overset{x_c}{}}k(x)𝑑x+\frac{2\pi x_b}{a}=\pi \left(n+\gamma \right).$$
(4)
Now we can obtain the DOS by differentiating (4), remembering that $`x_b`$ and $`x_c`$ depend on $`ϵ`$ according to (2,3); we obtain a lengthy analytic expression for $`\rho (ϵ)`$. It is plotted in Fig. 3(a), along with the DOS obtained from numerical diagonalisation of (1), with appropriate provision for finite-size scaling (see below): clearly the agreement is excellent. This was foreseeable, since aside from a vanishing number of low-lying states, all the wave functions have slowly varying $`k(x)`$ (or, in the Bragg-forbidden region, not varying at all); thus the WKB method is expected to be accurate at all energies. Note that $`\rho (ϵ)`$ is qualitatively different from the DOS of the $`\kappa =0`$ tight-binding model, while it agrees quantitatively with $`\rho 1/\mathrm{}\omega ^{}`$ at $`ϵ=2t`$ and with $`\rho \sqrt{2/\kappa a^2}ϵ^{1/2}`$ for $`ϵ/t1`$, the limiting expressions obtained above.

Calculating the density of states for $`d>1`$ is simple, because the Hamiltonian is separable. Consider, for example, the case of $`d=2`$. It is not hard to show that $`\rho _{2\mathrm{D}}(ϵ)=_{\mathrm{}}^ϵ\rho (ϵ\chi )\rho (\chi )𝑑\chi `$, where $`\rho _{2\mathrm{D}}`$ is the DOS of the 2D system, and $`\rho `$ is the DOS of the 1D Hamiltonian (1). The result of applying this formula is presented in Fig. 3(b), and compared with numerical results; again, the agreement is excellent. Note that the steep feature at $`ϵ=0`$ and the kink at $`ϵ=4t`$ both result from the logarithmic singularity in the 1D DOS. As $`ϵ4t`$, we recover $`\rho _{2\mathrm{D}}0`$, as expected from the analytic arguments above. Furthermore, the large-$`ϵ`$ asymptote is exactly the constant predicted by that analysis. Higher-dimensional densities of states may be generated in a similar way.

*Finite-size scaling: numerics.* In the foregoing we obtained the DOS in the continuum limit $`l/a\mathrm{}`$. However in experiments $`l/a`$ is finite. We should therefore check that our results resemble the behaviour of the system when it has a steeper trap. Thus we have studied the single-particle spectrum and wave functions by numerical diagonalisation of (1) 111We used the diagonalisation function of Eaton et al. .. By considering increasingly shallow traps, the numerical calculations also allow us to confirm the validity of WKB for the continuum limit. We have applied the same method in $`d=2`$, starting with the 2D version of (1). The results are analogous so we shall, as above, describe in detail the 1D case only.

In order to have a finite-dimensional Hamiltonian matrix, we must describe the lattice by a finite number of sites, $`N`$. This leads to numerical artefacts for $`ϵ>ϵ_N`$, defined by $`N2x_c(ϵ_N)/a`$, as the particle visits the artificial limits of the lattice. On the other hand we expect this to be essentially equivalent to the original problem for $`ϵϵ_N`$. Calculations for increasingly large $`N`$ confirm this (see the inset of Fig. 2).

Fig. 1 shows some numerically determined wave functions, for fixed trap strength. In spite of the relatively small system size ($`l/a=10`$), the states seem well described by the continuum-limit solution. At low energies, $`ϵ0`$, the wave functions resemble those of a continuum 1D harmonic oscillator. At higher energies, the period of these oscillations starts to become similar to the lattice constant, and commensuration effects emerge. At $`ϵ=2t`$, as that period reaches its minimum, the Bragg reflection points appear. Immediately above the top of the band, $`ϵ2t`$, the modulation of the wave number is evident, reaching its minimum and maximum near the classical and Bragg turning points, respectively. The rapid oscillations continue in the Bragg-forbidden region, where the amplitude decays exponentially. As the energy rises further, the distance between the Bragg and classical turning points shrinks, so that fewer and fewer oscillations take place between these two points. Eventually the particle is forced to localise on single sites (not shown). Note that the Bragg and classical turning points are very accurately described by Eqs. (2,3), even for the relatively small value of $`l/a`$ considered here. Furthermore, the number of high-energy (localised) solutions is reproduced correctly. This is because the condition that the wave function’s phase be single-valued reduces, at $`k=\pi `$, to the condition that the localised states be separated by an integer number of lattice spacings ($`x=na`$).

We now turn to the DOS. Obviously for finite $`l/a`$ the spectrum is not a continuum so the DOS is not a smooth curve. However, one can smooth it by the following procedure: we divide an energy range into intervals, of width $`\mathrm{\Delta }ϵ`$, and count the states within each interval. We then re-scale the vertical axis of the resulting histogram so that each column represents the number of states per unit energy, on average, in the corresponding interval. The result is a ‘binned’ DOS that looks smooth only if the intervals are sufficiently wide. Fig. 2 shows the result for different values of $`\kappa a^2/t`$. Note that we have further re-scaled the vertical axis by an overall factor of $`a/l`$. This leads to the collapse of all the data onto a single curve, suggesting that the continuum-limit DOS describes the overall distribution of energy levels even for quite steep traps, with only $`l/a100\text{ sites}`$. The inset illustrates the numerical cutoff artefact mentioned above. Finally, Fig. 3 compares the DOS of two finite (but fairly large) systems to the WKB predictions obtained above in 1D and 2D respectively. Evidently our WKB approach captures the continuum limit rather well.

*Conclusion.* We have calculated the DOS for a single atom in an optical lattice; this should be regarded as the logical first step towards a detailed theory of the experimentally realised many-particle systems. Our results are based on WKB theory, and refer to the limit $`la`$ or, equivalently, when the trapping potential becomes flat: $`\kappa 0`$. Numerical diagonalisation reveals this theory to be extremely accurate in that case, and moreover shows that, for finite-size systems, the binned DOS has the same overall features. Our main result is that the DOS, in this limit, is radically different from what is obtained for a homogeneous lattice (i.e. for $`\kappa =0`$ rather than $`\kappa 0`$): the square-root singularities in the 1D case are replaced by a logarithmic one, and the logarithmic van Hove singularity in 2D disappears altogether. Moreover, our theory provides a detailed picture of how this comes about. The crucial feature is the new turning points, associated with Bragg reflection, that appear at energies above the top of the conduction band. The possibility of inducing Bragg reflection by using a time-dependent external potential was considered in Scott et al. (2003). We have shown that Bragg reflection is, in fact, essential to understand the equilibrium single-particle spectrum of the optical lattice.

###### Acknowledgements.

We thank J. M. F. Gunn, A. F. Ho, M. W. Long and A. J. Schofield for useful discussions. We acknowledge financial support from EPSRC (CH) and the Leverhulme Trust (JQ). Note added.— It has been drawn to our attention that some of the features discussed here have recently been pointed out in an exact diagonalisation study Rigol and Muramatsu .


text:
# On Using Expansions to the Base of -2

## 1 Introduction

Recently, groups associated to elliptic and hyperelliptic curves received a lot of attention for cryptographic applications, and further kinds of curves were proposed and their arithmetic studied intensively. They allow smaller operands compared to RSA and DL in finite fields making them attractive for restricted devices. The performance on such small units is good . More general curves were suggested and the group operations optimized for cryptographic applications.

To compute scalar multiples $`mD`$, binary expansions of $`m`$ are used and the computation of $`mD`$ is split up as a sequence of additions and doublings. To achieve faster computations one uses windowing methods and signed representations (for a broad overview see Knuth ).

Our idea speeds up scalar multiplication in groups for which computing $`2D`$ and $`(D_1+D_2)`$ is faster than computing $`2D`$ and $`D_1+D_2`$, respectively. For elliptic and hyperelliptic curves the negative of an element can be obtained almost for free. Hence, these groups will most probably not benefit from our new idea. But, the situation is different for Picard curves or more general genus 3 quartic curves: using the $`2`$-adic expansion instead of a $`2`$-adic we reduce the complexity by at least 7%.

Of course, our considerations are not restricted to cryptography but allow speeding up scalar multiplication in groups in which “adding up to the neutral element” is easier computed than addition and so they could be of use in computer algebra systems, too.

In this note we first describe the idea of $`2`$-adic expansions and show how to apply them. Then we sketch the applications we have in mind, and finally show the time saving for Picard curves.

To fix notation, let $`G`$ be a finite abelian group of order $`\mathrm{}`$ and let $`D`$ be a generator of $`G`$. Furthermore, we assume that computing $`2D`$ or $`(D_1+D_2)`$ is actually faster than computing $`2D`$ or $`D_1+D_2`$ respectively.

## 2 $`2`$-adic Expansions

Assume that we want to compute $`mD,m<\mathrm{}`$ and put $`l(m)`$ the length and $`w(m)`$ the number of nonzero bits of the used expansion of $`m`$. Since we need both $`D`$ and $`D`$, we can allow signed digit representations.

If $`l(m)+w(m)`$ is even we start with $`D`$, otherwise with $`D`$. While the doubling is always replaced by the computation of $`2`$-times the intermediate result $`E`$ we need to pay a little more attention on how to perform the former additions and subtractions.

###### Algorithm 1

INPUT: $`m\mathrm{IN}`$, $`m=_{i=0}^{l(m)1}m_i2^i,m_i\{0,\pm 1\}`$, $`DG`$
OUTPUT: $`E:=mD`$

1. precompute and store $`D`$
2. compute $`l(m),w(m)`$;
3. put $`E:=(1)^fD`$, where $`f:=l(m)+w(m)mod2`$;
4. for $`i=l(m)2`$ to $`0`$ do

1. $`E:=2(E)`$;
    2. $`f:=1f`$;
    3. if $`m_i0`$

    1. $`E:=(E+(1)^fm_iD)`$;
        2. $`f:=1f`$;
5. output($`E`$);

The correctness follows from the fact that $`l(m)+w(m)2`$ is the total number of minus signs in front of the initial $`E`$. During the process $`f`$ keeps track of the parity of the number of sign changes. In step 4(c)i, no multiplication is required. $`f`$ assumes only values in $`\{0,1\}`$ and $`m_i\{\pm 1\}`$.

Using this idea introduces only little bookkeeping overhead, namely the additional variable $`f`$. So such a system is really practicable – and useful if the operations involving the negative signs are faster.

Remarks:

1. If the expansion of $`m`$ is not computed beforehand, one can always start with $`E=D,f=0`$, irrespective of the parity of $`l(m)+w(m)`$. The loop 4. is performed as above. Before Step 5. one checks whether $`f=1`$, in which case one outputs $`E`$ instead. So one avoids precomputing the expansion, at the price of a second negation with probability 1/2.
2. If only for one of addition or doubling the negative is faster, similar considerations hold if one only replaces that operation.
3. Of course the method can be combined with signed sliding windowing methods. In the applications we have in mind $`rD`$ and $`rD`$ can be computed with only a few more operations than $`rD`$ alone.

## 3 Applications

In this section we need to state some details from mathematics to show that there actually are applications of our idea. For an introduction to hyperelliptic curves see . The following holds for arbitrary curves.

Let $`C`$ be a curve of genus $`g`$. The group used for cryptographic applications is a subgroup of the divisor class group of $`C`$: we briefly recall its main properties. Let $`P_{\mathrm{}}C`$ be fixed. A divisor is a formal sum of points. We are interested in the degree zero divisors given by sums

$$DnP_{\mathrm{}}=\underset{i=0}{\overset{n}{}}P_inP_{\mathrm{}},P_iC\{P_{\mathrm{}}\}.$$
(1)
The principal divisors are the divisors of functions. The divisor class group is the group of the degree zero divisors modulo the principal ones. In each divisor class there exists a unique element (1) with $`ng`$ minimal.

To add two classes $`c_1,c_2`$ one formally adds the representing divisors: $`D_1+D_2(n_1+n_2)P_{\mathrm{}}`$. Then one determines a function $`f`$ passing through the points on $`D_1+D_2`$ with poles only in $`kP_{\mathrm{}}`$ ($`k`$ minimal) with multiplicities taken into account. Let $`D_3`$ be the divisor represented by the points of intersection of $`f`$ with the curve which are not in $`D_1`$ and $`D_2`$. Put $`n_3=kn_1n_2`$ and let $`c_3`$ be the class of $`D_3n_3P_{\mathrm{}}`$. Since $`D_1+D_2+D_3kP_{\mathrm{}}`$ add up to a principal divisor we get $`c_1+c_2+c_3=0`$.

Usually one proceeds to find the negation $`c_4=c_3`$ to get $`c_1+c_2=c_4`$. Our new proposal allows to skip this last step. Doublings just work the same with the function passing through the points of the representing divisor with doubled multiplicity.

For hyperelliptic curves, taking the negative is very simple. The formulae for genus $`2`$ and $`3`$  reveal that computing $`2D`$ instead of $`2D`$ saves only some additions in the underlying field. Therefore, we do not expect the $`2`$-adic expansion to lead to a saving.

But the situation is completely different for non-hyperelliptic curves. For char$`(\mathrm{IF}_q)3`$ a Picard curve can be given by:

$$zy^3=z^4f_4(x/z),f_4\mathrm{IF}_q[x]$$

where $`f_4`$ is monic, square-free and of degree 4.

The arithmetic on Picard curves is detailed in  (see also ). An addition needs 144M, 12S, and 2I and a doubling 158M, 16S, and 2I in the generic case. Applying Algorithm 1 reduces the costs to 133M, 9S, 2I or 147M, 13S, 2I respectively. Some field additions are saved as well. Thus, here the saving is at least 7.5% or 7%, respectively, assuming a ratio of $`10:1`$ for inversions and $`2:3`$ for squarings in relation to multiplications.

An ordinary genus 3 quartic over $`\mathrm{IF}_q`$ is given by a projective equation:

$$\left(a_1x^2+a_2y^2+a_3z^2+a_4xy+a_5xz+a_6yz\right)^2=l(x,y,z)xyz,$$

where $`l\mathrm{IF}_q[x,y,z]`$ is linear and $`a_i\mathrm{IF}_q`$. On these curves computing the negation is even more complicated than on Picard curves (see ). Therefore, the saving due to $`2`$-adic expansions is more dramatic.


text:
# Kondo Effect and Surface-State Electrons

## Abstract

We have used low temperature scanning tunneling spectroscopy and atomic manipulation to study the role of surface-state electrons in the Kondo effect of an isolated cobalt atom adsorbed on Ag(111). We show that the observed Kondo signature remains unchanged in close proximity of a monoatomic step, where the local density of states of the surface-state electrons is strongly perturbed. This result indicates a minor role for surface-state electrons in the Kondo effect of cobalt, compared to bulk electrons. A possible explanation for our findings is presented.

Kondo effect, surface-state electrons, scanning tunneling spectroscopy

Recently, scanning tunneling microscopy (STM) has been used to investigate the electronic properties of single magnetic impurities on the surface of metals, opening a new avenue of research into the Kondo effect Li ; Madhavan . The Kondo effect occurs when a magnetic impurity in a metal couples to surrounding conduction electrons. This coupling, which involves spin-flip scattering events on the site of the impurity, causes conduction electrons to form a correlated ground state that screens the magnetic moment of the impurity at temperatures below a characteristic Kondo temperature ($`T_K`$) Kondo HF . It also leads to a narrow resonance in the local density of states (LDOS) of the impurity at the Fermi energy ($`E_F`$), which is detected in STM and, in particular, in scanning tunneling spectroscopy (STS).

Owing to the atomic resolution in space and the better than meV resolution in energy, STM and STS can perform local investigations of the Kondo effect, directly at the adsorption site of single magnetic atoms on the surface of metals. However, for a comparison to bulk measurements, differences in the bulk and in the surface electronic structures must be taken into account. For instance, well known bulk Kondo systems yield Kondo temperatures substantially higher than the ones reported by STM Madhavan ; Manoharan ; Knorr . Furthermore, most STM and STS measurements have been performed on (111) facets of noble metals, where bulk electrons co-exist with Shockley surface-state electrons, which form a quasi-two-dimensional electron gas trapped between the surface barrier potential and a band gap in the crystal SurfState .

So far, there still is a lack of experimental and theoretical data concerning the role of the surface-state electrons in the Kondo effect. In this brief report we specifically focus on this topic. We report a STM and a STS study for cobalt adatoms on Ag(111), where the combination of controlled cobalt manipulation over the surface and tunneling-spectroscopy of cobalt adatoms enables us to tune the LDOS of the surface-state electrons and thus probe its influence on the Kondo temperature of Co/Ag(111). We also provide a possible scenario to interpret our data, which, we hope, will trigger further theoretical studies.

The measurements were performed in a home-built ultrahigh vacuum STM at a working temperature of $`T=4.6`$ K and the Ag(111) surface was cleaned by Ar<sup>+</sup> sputter/anneal cycles. Single cobalt atoms were evaporated onto the cold Ag substrate by heating a degassed cobalt wire ($`>99.99\%`$) wound around a pure W wire ($`>99.95\%`$). The evaporation, through an opening of the He shield of the cryostat, yielded a coverage of $`3\times 10^3`$ ML. No appreciable increase of other impurities was detected. After the evaporation, the thermalization of Ag(111) to $`4.6`$ K ensured a negligible thermal diffusion of cobalt adatoms within measurement times. Spectra of the differential conductance of the tunneling junction, $`dI/dV(V)`$, were recorded via lock-in detection (AC modulation was $`1`$ mV (rms) in amplitude and $`10`$ kHz in frequency), where $`V`$ is the sample bias measured with respect to the tip. The images and the spectra were recorded with a etched tungsten tip; the tip was further treated in situ by soft indentations into the Ag(111) surface, until the cobalt adatoms were imaged spherically (Fig. 1), and the $`dI/dV`$ had no structure near the zero bias voltage, i.e. the Fermi energy. The cobalt adatoms have then a Gaussian-like profile $`0.7`$ Å high with an apparent diameter of $`6`$ Å (full width at half maximum).

A typical $`dI/dV`$ spectrum acquired at the center of a cobalt adatom is presented in Fig. 2a. The dip in the LDOS near $`E_F`$ is the signature in STS of the Kondo effect of Co/Ag(111). When the tip is moved laterally off the cobalt adatom, the dip continuously decreases in amplitude and vanishes at $`8`$ Å away from the center of the cobalt atom; at a distance of $`50`$ Å, the $`dI/dV`$ spectrum of the Ag(111) surface-state is recovered (Fig. 2b). A similar dip at $`E_F`$ has previously been observed in STS for other Kondo systems like Ce/Ag(111) and Co/Cu(111) Li ; Manoharan ; Knorr , in contrast to the asymmetric line shape observed for Co/Au(111), Ti/Ag(100) and Co/Cu(100) Madhavan ; Nagaoka ; Knorr . The $`dI/dV`$ spectrum can be understood following the pioneering framework by Fano Fano . When $`TT_K`$, spin-flip processes are frozen out and the Kondo many-body problem reduces to a single-particle resonance near $`E_F`$ of width $`2k_BT_K`$ Noziers ; Kondo HF , which, in principle, should be detected in a $`dI/dV`$ spectrum. However, tunneling from the tip into the Kondo resonance occurs through two channels: a direct one, and an indirect one via the $`sp`$ substrate conduction band locally hybridized with the cobalt $`d`$ levels (which are giving rise to the magnetic moment of the cobalt atom). The interference between the two channels $``$ known as Fano interference $``$ produces a symmetric or an asymmetric line in the $`dI/dV`$ spectrum, depending on the relative weight of the two channels. Modelling of the dip of Fig. 2a with a Fano line Ujsaghy ; Ellipse2 , allows then to extract from the $`dI/dV`$ spectrum information concerning the Kondo effect for Co/Ag(111). The Fano fits on various cobalt spectra yielded the following results: from the line width we extract a temperature $`T_K=83(10)`$ K higher than our working temperature of 4.6 K, a Fano parameter $`q=0.0(1)`$ $``$ indicating that the predominant tunneling channel is through the $`sp`$ band $``$ and a shift of the resonance to $`5.8(4)`$ meV above $`E_F`$. The reported values are in good agreement with a previous STS study of Co/Ag(111), where a temperature of $`T_K=92(6)`$ K was found by Schneider et al. Schneider .

Figures 3 and 4 illustrate the main experimental findings of this report. To explore the role played by the Ag(111) surface-state electrons in the Kondo effect of Co/Ag(111), we combined manipulation and spectroscopy experiments to probe how the Kondo effect of Co/Ag(111) is affected in close proximity of a monoatomic step of Ag(111) where the surface-state LDOS ($`\rho _s`$) is strongly perturbed. Figure 3 illustrates a typical manipulation procedure we performed. By placing the STM tip over a cobalt adatom and by applying $`I=2.8`$ nA and $`V=2.8`$ V to the tunneling junction (Fig. 3a), the cobalt adatom is trapped below the tip and gently dragged to a desired location on the surface $``$ in Fig. 3b, at $`8`$ Å from a monoatomic step of Ag(111). Figure 4a presents $`dI/dV`$ data acquired over a cobalt adatom positioned below the step of Fig. 3, at distances $`r`$ from the step edge of $`8`$, $`13`$ and $`24`$ Å. The Fano fits to the three spectra do not reveal any appreciable change in the resonance, indicating that $`T_K`$ is constant with $`r`$. Near a step, however, $`\rho _s`$ is strongly perturbed at $`E_F`$ because of the back-scattering of the surface state by the step $``$ contrary to the bulk LDOS ($`\rho _b`$). Neglecting thermal broadening, the spatial variation of $`\rho _s(E_F)`$ is approximately described by a Bessel function $`1J_0(2k_Fr)`$ SurfState ; LDOS step , where $`k_F=0.084`$ $`\text{Å}^1`$ for Ag(111). The dependency of $`\rho _s(E_F)`$ on $`r`$ is presented in Fig. 4b where we have also indicated the positions where the spectra of Fig. 4a were recorded. As shown, when moving from $`r=24`$ Å to $`r=8`$ Å, $`\rho _s(E_F)`$ decreases by more then a factor three, to be compared to the constant value of $`T_K`$ observed for Co/Ag(111) at a Ag(111) step. The only significant spectroscopic change we observe in the Kondo effect of Co/Ag(111) is for artificially fabricated cobalt dimers, where, in agreement with the findings of Ref. Chen, , an abrupt disappearance of the Kondo resonance occurs.

To summarize the experimental data, we have studied the Kondo effect of cobalt atoms adsorbed on the Ag(111) surface, which also exhibits surface-state electrons. The STS data acquired near a Ag(111) step indicate that the Kondo temperature of Co/Ag(111) is robust on substantial changes in the LDOS of the surface-state electrons, suggesting a minor role for these electrons in this Kondo system.

To discuss these findings, we recall some fundamental aspects of the Kondo problem. The Kondo effect is essentially determined by one energy scale: $`k_BT_K`$. Within the Coqblin-Schrieffer model Kondo HF , the Kondo temperature is

$$T_K=De^{1/2J\rho },$$
(1)
where $`J`$ is the antiferromagnetic coupling strength between conduction electrons and the magnetic moment of a cobalt atom, and $`\rho `$ the LDOS at $`E_F`$ ($`D`$ is the bandwidth of silver). Although a quantitative interpretation of our experimental findings would need a microscopic theory describing the Kondo scattering of bulk and of surface-state electrons, which is beyond the scope of this report, we may draw, based on Eq. (1), some conclusions concerning the involvement of surface-state electrons in the Kondo effect.

First, as seen in Eq. (1), the Kondo temperature of a magnetic scatterer depends exponentially on $`J\rho `$. Since experimentally $`T_K`$ is found to be insensitive to a variation of $`\rho _s`$ by a factor of three, we conclude that the Kondo effect of cobalt atoms adsorbed on Ag(111) involves mainly bulk electrons: $`J\rho J_b\rho _b`$. This conclusion agrees with those of Ref. Knorr, , where the Kondo temperature of the Co/Cu(111) and Co/Cu(100) systems was shown to vary with the number of next nearest Cu neighbors ($`n`$) of the cobalt atom, hence with $`\rho _bn`$.

Next, we focus more specifically on the surface-state electron coupling $`J_s`$. At the tip-surface distances where usual $`dI/dV`$ spectra are acquired ($`10`$ Å over the surface), about 2/3 of the current into the Ag(111) surface is due to tunneling into the surface-state LDOS (see Fig. 2b), and in particular $`\rho _s/\rho _b2`$ at $`E_F`$. Assuming that this ratio between the two LDOS holds also at the cobalt site, i.e. that $`\rho _s\rho _b`$ at the impurity site, the feeble contribution of the surface-state electrons to the Kondo effect is to be found then in a weaker coupling to the magnetic moment of the cobalt atom compared to bulk electrons: $`J_sJ_b`$. While $`J_s`$ is small, it is clearly non-zero, as demonstrated by Manoharan et al. in their quantum mirage experiment Manoharan , where they exploited the scattering of the Cu(111) surface-state electrons in an elliptical quantum corral to project the Fano line of a cobalt adatom in one focus to the other empty focus, about $`80`$ Å away. However, only ellipses of specific dimensions will have a sufficient surface-state electron amplitude at the focal adatom to yield a detectable Kondo signal at the opposite focus Ellipse ; Ellipse2 , in agreement with our finding of a minor involvement of surface-state electrons in the surface Kondo effect.

In conclusion, we have studied the Kondo effect of cobalt atoms adsorbed on the surface of Ag(111) by STM and by STS. The Kondo signature in the $`dI/dV`$ of Co/Ag(111) is a dipped Fano resonance near $`E_F`$ from which is extracted a Kondo temperature of $`T_K=83(10)`$ K. A combination of cobalt-manipulation and tunneling spectroscopy, suggests, on an experimental basis, that the role of surface-state electrons in the Kondo effect of Co/Ag(111) is minor, possibly because of a weaker coupling $`J_s`$ compared to the one of the bulk electrons. This appears to be a general property of Kondo systems and, we hope, will trigger further theoretical studies for a quantitative understanding of our experimental data.

We gratefully acknowledge A. Schiller and R. Bulla for fruitful discussions, and thank the Deutsche Forschungsgemeinschaft for financial support.


text:
# Quantum and classical features in the explanation of collisional decoherence

## I Introduction

Paradoxical as it may seem, one of the most intriguing features of the quantum world is certainly its relationship to the classical one KieferNEW . Apart from the conceptual charm this question is of practical relevance in many respects. Full understanding of this issue would give an answer to which quantum phenomena might appear or persist on a meso- or macroscopic scale, to the degree of feasibility of quantum computation, to the reliance of classical pictures and classical computations, often simpler and more intuitive. First significant steps in this direction are certainly given by recent experiments in which loss of coherence due to the interaction with the environment can be quantitatively studied. This phenomenon goes under the name of decoherence, and its theoretical description is based on an understanding of open system dynamics Petruccione .

In this brief paper we will focus on recent beautiful experiments concerned with the observation of collisional decoherence in matter-wave interferometry ZeilingerQBM-exp , sketchily giving a most simple theoretical argument for the loss of visibility, thus showing the connection with the classical linear Boltzmann equation and further giving an exact expression for the macroscopic scattering cross-section.

## II Analysis of collisional decoherence

In the aforementioned experiments the visibility of the interference fringes, obtained letting a beam of fullerenes go through a Talbot-Laue interferometer, is progressively reduced by raising the pressure of the background gas present in the experimental apparatus. The visibility decreases according to the formula

$$V=V_0\mathrm{exp}(\mathrm{\Gamma }t)$$
(1)
where $`V_0`$ is the reference value, $`t`$ the time of flight through the apparatus, $`\mathrm{\Gamma }`$ a decoherence rate given by

$$\mathrm{\Gamma }=nv_0\sigma _{\mathrm{macro}}(v_0)$$
(2)
with $`n`$ the density of the gas, $`M`$ the mass of the fullerenes, $`v_0=p_0/M`$ the modulus of their incoming velocity, supposed to be essentially unchanged despite the collisions with the background particles, $`\sigma _{\mathrm{macro}}(v_0)`$ the total macroscopic scattering cross-section off the background gas, parametrically depending on the modulus of the incoming velocity or equivalently on momentum. To give a theoretical justification of the formula (1) besides a detailed characterization of the experimental apparatus (in this connection see also garda03post ) one essentially needs a quantum mechanical description of the fullerene dynamics due to collisions with a surrounding gas. This is achieved by means of a quantum master-equation which, as the underlying physics in terms of collisions obviously tells, is actually a quantum counterpart of the classical linear Boltzmann equation. The history of this master-equation is relatively long, and we refer the reader to vienna  for a discussion of this point. A general result in this direction, obviously still open to improvements is given in art7 , to which we refer for further detail and references. The master-equation takes the form

$`{\displaystyle \frac{d\widehat{\rho }}{dt}}={\displaystyle \frac{i}{\mathrm{}}}[\widehat{𝖧}_0,\widehat{\rho }]+{\displaystyle \frac{n}{M^2}}{\displaystyle d^3𝒒\mathrm{\Sigma }(q)}`$ (3)
$`\times \left[e^{\frac{i}{\mathrm{}}𝒒\widehat{𝗑}}\sqrt{S(𝒒,\widehat{𝗉})}\widehat{\rho }\sqrt{S(𝒒,\widehat{𝗉})}e^{\frac{i}{\mathrm{}}𝒒\widehat{𝗑}}{\displaystyle \frac{1}{2}}\{S(𝒒,\widehat{𝗉}),\widehat{\rho }\}\right],`$
with $`\widehat{𝗑}`$ and $`\widehat{𝗉}`$ position and momentum operator for the test particle (presently the fullerene), $`\widehat{𝖧}_0=\widehat{𝗉}^2/2M`$, $`\mathrm{\Sigma }(q)=M^2(2\pi \mathrm{})^3\frac{2\pi }{\mathrm{}}|\stackrel{~}{t}(q)|^2`$ ($`\stackrel{~}{t}(q)`$ being the Fourier transform of the interaction potential), $`S(𝒒,\widehat{𝗉})`$ the operator-valued dynamic structure factor, given by the Fourier transform of the time dependent spatial autocorrelation function of the gas

$`S(𝒒,E(𝒒,𝒑))`$ $`=`$ $`{\displaystyle \frac{1}{2\pi \mathrm{}}}{\displaystyle 𝑑td^3𝒙e^{\frac{i}{\mathrm{}}(E(𝒒,𝒑)t𝒒𝒙)}}`$
$`\times {\displaystyle \frac{1}{N}}{\displaystyle d^3𝒚N(𝒚)N(𝒙+𝒚,t)}.`$
with $`𝒒`$ and $`EE(𝒒,𝒑)=\frac{q^2}{2M}+\frac{𝒑𝒒}{M}`$ momentum and energy transferred to the fullerene molecule in each single collision. The two-point correlation function $`S(𝒒,E)`$ is related as shown by van Hove vanHove  to the macroscopic differential scattering cross-section of a test particle off a medium through the formula

$$\frac{d^2\sigma _{\mathrm{macro}}}{d\mathrm{\Omega }dE}=\frac{p^{}}{p}\mathrm{\Sigma }(q)S(𝒒,E).$$
(5)
Eq. (3) correctly avoids the extra $`2\pi `$ factor appearing for example in the well-known result of Gallis and Fleming Gallis90 , as stressed in ZeilingerQBM-th  and also found in HalliwellQBM , in both of which a particular example of dynamic structure factor also implicitly appears. The redundancy of this $`2\pi `$ factor can be best understood in physical terms as in Petruccione , where the authors already stressed that this factor would lead to a decoherence rate given by $`2\pi `$ times the total scattering rate. Experiments do in fact rule out this factor. With respect to other results the quantum linear Boltzmann equation (3) encompasses the momentum of the test particle as a dynamical variable, thus allowing for a correct description of energy transfer and for the expected thermal stationary solution art5 . Since on the time scale of the experiment thermalization does not play a role, as a first simplification the operator $`\widehat{𝗉}`$ in $`S(𝒒,\widehat{𝗉})`$ can be replaced by $`𝒑_0`$, the incoming momentum of the fullerene. Eq. (3) thus becomes, also exploiting (5)

$$\frac{d\widehat{\rho }}{dt}=\frac{i}{\mathrm{}}[\widehat{𝖧}_0,\widehat{\rho }]+\frac{n}{M^2}d^3𝒒\mathrm{\Sigma }(q)S(𝒒,𝒑_0)e^{\frac{i}{\mathrm{}}𝒒\widehat{𝗑}}\widehat{\rho }e^{\frac{i}{\mathrm{}}𝒒\widehat{𝗑}}n\frac{p_0}{M}\sigma _{\mathrm{macro}}(p_0)\widehat{\rho }.$$
(6)
Note that the dependence on the incoming momentum of the test particle is still crucial in order to explain the velocity dependence in the macroscopic scattering cross-section appearing in (2). The quantum nature of (6), apart from the free evolution responsible for the interference phenomena, only appears in the second term at the r.h.s., the loss term simply being a multiplicative factor as in the classical linear Boltzmann equation, apart from the obvious fact that the microscopic scattering cross-section describing the single collision is based on a quantum or at least a semi-classical calculation Williams . The second term at the r.h.s. of (6) however does not contribute to the interference fringes, because of its rapid oscillations leading to a destruction of the off-diagonal matrix elements as described e.g. in AlickiPRA ; Petruccione . Therefore eq. (6) predicts an exponential decay of the intensity relevant for the interference fringes according to the rate $`n\frac{p_0}{M}\sigma _{\mathrm{macro}}(p_0)`$, which can be simply calculated from the loss term in the classical linear Boltzmann equation. Following Williams , and expressing the scattering cross-section as a function of the incoming velocity, we can write

$$nv_0\sigma _{\mathrm{macro}}(𝒗_0)=nv_0d^3uf_{\mathrm{MB}}(𝒖)\frac{|𝒗_0𝒖|}{v_0}\sigma _{\mathrm{micro}}(𝒗_0𝒖),$$
(7)
where $`f_{\mathrm{MB}}(𝒖)`$ is the thermal Maxwell-Boltzmann distribution of velocities in the background gas. For a microscopic scattering cross-section of the form

$$\sigma _{\mathrm{micro}}(𝒗_0𝒖)=K|𝒗_0𝒖|^\alpha $$
(8)
the integral in (7) can be exactly calculated. Introducing the most probable velocity $`v_{\mathrm{mp}}=\sqrt{2/\beta m}`$ of gas particles with temperature $`1/\beta `$ and mass $`m`$ one has

$`\sigma _{\mathrm{macro}}(v_0)`$ $`={\displaystyle \frac{K}{v_0}}{\displaystyle d^3uu^{\alpha +1}f_{\mathrm{MB}}(𝒖+𝒗_0)}`$ (9)
$`=K{\displaystyle \frac{2}{\sqrt{\pi }}}{\displaystyle \frac{v_{\mathrm{mp}}^{\alpha +2}}{v_0}}e^{\left(\frac{v_0}{v_{\mathrm{mp}}}\right)^2}{\displaystyle _0^+\mathrm{}}𝑑tt^{\alpha +2}e^{t^2}\mathrm{sinh}\left(2t{\displaystyle \frac{v_0}{v_{\mathrm{mp}}}}\right)`$
and exploiting the result Ryzhik

$$_0^+\mathrm{}𝑑xx^{2\mu 1}e^{x^2}\mathrm{sinh}(\gamma x)=\frac{\gamma }{2}\mathrm{\Gamma }\left(\mu +\frac{1}{2}\right)e^{\frac{1}{4}\gamma ^2}\mathrm{\Phi }(1\mu ,\frac{3}{2};\frac{1}{4}\gamma ^2)$$

valid for $`\mu >\frac{1}{2}`$, with $`\mathrm{\Phi }(\alpha ,\beta ;\gamma )`$ the confluent hypergeometric function, one immediately obtains the exact expression

$$\sigma _{\mathrm{macro}}(v_0)=K\frac{2}{\sqrt{\pi }}\mathrm{\Gamma }\left(\frac{\alpha }{2}+2\right)\frac{v_{\mathrm{mp}}^{\alpha +1}}{v_0}\mathrm{\Phi }(\left(\frac{\alpha }{2}+\frac{1}{2}\right),\frac{3}{2};\left(\frac{v_0}{v_{\mathrm{mp}}}\right)^2).$$
(10)
For the case of interest $`\sigma _{\mathrm{micro}}`$, calculated on the basis of an isotropic potential $`U(r)=\frac{C_6}{r^6}`$, corresponds to the particular choice ZeilingerQBM-exp ; Maitland

$$K=\frac{\left(\pi ^6\frac{3}{8}\right)^{\frac{2}{5}}}{\mathrm{sin}\left(\frac{\pi }{5}\right)\mathrm{\Gamma }\left(\frac{2}{5}\right)}\left(\frac{C_6}{\mathrm{}}\right)^{\frac{2}{5}}\alpha =\frac{2}{5}$$

thus recovering, when considering only the first contributions in the expansion in terms of the small parameter $`v_0/v_{\mathrm{mp}}`$, the result given in ZeilingerQBM-th1

$$\sigma _{\mathrm{macro}}(p_0)=\frac{\left(\pi ^6\frac{3}{8}\right)^{\frac{2}{5}}}{\mathrm{sin}\left(\frac{\pi }{5}\right)\mathrm{\Gamma }\left(\frac{2}{5}\right)}\left(\frac{C_6}{\mathrm{}}\right)^{\frac{2}{5}}\frac{2}{\sqrt{\pi }}\mathrm{\Gamma }\left(\frac{9}{5}\right)\frac{v_{\mathrm{mp}}^{\frac{3}{5}}}{v_0}\left[1+\frac{1}{5}\left(\frac{v_0}{v_{\mathrm{mp}}}\right)^2+𝒪\left(\left(\frac{v_0}{v_{\mathrm{mp}}}\right)^4\right)\right].$$
(11)
A full theoretical explanation of the observed data, also keeping the detailed features of the interferometer into account, has been hinted in ZeilingerQBM-exp ; ZeilingerQBM-th1  and will appear in a further publication by the same group ZeilingerQBM-th2 . The argument presented here, however, gives a remarkably simple explanation of the experimental results by impinging on both quantum and classical aspects. On the one hand one needs a true quantum linear Boltzmann equation in order to understand how the gain term in the equation does not contribute to the interference pattern, due to a suppression of off-diagonal terms, on the other hand the loss term which quantitatively describes the reduction of visibility when treated classically leads to the formula (7) used to fit the experimental data. The main difference between the approach presented here and the one outlined in ZeilingerQBM-exp ; ZeilingerQBM-th1  lies in the fact that here we start with a quantum linear Boltzmann equation already containing, through the dynamic structure factor and the dependence on the momentum of the incoming particle, all the informations about the macroscopic scattering cross-section relevant for the loss of visibility. When considering a master-equation like the one by Gallis and Fleming, instead, the macroscopic or effective scattering cross-section depending on the velocity of the incoming fullerene has to be introduced as an additional information.

###### Acknowledgements.

The author would like to thank Prof. L. Lanz and Dr. K. Hornberger for useful discussions. The work was supported by MIUR under Cofinanziamento and FIRB.


text: # Untitled Document