text:
# Relativistic Unitarized Quark/Meson Model in Momentum Space

## 1 INTRODUCTION

Mesonic $`q\overline{q}`$ states are the simplest color-singlet configurations in QCD. However, physical mesons appear to be much more complicated objects, in view of the quite disparate spectra of especially the light scalars and pseudoscalars, and the large variety of states ranging from OZI-stable quarkonia to extremely broad, often badly established resonances. Because of these complexities, theoretical approaches usually focus on parts of the meson data only. Thus, one has quarkonium potential models (see e.g. Ref. ), relativistic constituent quark models (see e.g. Ref. ), chiral quark models (see e.g. Ref. ), heavy-quark effective theory , the linear $`\sigma `$ model (L$`\sigma `$M) , Nambu–Jona-Lasinio models , (unitarized) chiral perturbation theory ), meson-exchange models (see e.g. Ref. ), and so on. While it lies outside the scope of this short paper to describe in detail the respective merits and drawbacks of all these approaches, on can generally state that a good reproduction of either spectra or scattering data is achieved, often at the cost of a non-negligible number of free fit parameters, but without accomplishing a unified description of mesonic spectra and meson-meson scattering. In our view, however, spectra and scattering of mesons are inexorably intertwined concepts, considering that most mesons are resonances, decaying strongly into final states of two or more lighter mesons.

Therefore, we believe one should treat the mechanisms responsible for quark confinement and for mesonic decay on an equal footing. This we have realized already in an essentially non-relativistic (NR) framework, employing the coupled-channel Schrödinger formalism, with several confined $`q\overline{q}`$ channels and many free two-meson channels. For a detailed discussion of the configuration-space model, we refer to a separate contribution  to this workshop, so let us just emphasize here some important non-perturbative features. From the spectroscopic point of view, meson-loop effects on the $`q\overline{q}`$ spectra are included to all (ladder) orders, even those resulting from closed meson-meson channels. Conversely, the strong influence of $`s`$-channel $`q\overline{q}`$ states on non-exotic meson-meson scattering is taken into account. Furthermore, bound states, resonances, and scattering observables all follow directly from an analytic $`S`$-matrix derived in closed form, thus dispensing with any kind of perturbative methods, which could be highly untrustworthy for the often very large couplings involved. As a matter of fact, perhaps the most remarkable conquest of the model is its parameter-free prediction of the light scalar mesons , including the $`\kappa `$(800) ($`K_0^{}`$) , which has very recently been experimentally confirmed . Precisely for the light scalars, perturbation theory manifestly breaks down .

Other coupled-channel, also called unitarized, approaches are due to Eichten et al. (Ref. , first paper), Bicudo & Ribeiro (Ref. , fourth paper), and Törnqvist with collaborators (Ref.  and references therein). The conclusions of the latter two models are in qualitative agreement with ours, in the sense that large shifts are found owing to unitarization, even for (bare) states below the lowest decay threshold. However, in the case of the scalar mesons, important differences between our model and especially the one of Ref.  come to light, as the $`\kappa `$(800) is not found in the latter reference (see Ref.  for a detailed comparison).

Notwithstanding its successes, our coordinate-space model has, of course, a number of shortcomings. First of all, the light pseudoscalar mesons are quite badly off, except for the kaon. This may be due to the largely NR formalism without Dirac spin structure , the disregard of chiral symmetry, and possibly also the neglect of scalar-pseudoscalar two-meson channels. Moreover, for most processes, the experimental meson-meson phase shifts are only roughly and not very accurately reproduced, which is no wonder in view of the neglect of final-state interactions from direct $`t`$-channel meson exchange, and the absence of any “fitting freedom”. Then there are problems with pseudothresholds in heavy-light systems, owing to the non-covariance of minimally relativistic equations in configuration space. Finally, too many, viz. 3, parameters — albeit a very modest number as compared to other approaches — are needed for the $`{}_{}{}^{3}P_{0}^{}`$ transition potential, due to the lack of a microscopic description in terms of quark exchange.

Nevertheless, we are convinced one can deal with these shortcomings by reformulating the model in momentum space, in a covariant three-dimensional (3D) quasipotential framework. Here, one could object why we do not attempt to tackle the full Bethe-Salpeter (BS) equation. The reason is twofold: first of all, the numerical effort to solve the resulting coupled 4D integro-differential equations with many channels would be enormous; but perhaps more importantly, the confinement mechanism employing potentials is an inherently 3D concept, possibly leading to a pathological, non-confining behavior when generalized to four dimensions (see e.g. Ref. , third paper). Anyhow, a covariant formulation would also allow to include, at a later stage, the Dirac spin structure of the quarks, and address the issue of dynamical chiral-symmetry breaking and generation of constituent mass. The resulting multichannel equations in momentum space we shall show to be explicitly solvable, provided we stick to harmonic confinement and work on a spline basis.

This paper is organized as follows. In Sec. 2 the good old harmonic oscillator (HO) is solved in momentum space using the spline technique, for two different 3D relativistic equations. In Sec. 3 a two-channel toy model, with one harmonically confined $`q\overline{q}`$ channel coupled via quark exchange to one free two-meson channel, is formally solved in closed form, employing again a spline expansion. Conclusions and an outlook on how to further develop the model are presented in Sec. 4.

## 2 RELATIVISTIC HARMONIC OSCILLATOR

In this section, we review the well-known 3D HO, first in coordinate space, and then in momentum space. In the latter representation, a 3D relativistic generalization is straightforward.

### 2.1 Non-Relativistic Harmonic Oscillator

The equal-mass radial Schrödinger equation (SE) for a local potential reads

$$\left\{\frac{1}{m}(\frac{d^2}{dr^2}+\frac{l(l+1)}{r^2})+V(r)\right\}u(r)=Eu(r).$$
(1)
If we take the HO potential

$$V(r)=\frac{1}{2}\mu \omega ^2r^2,\mu =\frac{1}{2}m,$$
(2)
define $`x\sqrt{\mu \omega }r`$, and assume for the moment $`l=0`$, we get

$$\left\{\frac{d^2}{dx^2}+x^2\right\}u(x)=\frac{E}{\frac{1}{2}\omega },u(x).$$
(3)
The spectrum of the 3D HO potential is very well known, viz. ($`l=0`$)

$$E=\omega (2n+\frac{3}{2}),n=0,1,2,\mathrm{},$$
(4)
leading to the eigenvalues $`3,7,11,\mathrm{}`$ for Eq. (3).

Alternatively, we can work in momentum space. The homogeneous Lippmann-Schwinger (LS) equation for the bound-state vertex function reads

$$\mathrm{\Gamma }(\stackrel{}{p})=\frac{d^3p^{}}{(2\pi )^3}V(\stackrel{}{p},\stackrel{}{p}^{})G_0(\stackrel{}{p}^{};E)\mathrm{\Gamma }(\stackrel{}{p}^{}),$$
(5)
where the free two-body Green’s function is given by

$$G_0^1(\stackrel{}{p}^{})=E\frac{\stackrel{}{p}^2}{m}+iϵ.$$
(6)
Now define the wave function $`\mathrm{\Psi }G_0\mathrm{\Gamma }`$. Then Eq. (5) becomes

$$G_0^1(\stackrel{}{p})\mathrm{\Psi }(\stackrel{}{p})=\frac{d^3p^{}}{(2\pi )^3}V(\stackrel{}{p},\stackrel{}{p}^{})\mathrm{\Psi }(\stackrel{}{p}^{}),$$
(7)
or equivalently

$$\frac{\stackrel{}{p}^2}{m}\mathrm{\Psi }(\stackrel{}{p})+\frac{d^3p^{}}{(2\pi )^3}V(\stackrel{}{p},\stackrel{}{p}^{})\mathrm{\Psi }(\stackrel{}{p}^{})=E\mathrm{\Psi }(\stackrel{}{p}).$$
(8)
This is just the SE in momentum space. In this representation, the HO potential is given by

$$V(\stackrel{}{p},\stackrel{}{p}^{})=\frac{1}{2}\mu \omega ^2\delta ^{(3)}(\stackrel{}{p}\stackrel{}{p}^{})\mathrm{\Delta }_\stackrel{}{p},$$
(9)
which is local as in coordinate space. Substituting into (8) yields

$$\left\{\frac{\stackrel{}{p}^2}{m}\frac{1}{2}\mu \omega ^2\mathrm{\Delta }_\stackrel{}{p}\right\}\mathrm{\Psi }(\stackrel{}{p})=E\mathrm{\Psi }(\stackrel{}{p}).$$
(10)
Taking again $`l=0`$ and introducing $`v(\text{p})=\text{p}\mathrm{\Psi }(\text{p})`$, with $`\text{p}|\stackrel{}{p}|`$, we obtain

$$\left\{\frac{\text{p}^2}{m}\frac{1}{2}\mu \omega ^2\frac{d^2}{d\text{p}^2}\right\}v(\text{p})=Ev(\text{p}).$$
(11)
If we define the dimensionless variable $`y\text{p}/\sqrt{\mu \omega }`$ and write $`m=2\mu `$, we get

$$\left\{y^2\frac{d^2}{dy^2}\right\}v(y)=\frac{E}{\frac{1}{2}\omega }v(y).$$
(12)
This is identical to Eq. (3). Note that now the roles of the kinetic and the potential term are interchanged.

### 2.2 Blankenbecler–Sugar–Logunov–Tavkhelidze (BSLT) Equation

The momentum-space equal-mass (scalar) BSLT  equation for the wave function has exactly the same form as in the LS case of Eq. (7), but now with center-of-mass (CM) Green’s function

$$G_0^1(\stackrel{}{p})=\frac{E_p}{m^2}(\frac{s}{4}E_p^2+iϵ),$$
(13)
where $`E_p=\sqrt{\stackrel{}{p}^2+m^2}`$, and $`s`$ is the total CM energy squared $`s=P^2=(p_1+p_2)^2=(2m+E)^2`$. Note that $`G_0^{\text{BSLT}}`$ can be derived in a completely covariant fashion, and written down explicitly for an arbitrary frame . In the NR limit we have $`E_pm`$ and $`s4m^2+4mE`$, recovering the LS Green’s function (6). With the new $`G_0^1`$, we get instead of Eq. (8)

$$\frac{\stackrel{}{p}^2}{m}\mathrm{\Psi }(\stackrel{}{p})+\frac{m}{E_p}\frac{d^3p^{}}{(2\pi )^3}V(\stackrel{}{p},\stackrel{}{p}^{})\mathrm{\Psi }(\stackrel{}{p}^{})=\frac{s4m^2}{4m}\mathrm{\Psi }(\stackrel{}{p}).$$
(14)
In terms of the above-defined dimensionless variable $`y`$, the BSLT equation becomes

$$\left\{y^2\frac{1}{\sqrt{1+{\displaystyle \frac{\omega }{2m}}y^2}}\frac{d^2}{dy^2}\right\}v(y)=\frac{s4m^2}{2m\omega }v(y).$$
(15)
From this equation we readily see that deviations from the LS case are determined by the parameter $`\omega /2m`$. Furthermore, in the limit $`\omega /2m0`$, the eigenvalues $`3,7,11,\mathrm{}`$ get recovered. However, there are additional relativistic effects due to $`s4m^2+4mE`$.

### 2.3 Equal-Time (ET) Equation

The ET, also called Mandelzweig–Wallace , Green’s function can be obtained from the BSLT (or Salpeter) equation by adding an extra term $`G_0^C`$ to the two-body Green’s function $`G_0`$ which approximates crossed-ladder contributions. This extra term becomes exact in the eikonal limit, or when one of the particles gets infinitely heavy. Therefore, the covariant ET equation manifestly has the correct one-body limit, i.e., one recovers either the Klein–Gordon or the Dirac equation, for scalar or spin-1/2 particles, respectively. In the scalar case, the CM ET Green’s function is given by

$$G_0^1(\stackrel{}{p})=\frac{E_p^3}{m^2}\frac{{\displaystyle \frac{s}{4}}E_p^2+iϵ}{2E_p^2{\displaystyle \frac{s}{4}}}.$$
(16)
The final form of the $`S`$-wave ET equation for the HO potential in momentum space reads

$$\left\{y^2\frac{1}{\sqrt{1+{\displaystyle \frac{\omega }{2m}}y^2}}\left[2\frac{1}{1+{\displaystyle \frac{\omega }{2m}}y^2}\right]\frac{d^2}{dy^2}\right\}v(y)=\frac{s4m^2}{2m\omega }\left\{1\frac{{\displaystyle \frac{\omega }{2m}}}{\left(1+{\displaystyle \frac{\omega }{2m}}\right)^{\frac{3}{2}}}\frac{d^2}{dy^2}\right\}v(y).$$
(17)
### 2.4 Non-Zero Orbital Angular Momentum

For $`l0`$, the dimensionless LS equation becomes

$$\left\{\frac{d^2}{dx^2}+\frac{l(l+1)}{x^2}+x^2\right\}u(x)=\frac{E}{\frac{1}{2}\omega }u(x).$$
(18)
In order to numerically solve the LS and corresponding BSLT, ET equations, we substitute $`u(x)x^l\stackrel{~}{u}(x)`$. Then

$$\left\{\frac{d^2}{dx^2}+\frac{2l}{x^2}(1x\frac{d}{dx})+x^2\right\}\stackrel{~}{u}(x)=\frac{E}{\frac{1}{2}\omega }\stackrel{~}{u}(x).$$
(19)
The function $`\stackrel{~}{u}(x)`$ has the right boundary conditions to allow an easy numerical solution using an expansion on a spline basis (for details, see Ref.  and second paper of Ref. ). Thus, the NR HO spectrum is recovered with great accuracy, viz.

$$E=\omega (2n+l+\frac{3}{2}),n,l=0,1,2,\mathrm{}.$$
(20)
In the case of the BSLT and ET equations, the same can be done by replacing everywhere

$$\frac{d^2}{dy^2}\frac{2l}{y^2}(1y\frac{d}{dy})\frac{d^2}{dy^2}.$$
(21)
The LS, BSLT, and ET results for $`l=0`$ are depicted in Fig. 1. We see that the relativistic corrections are relatively modest for the BSLT equation, but can become huge in the ET case, especially for excited states. The latter large effect has to do with the sign change of the ET propagator (16) as a function of the relative momentum, and has also been observed in the Dirac case (see Ref. , second paper). However, this overshooting of the ET equation will strongly depend on the chosen Dirac structure of the confining force in a realistic model, and hence should not cause any worries at this stage.

## 3 UNITARIZED MODEL IN MOMENTUM SPACE

Consider now a simple two-channel model of one confined $`q\overline{q}`$ system coupled to one free meson-meson channel. In a unitarized picture, this can formally be described by a $`2\times 2`$ $`T`$ matrix, with the dynamical equation (LS, BSLT, or ET)

$$\left(\begin{array}{cc}T_{00}& T_{01}\\ T_{10}& T_{11}\end{array}\right)=\left(\begin{array}{cc}V_{00}& V_{01}\\ V_{10}& V_{11}\end{array}\right)+\left(\begin{array}{cc}V_{00}& V_{01}\\ V_{10}& V_{11}\end{array}\right)\left(\begin{array}{cc}G_0& 0\\ 0& G_1\end{array}\right)\left(\begin{array}{cc}T_{00}& T_{01}\\ T_{10}& T_{11}\end{array}\right).$$
(22)
Equation (22) is diagrammatically represented in Fig. 2.

Dashed lines are mesons, solid lines with arrows are (anti)quarks, and curly exchanges symbolize confinement. Furthermore, open circles represent meson-quark-antiquark vertices, and dots stand for three-meson vertices. Note that the communication between the $`q\overline{q}`$ and two-meson channels takes place through quark exchange. This will give rise to a generally non-local $`{}_{}{}^{3}P_{0}^{}`$ transition potential in configuration space, as soon as the quarks and mesons have different masses. The precise form of such a potential will also depend on possible form factors to be attributed to the various vertices.

In order to solve Eq. (22), we must realize that no asymptotic scattering of quarks is possible, so $`T_{00}`$, $`T_{01}`$, and $`T_{10}`$ must in the end be eliminated from the equations: only $`T_{11}`$ can be on shell. If we now take again harmonic confinement, Eq. (22) becomes a set of coupled integro-differential equations. Omitting for clarity angular-momentum indices, we get e.g. for $`T_{01}`$

$$T_{01}(\text{p},\text{p}^{})=V_{01}(\text{p},\text{p}^{})\frac{1}{2}\mu \omega ^2\mathrm{\Delta }_\text{p}G_0(\text{p})T_{01}(\text{p},\text{p}^{})+\frac{\text{k}^2d\text{k}}{4\pi ^2}V_{01}(\text{p},\text{k})G_1(\text{k})T_{11}(\text{k},\text{p}^{}).$$
(23)
Note that $`\text{p}^{}`$ is an essentially dummy variable in the equations, so we can choose a fixed value for it, say $`\overline{\text{p}}`$. If we now expand $`T_{01}`$ and $`T_{11}`$ on a spline basis and discretize the initial momentum p, then $`T_{01}`$ can be expressed in terms of $`T_{11}`$ as follows.

$$T_{01}(\text{p}_i,\overline{\text{p}})=\underset{j=1}{\overset{N}{}}t_{01}^js_{01}^j(\text{p}_i),v_{01}^iV_{01}(\text{p}_i,\overline{\text{p}});$$
(24)
$$T_{11}(\text{k},\overline{\text{p}})=\underset{j=1}{\overset{N}{}}t_{11}^js_{11}^j(\text{k});$$
(25)
$$A_{01}^{ij}\left\{1+\frac{1}{2}\mu \omega ^2\mathrm{\Delta }_{\text{p}_i}G_0(\text{p}_i)\right\}s_{01}^j(\text{p}_i);$$
(26)
$$A_{11}^{ij}\frac{\text{k}^2d\text{k}}{4\pi ^2}V_{01}(\text{p}_i,\text{k})G_1(\text{k})s_{11}^j(\text{k}).$$
(27)
Combining Eqs. (24)–(27), we can write

$$𝐀_{01}𝐭_{01}=𝐯_{01}+𝐀_{11}𝐭_{11}𝐭_{01}=𝐀_{01}^1\{𝐯_{01}+𝐀_{11}𝐭_{11}\},$$
(28)
where the boldface capitals and small letters are obvious shorthands for matrices and vectors, respectively. Turning next to $`T_{11}`$, we get the integral equation

$$T_{11}(\text{p},\overline{\text{p}})=V_{11}(\text{p},\overline{\text{p}})+\frac{\text{k}^2d\text{k}}{4\pi ^2}V_{11}(\text{p},\text{k})G_1(\text{k})T_{11}(\text{k},\overline{\text{p}})+\frac{\text{k}^2d\text{k}}{4\pi ^2}V_{10}(\text{p},\text{k})G_0(\text{k})T_{01}(\text{k},\overline{\text{p}}).$$
(29)
Similarly to Eqs. (24)–(27), we define now

$$B_{11}^{ij}s_{11}^j(\text{p}_i)\frac{\text{k}^2d\text{k}}{4\pi ^2}V_{11}(\text{p}_i,\text{k})G_1(\text{k})s_{11}^j(\text{k});$$
(30)
$$B_{01}^{ij}\frac{\text{k}^2d\text{k}}{4\pi ^2}V_{10}(\text{p}_i,\text{k})G_0(\text{k})s_{01}^j(\text{k}),v_{11}^iV_{11}(\text{p}_i,\overline{\text{p}}).$$
(31)
Then we can express $`𝐭_{11}`$ in terms of itself, using Eq. (28):

$$𝐁_{11}𝐭_{11}=𝐯_{11}+𝐁_{01}𝐭_{01}\stackrel{=}{\left(28\right)}𝐯_{11}+𝐁_{01}𝐀_{01}^1\{𝐯_{01}+𝐀_{11}𝐭_{11}\}.$$
(32)
This finally leads to

$$\{𝐁_{11}𝐁_{01}𝐀_{01}^1𝐀_{11}\}𝐭_{11}=𝐯_{11}+𝐁_{01}𝐀_{01}^1𝐯_{01},$$
(33)
$$𝐭_{11}=\{𝐁_{11}𝐁_{01}𝐀_{01}^1𝐀_{11}\}^1\{𝐯_{11}+𝐁_{01}𝐀_{01}^1𝐯_{01}\}.$$
(34)
Thus, we have obtained a closed-form expression for the $`T`$ matrix (hence the $`S`$ matrix), containing the complete information of the confinement-plus-scattering process. Needless to say that all these algebraic manipulations only make sense if the involved integrals actually exist and the spline expansions converge. Yet, we are confident this will be the case when appropriate boundary conditions are chosen, since the same spline expansion we are using here has been successfully applied to few-body scattering calculations in momentum space . Nevertheless, in the Dirac case, to be dealt with in the future, form factors will have to be included for the meson-quark-antiquark vertices so as to make the integrals involving three quark propagators convergent. In any case, such form factors are physically meaningful, and should therefore already be considered in the approximation with spinless fermions, whenever doing phenomenology.

Bound states and resonances are given by zeroes in the determinant of the inverted-matrix factor right after the equal sign in Eq. (34). Resonance poles in the second Riemann sheet can be searched for by analytic continuation to complex energies, using e.g. contour-rotation methods. Moreover, from the structure of the equations one readily sees that the inclusion of final-state interactions, due to $`t`$-channel meson exchange in the meson-meson channel, simply amounts to taking a non-vanishing $`V_{11}`$ in Eq. (29), which will not lead to any additional numerical effort. Finally, the generalization to many scattering channels and also more $`q\overline{q}`$ channels is straightforward, at the expense of bigger matrices.

## 4 CONCLUSIONS AND OUTLOOK

In the foregoing, we have demonstrated that the formulation of a unitarized quark/meson model with harmonic confinement is relatively easy in momentum space. Thus, the $`{}_{}{}^{3}P_{0}^{}`$ mechanism is dynamically described through $`P`$-wave quark exchange. Furthermore, the inclusion of final-state interactions via $`t`$-channel meson exchange is straightforward, and does not complicate the structure of the equations. The numerical resolution of these equations using splines looks very convenient, which should result in (approximate) analytic solutions, thereby facilitating an analytic continuation to complex energies in order to search for resonance poles.

Our first priority is, of course, to achieve a total control of the numerics, for real as well as complex energy. Then a simple toy model will be worked out so as to reproduce known coordinate-space results. Furthermore, cases problematic in the configuation-space formulation will be studied, such as heavy-light mesons. Also, the influence of $`t`$-channel meson exchange on meson-meson phase shifts will be investigated. In a somewhat more distant perspective, we intend to do a lot of phenomenology, and use the feedback to further refine the model.

Acknowledgments The authors thank F. Kleefeld and J. A. Tjon for valuable discussions, and P. Nogueira for expert graphical assistance. This work was partly supported by the Fundação para a Ciência e a Tecnologia (FCT) of the Ministério do Ensino Superior, Ciência e Tecnologia of Portugal, under contract number CERN/P/FIS/43697/2001.


format: text

text:
# Using nonlocal electrostatics for solvation free energy computations: ions and small molecules

## 1 Introduction

The unique properties of water, namely its polar nature, high dielectric constant, and its ability to form hydrogen bonds ball , are responsible for the existence of life as we know it. At the same time, these very properties are the main obstacle in modelling hydration effects. An accurate model of ion and small molecule hydration will therefore have important applications in computational chemistry, chemical engineering, and drug design.

In the latter, the interest stems from the fact that drug-like small molecules are usually moderately polar to achieve the required bioavailability. Therefore their solvation free energy is generally dominated by the electrostatic contribution, while the nonpolar contribution is more or less negligible. One of the most successful models for predicting solvation free energies of small molecules (and proteins) therefore employs the Poisson or Poisson-Boltzmann equations hon93 ; hon95 . Within this approach, the solvent is typically treated as a structureless medium of a given dielectric constant $`\epsilon `$ (for water, $`\epsilon 80`$.) The solute itself is treated as a low-$`\epsilon `$ “cavity” jac95 .

Clearly, this approach fails to take into account many relevant details of the solvent structure. The presence of long-ranged electrostatic fields induces static correlations among the polar molecules that may vary considerably in space depending on the molecular build-up of the solute. These can be taken into account by performing Molecular Dynamics simulations, thus leaving the continuum description at the expense of high computational cost even for small molecules levy .

Including more structural detail into the continuum description on the one hand without increasing computational complexity on the other is therefore highly desirable. At the same time, however, also the avoidance of additional parameter sets, which is common practice in the use of continuum electrostatics in computational chemistry, is a key requirement.

In this work, we follow one particular route to achieve this goal. We develop an approach for the computation of the electrostatic (polar) contribution to the free energy of solvation based on nonlocal electrostatics nonlocal ; kornyshev1 ; kornyshev2 ; sutmann ; basilevsky1 ; basilevsky2 . This continuum approach, originating in the physics literature as a generalization of classical electrostatics to account for media with spatial dispersion, is developed here within a physically transparent minimal parameter set and applied to the simplest systems of single ions and small molecules for which experimental and theoretical results are available for quantitative comparison. In particular we address the issue of a systematic choice of ion radii which most commonly is based on essentially empirical parametrizations lee . We conclude with an outlook on extensions of our approach for more complex, and biologically relevant, molecules.

## 2 Nonlocal electrostatics of the solvent

The electrostatic potential $`\varphi `$ of a charged molecule in solution is, within classical (local) electrostatics, given by the Poisson equation jackson

$$\left[\epsilon (𝐫)\epsilon _0\varphi (𝐫)\right]=\varrho (𝐫)$$
(1)
where $`\epsilon _0`$ denotes the dielectric constant of vacuum, and $`\epsilon (𝐫)`$ is a local dielectric function which is taken as the constants $`\epsilon _{water}78`$ of the solution (water), and a much smaller value within a solute molecule (for proteins, e.g., a typical value taken is $`\epsilon 25`$). The rhs of eq.(1) is the charge density giving rise to the potential.

This arguably simple description is prone to complications: first, the transition region between solute and solution is ill-described and needs a further (rather delicate) modeling. Second, very little is said about structural effects due to the orientation of the polar water molecules near highly charged regions of a given solute molecule.

The simplest way to introduce structural effects into a continuum description of the solvent is to account for correlations due to polarization effects between solution molecules characterized by a correlation length $`\lambda `$. This gives rise to a nonlocal generalization of eq.(1) given by nonlocal

$$_𝐫𝑑V^{}\epsilon (|𝐫𝐫^{}|)\epsilon _0_𝐫^{}\varphi (𝐫^{})=\varrho (𝐫)$$
(2)
where isotropy of the liquid medium is supposed. Eq.(2) needs to be accompanied by suitable boundary conditions (see basilevsky1 ). These conditions can be simplified considerably by assuming the atomic constituents as conducting spheres, which is exact for the case of ions, and approximate for the case of small molecules. The solution of the nonlocal solvation problem with a treatment of the full boundary conditions is beyond the scope of this work and will be reported elsewhere we .

The rationale behind eq.(2) is the following. Since the electrostatic fields are long-ranged, the solvent molecules will feel the presence of the fields of other molecules over characteristic distances (that will be of particular relevance below). Maintaining locality of the physical fields (i.e., electrostatic potentials and electric fields), it is the response functions into which this purely static correlation effect needs to be embedded. We note that the dielectric function is related to the dielectric response function via $`\epsilon =1+\chi `$, irrespective of the local or nonlocal character of the theoretical description.

Due to the assumption of spatial isotropy, eq.(2) lends itself to a treatment in Fourier space, reducing the linear integro-differential equation eq.(2) to an algebraic equation involving the Fourier-transformed function $`\epsilon (|𝐤|)`$ with wavevector $`𝐤`$. Explicit functional dependences of either $`\epsilon (|𝐫𝐫^{}|)`$ or $`\epsilon (|𝐤|)`$ have been derived from various approaches, e.g. within a Ginzburg-Landau theory for the polarization fields sutmann .

Motivated by these works we have formulated a family of nonlocal dielectric functions fulfilling the following minimal requirements anhi-diplom  :

i) for high $`k`$-values (equivalently, on small spatial scales) the solvent molecules cannot follow the polarization forces.

Thus, a limiting value for the macroscopic dielectric function has to be reached in this limit. Typical choices in the literature are $`lim_k\mathrm{}\epsilon (k)\epsilon _{\mathrm{}}1.8`$ sutmann ; hasted , or $`\epsilon _{\mathrm{}}=1`$ landau . We have considered both cases and do not find a significant dependence on this choice;

ii) in the opposite limit, $`k0`$, the value for $`\epsilon `$ should equal the macroscopic value;

iii) the physical length scale of the polarization fluctuations is characterized by a correlation length $`\lambda `$. The local electrostatic limit is given by $`lim_{\lambda 0}\epsilon (|𝐫𝐫^{}|)=\epsilon _{loc}\delta (𝐫𝐫^{})`$;

iv) causality conditions (Kramers-Kronig relations) need to be fulfilled landau .

While the conditions i)-iv) are clearly not sufficient to determine the dielectric function unambiguously, we find that the following class of functions in Fourier space fulfill these criteria sufficiently well fourier

$$\epsilon (k)=\frac{1}{(2\pi )^{3/2}}\left[\epsilon _{\mathrm{}}\frac{\epsilon _{\mathrm{}}\epsilon _{loc}}{(1+\lambda ^2k^2)^n}\right],$$
(3)
which depend on only two parameters, $`n`$, and $`\lambda `$. For $`n=2`$, the dielectric function decays exponentially in real-space, while $`\epsilon (k)`$ is a Lorentzian for $`n=1`$, decaying in real-space as $`\mathrm{exp}(r/\lambda )/r`$. Note that formally the case $`n=1`$ was discussed previously in kornyshev1 ; kornyshev2 , based, however, on a different physical interpretation (see also ritschel  for a similar approach in the context of nuclear physics). We have also tested alternative choices, e.g. a Gaussian model and a nonlocal model derived from a Ginzburg-Landau theory for polarization fields sutmann .

The comparison of our nonlocal models with other approaches and experimental data is made possible by using the concept of an effective local dielectric function $`\widehat{\epsilon }(r)`$. It is usually defined for the potential of a point charge $`q`$,

$$\varphi (r)=\frac{1}{4\pi \epsilon _0}\frac{q}{\widehat{\epsilon }(r)r}$$
(4)
where

$$\widehat{\epsilon }(r)\left[\frac{d𝐤}{(2\pi )^{5/2}}\frac{\mathrm{sin}(kr)}{kr}\frac{1}{\epsilon (k)}\right]^1.$$
(5)
For a general charge distribution $`\widehat{\epsilon }(r)`$ can be defined as well we . We note that our determination of the effective dielectric constant is consistently performed within the nonlocal continuum theory and does not rely on data-based fits meh84 ; mallik .

Expression (5) allows to test our approach against previous results. We took empirical models for the radial dependence of $`\epsilon `$ derived for experimental results from Mehler and Eichele (ME) meh84  and Conway (CO) con51 . The CO-model is reproduced best by the Lorentzian ($`n=1`$) for $`\lambda =15`$Å, while the ME-model is reproduced best for $`\lambda =24.13`$Å. The purely exponential model ($`n=2`$) displays unphysical singularities in $`\widehat{\epsilon }(r)`$ at small $`r`$, while with a choice of $`\lambda =5`$Å, the CO-model is reproduced for larger distances, $`r>15`$Å. Figure 1 shows $`\widehat{\epsilon }(r)`$, as computed by eq.(5) for a single sodium ion, using the Lorentzian model with $`\lambda =15`$Å. The computation of $`\widehat{\epsilon }(r)`$ performed here makes use of a specifically adapted FFT applicable also to molecules we .

A Gaussian choice for $`\epsilon (k)`$, however, leads to strong oscillations in the effective dielectric function, although its overall shape resembles the CO-model for $`\lambda =5`$Å. Checking our approach against the nonlocal model by Sutmann et al. sutmann , we find that their theory fails to comply with our requirement iii), i.e., it does not reproduce the correct limiting value at large distances.

## 3 Single ions and the choice of their radii

We now turn to the application of our approach to small ions. We treat these ions as Born spheres. A standard problem in the definition of Born-type ions is the definition of the ion radius lee . Starting from the accepted interpretation that a solvated ion is surrounded by solvation shells, the first of these shells will be “as close as possible” to the ion. We can therefore identify the position of the centers of the oxygen in the first solvation shell with the first peak in the radial distribution function (rdf) which can be obtained either from scattering experiments or from molecular dynamics simulations.

We defined the radius of the water oxygen as half the position of the first peak in the oxygen-oxygen radial distribution function (rdf) of bulk water. The ion radii were then derived by subtracting this radius from the position of the first peak of the ion-oxygen rdf derived from a molecular dynamics simulation.

As a specific input for our calculations we take the values for the ion radii determined by Åqvist Aqv90 , obtained from a combined free energy perturbation/force field approach based on the GROMOS force field vGb , employing several different water models (flexible SPC, rigid SPC, and TIP3P). Based on this input, we have determined the solvation free energy from nonlocal electrostatics, which is given by

$$\mathrm{\Delta }G^{polar}=\frac{1}{2\epsilon _0}𝑑𝐤\left\{1\frac{1}{(2\pi )^{3/2}\epsilon (𝐤)}\right\}D^2(𝐤)$$
(6)
where $`𝐃(𝐤)`$ is the dielectric flux density.

The result of the calculation is shown in Figs. 2 for monovalent ions and in Fig. 3 for divalent ions, for all models to be compared here. Note that we chose the parameter $`\lambda `$ based on its fit to the effective dielectric function. The CO-model was reproduced by $`\lambda =5`$Å for $`n=2`$ (exponential model), by $`\lambda =15`$Å for $`n=1`$ (Lorentzian) and by $`\lambda =5`$Å for the Gaussian model. The ME-model could only be reproduced by $`n=1`$ with $`\lambda =24.13`$Å. The value for the exponential model is in accord with the findings in Ref. basilevsky1 . Also note, as shown in Table 1, that the values we obtain are only marginally corrected (if at all) by effects due to nonlinear saturation, which contributes in principle as well to a reduction of the dielectric constant near the ion.

The model by Sutmann was used with two different limiting values (Sutmann 1 with $`\epsilon _{\mathrm{}}=1`$, and Sutmann 2 with $`\epsilon _{\mathrm{}}=1.8`$). As the figures show, our results are consistently better than all the other theoretical curves.

The computation of electrostatic solvation free energies is also possible by employing eq.(5) in a standard (local) Poisson solver. We have implemented this in an available library (BALL) okhpl ; we . Fig. 4 compares the results of our computations for the monovalent ions with the nonlocal theory based on eq.(6) and from the effective local dielectric function, eq.(5), demonstrating the consistency of both approaches. Fig. 5 finally compares the theoretical results of the local and nonlocal approaches, obtained with the Poisson solver, with the experiment values. Evidently, the nonlocal approach yields results consistently superior to those obtained from the standard local theory.

We stress the significant advantage of our computations to work with first-principle radii without arbitrary adjustments. In Figs. 2-5, the Åqvist-radii were used Aqv90 . In addition, we have also tested the Shannon-radii derived from X-ray crystal data sha1 ; sha2 , without significant effect on our results. Our result thus give a basis to the general belief that the need to introduce effective radii is in fact a consequence of the local water structure around the ions sandberg . We believe that our nonlocal approach demonstrates that it is therefore preferable conceptually to introduce the length-scale governing the structural effects in the solvent rather than introducing necessarily artificial procedures to adjust the Born radii.

## 4 Small molecules

As a second application we apply our approach to determine the solvation free energy for small alcohols for which we can assume that the polar contribution exceeds the nonpolar one. Fortunately, accurate measurements of the solvation free energies of these molecules are available. We write the charge distribution of these molecules as a linear superposition of radially symmetric partial distributions translated by a vector $`𝐑`$, $`\varrho _i=\varrho (|𝐫+𝐑_i|)`$, i.e. $`\varrho (𝐫)=_{i=1}^N\varrho _i(|𝐫+𝐑_i|)`$. Its Fourier transform is then given by the expression $`\varrho (𝐤)=_{i=1}^N\varrho _i(k)\mathrm{exp}(i𝐑_𝐢𝐤).`$

Again we have to define the radii of the atoms. For this we chose to classify the atoms into classes depending on their chemical environment in the molecule (e.g., like the hydrogen atom in an OH-group), expecting that the radii of all atoms in a certain class are more or less similar. For methanol ($`\mathrm{C}H_3OH`$), e.g., we defined four classes of atoms and used their rdf’s with the oxygen of water to find the following set of radii: Methyl C: 2.135 Å, hydroxyl O: 2.014 Å, hydroxyl H: 1.115 Å, methyl H: 1.394 Å.

Compared to the calculation of the free energies for the ions this computation is slightly more involved as it requires a three-dimensional integration instead of a one-dimensional one. Within the framework of the simplified boundary conditions, the unperturbed fields of the atomic constituents are simply superimposed. To ensure short computation times, in our implementation we have used the VEGAS-Monte-Carlo integration scheme VEGAS  that is supplied with the GNU scientific library GSL gsl . Details of these computations will again be given elsewhere we .

The results for some small alcohols are shown in Table 1. For the comparison we have used the exponential with $`\lambda =5`$ Å, the Lorentzian with $`\lambda =15`$ Å, and the Lorentzian with $`\lambda =24.13`$ Å. All values in these tables are given in kJ/mol.

The interpretation of the results is more complicated than in the Born ion case. While there we could assume that the contribution of the nonpolar part of the free energy of solvation could be neglected (for the Born ions it can be estimated to be of the order of $`1020`$ kJ/mol), this is not the case for the alcohols. The electrostatic contribution is still the dominant part, but not the only significant one. We therefore applied a very simple model for the nonpolar contribution uhlig  in order to be able to compare our results to experimental data. The results given for the polar contribution can also be optimized by considering improved charge distributions we .

## 5 Conclusions and outlook

We have demonstrated that the approach to the computation of solvation free energies based on nonlocal electrostatics is able to reproduce experimental data for ions and small molecules with reasonable accuracy. The approach we put forward has the basic advantage to rely on essentially one parameter which has a transparent physical interpretation as the correlation length of the polarization fluctuations. While within the present paper this parameter is used as a fitting parameter, it should be clear that it might be also determined by experiment or simulation. A further significant advantage is that our approach does not make use of the commonly used adjustments of atomic or ionic radii. A challenge for the future will be the extension of our approach to determine electrostatic solvation free energies for more complex and biologically relevant molecules. Work in this direction is in progress we .

Acknowledgement. This work is supported by the DFG under Schwerpunktsprogramm “Informatikmethoden zur Analyse und Interpretation grosser genomischer Datenmengen” (grant LE952/2-1).


text:
# A source of cold atoms for a continuously loaded magnetic guide

## I Principle of operation

In order to inject atoms into the magnetic guide, we use a magneto-optical trap (MOT) based upon four laser beams in a tetrahedral configuration superimposed with a magnetic two-dimensional quadrupole field. The setup consists of a moving molasses  in the longitudinal direction combined with cooling and confining forces in the transverse directions and will be called injecting MOT in the following. A detailed description can be found in . The injecting MOT allows to capture atoms entering the trapping volume with a velocity $`v<v^{}`$, where $`v^{}`$ is a function of the available laser power, and to subsequently slow the atoms down to the final velocity $`v_\mathrm{b}`$ which only depends on the relative frequencies of the laser beams. In this way we produce an atomic beam with a temperature on the order of 50 $`\mu `$K and a velocity in the range of 30 cm/s to 3 m/s.

If the injecting MOT is loaded from a rubidium background gas with vapor pressure $`p`$, it captures atoms at a rate proportional to $`p`$. However, collisions between captured and thermal background atoms will limit the flux $`\mathrm{\Phi }`$ of outgoing atoms to

$$\mathrm{\Phi }p\mathrm{exp}(\frac{\alpha p}{v_\mathrm{b}}),$$
(2)
where $`\alpha `$ depends on the characteristic length of the trap and the collision cross section. This process sets an upper limit $`\mathrm{\Phi }_{\mathrm{max}}=\mathrm{max}_{\{p\}}(\mathrm{\Phi })`$ to the flux that can be achieved. It decreases as $`\mathrm{\Phi }_{\mathrm{max}}v_\mathrm{b}`$ with decreasing velocity $`v_\mathrm{b}`$. Since $`\gamma \mathrm{\Phi }/v_\mathrm{b}`$, it also limits the collision rate in the magnetic guide provided that the temperature and the confining force of the guide are unchanged. In our previous setup , we were able to produce a flux of $`10^9`$ atoms/s at $`v_\mathrm{b}=2.5`$ m/s and a pressure $`p=4\times 10^8`$ mbar. In order to overcome this limitation, we have to operate the injecting MOT in an ultra-high vacuum chamber.

For this purpose, many techniques have been developed. Thermal beams can be slowed down by means of a Zeeman slower . Alternatively, a collimated beam of atoms can be obtained by a vapor-loaded MOT with a leak at its center. For instance, one can drill a hole in one of the mirrors of a MOT ), or use a pyramidal mirror structure with a hole at its vertex , or use an extra pushing beam, which destabilizes the MOT at its center . A third method for producing a cold beam of atoms relies on a 2D-MOT . All these sources produce a relatively intense beam (from 106 to several $`10^{10}`$ atoms per second), with an average velocity $`\overline{v}`$ between 10 and 50 m/s. In our setup, we have implemented the 2D-MOT scheme. We note that the relation (2) also applies to the source producing the atomic beam. However, it is much less restrictive since the atomic velocities at its exit are orders of magnitude higher than in the injecting MOT.

## II Experimental setup

A schematic drawing of the experimental setup is shown in fig. 1. The first MOT is located in a rectangular quartz cell (80 mm $`\times `$ 45 mm $`\times `$ 45 mm) connected to a Rubidium reservoir by a valve (not shown in fig. 1) and joined to the main part of the vacuum system by a differential pumping tube. It is a vapor-loaded two-dimensional MOT based on a design that has been characterized in detail in . For a thermal atom to be captured by the MOT, firstly it needs to have a transverse velocity smaller than the MOT’s capture velocity. Secondly, the atom’s longitudinal velocity $`v_{}`$ has to be low enough so that the transit time $`\tau =l/v_{}`$ ($`l`$ being the length of the MOT) is sufficiently long for the atom to be transversally trapped and cooled. In this way the 2D-MOT produces two beams of atoms which leave the trap on axis in opposite directions.

Two pairs of elongated coils in anti-Helmholtz configuration produce a two-dimensional quadrupole magnetic field whose zero-field line is parallel to the optical table and centered with respect to the glass cell. Typical field gradients are on the order of 12 G/cm. The cooling laser light for the MOT is produced by a tapered amplifier providing up to 450 mW of light. The amplifier is injected by a grating-stabilized diode laser which is red-detuned by $`\mathrm{\Delta }2\mathrm{\Gamma }`$ from the $`5S_{1/2},F=25P_{3/2},F=3`$ transition. A second grating-stabilized diode laser tuned to the $`5S_{1/2},F=15P_{3/2},F=2`$ transition is used to injection-lock another laser diode which provides up to 30 mW of repumping light. The cooling and repumping laser beams are superposed and split into two beams which are subsequently circularly polarized. After widening the beams by telescopes to a beam waist of 13 mm, the beams are steered to the quartz cell which they pass perpendicular to each other and to the zero-field line of the magnetic field. They are then retro-reflected after having passed another quarter-wave plate, thus providing the third and forth beam of the MOT. We find that the resulting intensity imbalance between the ingoing and the retro-reflected beams does not significantly perturb the functioning of the trap. Typically, the trap is operated at a total rubidium background pressure of $`3\times 10^7`$ mbar. It is located close to the entrance of a differential pumping tube, with a diameter $`d=7`$ mm and a length $`l=450`$ mm, capable of maintaining a pressure ratio of more than two orders of magnitude between the first quartz cell and the cell of the injecting MOT. The total distance between the first and the second MOT is $`L800`$ mm.

We use two methods to study the properties of the 2D-MOT. By measuring the fluorescence of the injecting MOT, we obtain a signal which is proportional to the flux of atoms that is produced by the 2D-MOT and recaptured by the injecting MOT. Even though this quantity is the most important parameter for the next stage of the experiment, it does not allow to characterize the performance of the 2D-MOT in a quantitative and independent way. Therefore, to gather more information we perform absorption measurements on the beam of atoms that arrives in the glass cell of the recapturing MOT. For this purpose, we shine a probe laser beam perpendicular to the atomic beam into the glass cell of the recapturing MOT. The beam is resonant with the $`5S_{1/2},F=25P_{3/2},F=3`$ transition, has a diameter of 1 mm and is circularly polarized. A weak magnetic field in the direction of the probe beam ensures that the atoms crossing the probe beam are optically pumped to the outermost Zeeman state. By measuring the absorption of the probe beam, we obtain the column density of the atomic beam as a function of the probe beam position (in the direction perpendicular to the atomic beam). Information about the velocity distribution of the atoms is obtained by performing time-of-flight measurements. After suddenly turning off the repumping laser of the 2D-MOT, we measure the atomic velocity dispersion by monitoring the absorption of the probe beam as a function of time. Note that this method works since the velocity distribution $`\mathrm{\Delta }v`$ is on the order the mean velocity $`\overline{v}`$ and since the distance between the probe beam and the 2D-MOT is large compared to the size of the MOT. In order to increase the signal-to-noise ratio of the measurement, we dither the frequency of the probe laser at 100 kHz and use a lock-in amplifier to measure the absorption.

## III Experimental results

For the 2D-MOT with circular beam profiles described in the previous section, we find that the number $`N`$ of atoms recaptured by the injecting MOT increases approximately linearly as a function of the cooling laser power $`P`$ when $`P`$ was varied between 16 mW and 160 mW per beam. When we reduce the cooling volume of the trap by introducing knife edges which limit the beam size, $`N`$ decreases dramatically in the case where the length of the trap is reduced while it is much less affected if the transverse beam size is reduced. It has already been observed and explained in  that elliptical beam profiles are advantageous for maximizing the atomic flux. Here, instead of using cylindrical lenses to create elliptical beams, we split each of the MOT beams into two beams before passing them through the telescopes. Each pair of beams is now passed through the glass cell side to side to each other so as to produce a 2D-MOT with elongated beam profiles (see fig.1). The laser power of the tapered amplifier can be arbitrarily repartitioned among the four beams, and we find that a maximum of recaptured atoms is observed for a configuration where the intensity of the two beams that are closer to the trap exit is lower than the intensity of the other two beams. The measurements presented in the following have been performed with this optimized configuration.

Absorption measurements of the atomic beam are done at a distance $`L800`$ mm from the 2D-MOT. The maximum absorption is on the order of $`0.7\%`$. Even though the absorption signal is small, it allows to calibrate those measurements that take advantage of lock-in detection. The column density profile of the atom beam in fig. 2 is recorded by scanning the probe beam vertically (i. e., in a direction perpendicular to the atomic beam). The full width at half maximum (FWHM) of about 10 mm is superior to the diameter $`d`$ of the differential pumping tube so that one might wonder whether the transverse beam size could be limited by the tube. If the transverse size of the 2D-MOT was negligible and the spatial density of the atomic flux was homogeneous at the tube exit, the linear density profile would have the shape of a half circle with diameter $`d_c=dL/l12.4`$ mm, giving a FWHM of 10.8 mm. A finite transverse trap size and the velocity-dependence of the beam profile due to gravity contribute to enlarge the beam size.

The velocity distribution of the atomic beam is obtained by rapidly switching off the repumping beam of the 2D-MOT and measuring the time-dependent absorption signal $`s(t)`$ as shown in fig. 3.

The atomic flux density per velocity class $`\rho _\mathrm{\Phi }(v)`$ is related to $`s(t)`$ by

$$\rho _\mathrm{\Phi }(v)t\frac{\mathrm{d}s}{\mathrm{d}t}\text{with}t(v)=L/v.$$
(3)
After calibrating the signal and integrating over the linear density profile, we obtain the atomic flux density which is shown in fig. 4. The beam is found to have a mean velocity $`\overline{v}=dvv\rho _\mathrm{\Phi }/dv\rho _\mathrm{\Phi }=38`$ m/s and a rather large velocity spread $`\mathrm{\Delta }v=(v\overline{v})^2^{1/2}=17`$ m/s. The highest flux that we have measured was $`\mathrm{\Phi }10^{10}`$ atoms/s. This number includes atoms with velocities up to 100 m/s. The experimentally most relevant quantity is the atomic flux that the injecting MOT is able to recapture. The capture velocity that can be achieved strongly depends on the trap geometry as well as the laser power available for the injecting MOT and will typically lie in a range between 25 m/s and 40 m/s. The fraction of atoms with a velocity below 35 m/s is 50%, which corresponds to a useful flux of $`5\times 10^9`$ atoms/s.

In this article, we have reported on the performance of a 2D-MOT similar to the one described in . In our setup, we use two adjacent 2D-MOT which allows to maximize the atomic flux by redistributing the available laser power among the two traps. This source has been developed for loading a moving molasses MOT which serves to inject a cold and slow atomic beam into a magnetic guide.

## acknowledgements

This work was supported by the Bureau National de la Métrologie, the Délégation Générale de l’Armement, the Centre National de la Recherche Scientifique and the Région Ile de France. C. F. Roos acknowledges support from the European Union (contract HPMFCT-2000-00478).


text:
# Nanoscale self-affine surface smoothing by ion bombardment

## Abstract

Topography of silicon surfaces irradiated by a 2 MeV Si+ ion beam at normal incidence and ion fluences in the range $`10^{15}10^{16}`$ ions/cm2 has been investigated using scanning tunneling microscopy. At length scales below $``$ 50 nm, surface smoothing is observed; the smoothing is more prominent at smaller length scales. The smoothed surface is self-affine with a scaling exponent $`\alpha =0.53\pm 0.02`$.

PACS no. 61.16.Ch; 61.80.Jh; 68.35.Bs; 68.35.Ct

One of the fundamental problems in materials science is to understand the effects of particle radiation on solid surfaces. The evolution of solid surface topography during ion-beam irradiation is governed by the interplay between the dynamics of surface roughening due to sputtering and smoothing due to material transport during surface diffusion. These competing processes are responsible for the creation of characteristic surface features like quasiperiodic ripples and self-affine topographies. These have been observed in the ion energy regime where sputtering is dominant and ion incidence is tilted to the surface normal. Although there is a large number of observations of ripple formation there are only a few studies on the scaling of the surfaces evolved in ion bombardment. A common feature of most rough surfaces observed experimentally or in discrete models is that their roughness follows simple scaling laws. Surface root-mean-square roughness $`\sigma `$ is defined as $`\sigma =<[h(x,y)h]^2>^{1/2}`$, where $`h(x,y)`$ is the surface height at a point $`(x,y)`$ on the surface and $`h`$ is the average height. The surface is termed self-affine if $`\sigma `$ changes with the horizontal sampling length $`L`$ according to $`\sigma \mathrm{}L^\alpha `$, where $`0<\alpha <1`$ is the roughness exponent. The roughness exponent quantifies how roughness changes with length scale and its value is indicative of the surface texture.

For graphite bombarded with 5 keV Ar ions at an angle $`\theta =60^{}`$ with respect to the surface normal, Eklund et al. reported $`\alpha 0.20.4`$, consistent with the predictions of the Kardar-Parisi-Zhang (KPZ) equation in 2+1 dimensions. Krim et al observed a self-affine surface roughness generated by 5 keV Ar ion bombardment of an Fe thin film sample at $`\theta =25^{}`$, with a scaling exponent $`\alpha `$=0.53, with no theoretical model predicting this value. In all these cases an increase of surface roughness was observed due to ion bombardment. Since ion arrival on the surface is a stochastic process and sputtering events are spatially distributed and of variable magnitude, surfaces are generally roughened during bombardment. In all the studies mentioned above the conditions are such that the erosion of the surface due to sputtering in ion bombardment is dominant over surface atomic diffusion. However, if the surface atomic diffusion dominates over sputtering, surface smoothing rather than roughening can occur. Carter and Vishnyakov have shown that inclusion of a directed flux of atoms parallel to the surface, generated by ion bombardment, in a stochastic differential equation description of the dynamics of surface evolution during sputter-erosion can induce smoothing for near-normal ($`\theta 0`$) ion incidence. The flux of atoms parallel to the surface provides an effective diffusion causing surface smoothing which competes with the roughening caused by sputtering. For $`\theta 0`$, roughening is weak as sputtering yield is small and smoothing dominates. Indeed for an ion incidence angle $`\theta 0`$, surface smoothing have been observed in ion bombardment over a large range of ion energies. Although some observations of surface smoothing have been reported, to our knowledge there has been no scaling studies of ion-beam induced surface smoothing. In scaling studies for nonequilibrium film growth by deposition, a value of $`\alpha 0.35`$ is expected when surface mobility of deposited particles are not allowed and $`\alpha `$=0.66 is expected when surface mobility is allowed. For ion-induced roughening the observed value of $`\alpha `$=0.2$``$0.4 is in reasonable agreement with the exponent for growth without surface diffusion. For ion-beam induced smoothing, where surface diffusivity is important, one may expect a different value of the scaling exponent $`\alpha `$.

In this Letter, we present scanning tunneling microscopy (STM) characterization of surface smoothing in 2 MeV Si+ ion irradiation of Si surfaces at normal incidence ($`\theta =0`$). At length scales below $`50`$ nm we observe smoothing of the ion-bombarded surface. The observed value of the roughness exponent $`\alpha =0.53\pm 0.02`$ indicates the self-affine nature of the smoothed surface. The ion irradiated surface shows smoother surface texture at smaller length scales. We have chosen MeV ions for which sputtering yield is small. In comparison, the collision-induced atomic displacement and effective surface diffusivity is large. Together with normal incidence, these conditions are expected to cause smoothing. The observation of scale dependent smoothing with increased smoothing at smaller length scales has direct bearing on ion beam processing of nanostructures.

Si(100) substrates were irradiated with 2.0 MeV Si+ ions in the ion implantation beam line of our 3 MV tandem Pelletron accelerator. The ion beam was incident along the surface normal ($`\theta 0`$) and rastered on the sample in order to obtain a uniformly irradiated area. One half of the sample was masked and hence unirradiated. An ion beam flux of $`1\times 10^{12}`$ cm-2 sec-1 was used with fluences in the range $`10^{15}`$ to $`10^{16}`$ ions/cm2. The samples were kept at room temperature during ion irradiation. The pressure in the chamber was $`10^7`$ mbar. The sample was then taken out of the irradiation chamber and inserted into a STM chamber (pressure: 3$`\times 10^{10}`$ mbar) with an Omicron variable temperature STM operating at room temperature. STM height calibration was done by measuring atomic step heights on clean Si(111) and Si(100) surfaces. Roughness measurements were made on the pristine and the irradiated halves of the sample. We did not remove the thin ($``$1.5 nm) native oxide from the Si surface because the surface topography may be perturbed by the effect of Ehrlich-Schwoebel barriers in different crystallographic directions on a crystalline surface. In this regard the presence of the thin oxide layer is helpful and the effect of the anisotropic diffusion can be neglected.

In order to determine the roughness exponent from STM images we follow the procedure described in ref. 6. Typical STM images from the pristine and the irradiated (fluence $`4\times 10^{15}`$ ions/cm2) parts of a sample are shown in Fig.1. A large number of scans, each of size $`L`$, were recorded on the surface at random locations. The $`\sigma `$ values for the rms roughness given by the instrument for the individual scans were then averaged. This procedure was repeated for many different sizes and a set of average $`\sigma `$ vs. $`L`$ values was obtained (each $`\overline{\sigma }`$ is the average of six to fifteen measurements). Each $`\sigma `$ value was computed after the instrument plane fitting and subtraction procedure had been carried out. $`\overline{\sigma }`$ vs. $`L`$ log-log plots for both halves of the sample are shown in Fig.2. For the ion-bombarded area of the sample we observe surface smoothing and by fitting the linear part of the data we obtain $`\alpha =0.53\pm 0.02`$ below a length scale of $``$ 50 nm, indicating the self-affine nature of the irradiated surface. Below this length scale, the pristine half of the sample shows no linear region in the log-log plot of $`\overline{\sigma }`$ vs. $`L`$. Two vertical profiles $`h(x)`$ along the lines marked in Fig.1. are shown in the inset of Fig.2. It is also clear from these profiles that for the irradiated part of the sample the surface is much smoother at shorter length scales as indicated by the roughness data and the scaling exponent.

In earlier scaling studies on ion-bombarded surfaces, the conditions of ion energy and the angle of incidence were favorable for strong sputtering and sputter-erosion of surfaces caused roughening. In order to explain the dominance of smoothing over roughening in our case let us first compare the sputtering yields. From the conditions in refs. and , we estimate the sputtering yields of 3.7 atoms/ion and 3.9 atoms/ion, respectively, using the TRIM (transport of ions through matter) code. In our case the higher ion energy and the normal incidence $``$ both contribute to lowering the sputtering yield, which is $`<`$ 0.2 atom/ion. Thus the sputtering yield is smaller by almost a factor of 20. This indicates why surface erosion, main reason for roughness enhancement, is not significant in our case. In fact at large length scales surface roughness remains unaffected by ion bombardment. On the other hand, number of surface atoms that would contribute to effective surface mobility is large as discussed below. In ion-atom collisions in solids and at the surface, the elastic energy lost by an ion is transferred to a recoil atom, which itself collides with other atoms in the solid and so forth. In this way the ion creates what is called a collision cascade. The displaced atoms in this collision cascade may acquire a kinetic energy enough to escape from the solid surface $``$ a phenomenon known as sputtering. However, if the energy (component normal to surface) of the displaced atoms is smaller than the surface binding energy, the atoms may reach the surface but cannot leave the surface. They can however drift parallel to the surface. We show the results of a TRIM simulation of sputtering yield for our case in Fig.3. This shows the atoms reaching the surface vs. their energies normal to the surface. Atoms which have energies greater than the surface binding energy ($``$ 4.7 eV) will be sputtered. However, we notice that a large number of atoms reach the surface with low energy ($`<`$ 4.7 eV) with the number of atoms/eV peaking at $``$ 1 eV. These atoms will not leave the surface (not be sputtered). The role of these atoms is important in surface smoothing. These atoms have too low an energy (normal to surface) to escape the energy barrier at the surface and will translate parallel to the surface. This collision-induced atomic displacement and the consequent effective diffusivity parallel to the surface due to ballistic atomic transport can be the dominant surface relaxation mechanism. As a result smoothing may dominate roughening as discussed later in more details. Eklund at al studied submicron-scale surface roughening induced by ion bombardment and obtained an scaling exponent $`\alpha 0.20.4`$. This value of the exponent is reasonably explained by the anisotropic KPZ equation ($`\alpha =0.38`$) when the surface diffusion term is expected to contribute negligibly. On the other hand, there are no concrete predictions of the exponents for the case where ion beam induced surface smoothing or diffusivity is dominant. Neither we know any scaling theory which predicts $`\alpha 0.5`$. Assuming the possibility that the scaling theories applicable to nonequilibrium film growth may also be applicable to ion bombardment, so long as no eroded material is redeposited onto the surface, we compare the observed exponent with those expected for the deposition process, which are $`\alpha 0.35`$ when surface mobility of the deposited particles is ignored and $`\alpha =0.66`$ when surface mobility is allowed. In the first case the exponents are in good agreement for deposition and ion bombardment. In our case surface mobility is important and the observed value of $`\alpha =0.53`$ is closer to that for the deposition model that includes surface mobility. Incidentally, Krim et al. also observed $`\alpha =0.53`$ for ion bombardment of an Fe film on a MgO substrate where roughening, rather than smoothing, was dominant.

For ion irradiation, Carter and Vishnyakov derived an expression showing the relative magnitudes of the roughening (sputtering) term and the smoothing term due to recoiled atoms which qualitatively explains the domination of smoothing over roughening at normal and near-normal ($`\theta 0`$) incidence of the ion beam. However, there is no prediction of scaling exponent. For $`\theta 0`$ they predict that smoothing dominates roughening at all wave vectors. We find that at larger length scales ($`>50`$ nm) initial surface roughness remains practically unaffected by ion bombardment while smoothing becomes increasingly dominant at lower length scales below 50 nm.

In order to show the relative strength of the smoothing and the roughening terms, Carter and Vishnyakov extended the treatment given by Bradley and Harper, who showed (in 1+1 dimension) that, due to sputter-erosion alone, the deterministic defining equation for $`h(x,t)`$ can be written as

$`{\displaystyle \frac{h}{t}}=`$ $`{\displaystyle \frac{J}{N}}Y_0(\theta ){\displaystyle \frac{J}{N}}{\displaystyle \frac{}{\theta }}\left[Y_0(\theta )cos\theta \right]{\displaystyle \frac{h}{x}}`$ (2)
$`+{\displaystyle \frac{Ja}{N}}Y_0(\theta )\mathrm{\Gamma }_1(\theta ){\displaystyle \frac{^2h}{x^2}}`$
where $`J`$ is the mean ion flux incident at angle $`\theta ,N`$ is the solid atomic density, $`Y_0(\theta )`$ is the sputtering yield of a plane surface, $`a`$ is the mean depth of energy deposition by an ion, and $`\mathrm{\Gamma }_1(\theta )`$ is a function of $`\theta `$, and standard deviations $`\alpha `$ and $`\beta `$ of the bi-Gaussian ellipsoidal ion energy spatial deposition density function. For order of magnitude estimation the ellipsoidal distribution has been approximated by a spherical distribution with $`a=\alpha =\beta `$, in which case

$$\mathrm{\Gamma }_1(\theta )=sin^2\theta \left(cos^2\theta /2\right)\left(1+sin^2\theta \right)$$
(3)
In order to introduce the effective diffusion parallel to the surface they estimated the atomic flux parallel to the surface to modify the last term in Eq.(1):

$$\frac{J}{N}\left\{f(E)dcos2\theta Y_0(\theta )a\left[sin^2\theta \frac{cos^2\theta }{2}(1+sin^2\theta )\right]\right\}\frac{{}_{}{}^{2}h}{x^2}$$
(4)
where $`f(E)`$ is the no. of recoil atoms each ion generates in the solid and $`d`$ is the average distance traveled by the recoiled atoms. $`d`$ is of the order of a few interatomic distances. For $`\theta =0`$, the expression (3) \[ i.e, the last term in Eq(1) \] is negative and smoothing dominates roughening at all wave vectors. $`f(E)=k(E)/2E_d`$, where $`k(E)`$ is the fraction of ion energy deposited in elastic collisions and $`E_d`$ is a displacement energy. In the simulation results shown in Fig.3 we have used $`E_d`$=15 eV. The results shown in Fig.3 only qualitatively shows how a large number of hyperthermal recoil atoms, arriving at the surface but unable to escape the surface, can cause surface smoothing as implied by Eq.(1) along with expression (3). Expression (3) only qualitatively describes the effect of $`f(E)`$ in surface smoothing. For a quantitative understanding future theoretical work should include the effect of a distribution like that shown in Fig.3. So far theoretical works concentrated only on the low ion energy regime where sputter-erosion is dominant and the approximation ($`a=\alpha =\beta `$) used in deriving the expression for $`\mathrm{\Gamma }_1(\theta )`$ \[ Eq.(2) \] may be valid. However, for high ion energies, it is not valid. For example in our case, for 2 MeV Si ions in Si, $`a=1.94`$ $`\mu `$m, $`\alpha =248`$ nm and $`\beta =288`$ nm. In the existing theories it is assumed that energy released by the ions at a depth $`a`$ contributes an amount of energy to surface points that may induce surface atoms to break their bonds and leave the surface. This is true for low ion energies where $`a`$ is small. However, ion energy release deep inside the sample would hardly have any effect on surface atoms. Future theories must take this aspect into account.

In conclusion, we have observed nanoscale surface smoothing in ion bombardment. The smoothed surface is a self-affine fractal surface with a scaling exponent $`\alpha =0.53\pm 0.02`$. Below a length scale of $`50`$ nm, the smoothing is more dominant at smaller length scales. This phenomenon may be used in reducing surface roughness of nanostructural devices by ion beam processing as ion beams are widely used in device fabrication. Transport in nanostructures is expected to improve when roughness is minimized. For an understanding of the scaling exponent observed in surface smoothing further theoretical studies will be necessary.

Corresponding author: bhupen@iopb.res.in

FIGURE CAPTIONS:

Fig.1. STM images recorded on a pristine (top) and ion-bombarded (bottom) silicon surface. The scan size is 300 $`\times `$ 300 nm2 and the vertical scale (black to white) is 2.2 nm. Height profiles along the lines are shown in Fig.2.

Fig.2. Average root-mean-square roughness vs. scan size on the pristine and the ion irradiated surfaces. Each point represents an average of 6 to 15 scans recorded at random locations on the surface. Surface smoothing is observed at scan sizes below $`50\times 50`$ nm2. The least-squares fit (solid line) to the linear portion of the data for the irradiated sample gives the scaling exponent $`\alpha =0.53\pm 0.02`$. No linear part is observed for the pristine sample data. Two vertical profiles $`h(x)`$ measured along the lines marked in Fig.1, are shown in the inset (scales in nm): (a) pristine, (b) irradiated sample.

Fig.3. A Monte-Carlo simulation result showing the energy distribution of ion-beam induced displaced atoms reaching the surface. Atoms with energy $`>4.7`$ eV leave the surface (sputtered). The large number of atoms below 4.7 eV (surface binding energy) cannot leave the surface and contribute to an effective surface diffusion due to ballistic atomic transport leading to smoothing.

Fig.1 (bottom)

Fig.2.

Fig.3


text:
# 𝑇-odd quark distributions: QCD versus chiral models

## Abstract

The $`T`$-odd quark distrbutions are shown to vanish in the chiral sigma model in contrast to the opposite widespread opinion. This failure of the chiral sigma model is a feature of the model itself and has nothing to do with the recent progress in the clarification of the status of the $`T`$-odd distributions in QCD.

Recently certain progress has been made BHS-02 ; Collins-02 ; JY-02 ; BJY-02  in understanding the properties of the so-called $`T`$-odd quark distributions in QCD Sivers-90 ; ABM-95 ; BM-98 ; BM-99 . Naively the $`T`$-odd quark distribution functions in nucleon are defined as follows (notation of Ref. ABDM-02 ):

$$\mathrm{\Delta }_0^Tf(x,𝐤_{}^2)=\mathrm{Im}\frac{dy^{}d^2𝐲_{}}{2(2\pi )^3}e^{ixP^+y^{}+i𝐤_{}𝐲_{}}$$

$$\times P,|\overline{\psi }(0,y^{},𝐲_{})\gamma ^+\psi (0)|P,+,$$
(1)
$$\mathrm{\Delta }_T^0f(x,𝐤_{}^2)=\frac{dy^{}d^2𝐲_{}}{2(2\pi )^3}e^{ixP^+y^{}+i𝐤_{}𝐲_{}}$$

$$\times P|\overline{\psi }(0,y^{},𝐲_{})i\sigma ^{2+}\gamma _5\psi (0)|P.$$
(2)
Here $`\psi `$ is the quark field and $`|P,\pm `$is the nucleon state with the given momentum and helicity.

If one uses this *naive definition* of functions $`\mathrm{\Delta }_0^Tf`$ (introduced by Sivers Sivers-90 ) and $`\mathrm{\Delta }_T^0f`$ (studied in Ref. BM-98 ), then *in any theory* with nonbroken $`T`$-invariance these functions should vanish. However, as noticed in Ref. Collins-02  the *correct QCD definition* of functions $`\mathrm{\Delta }_0^Tf`$, $`\mathrm{\Delta }_T^0f`$ should include properly chosen Wilson lines. Moreover, the direction of the Wilson lines is sensitive to the considered hard process. The appearance of the Wilson lines in the definition of functions $`\mathrm{\Delta }_0^Tf`$, $`\mathrm{\Delta }_T^0f`$ breaks the naive $`T`$-invariance argument so that one can have nonvanishing functions $`\mathrm{\Delta }_0^Tf`$, $`\mathrm{\Delta }_T^0f`$ in QCD.

In Ref. ABDM-02  an attempt was made to study $`T`$-odd distributions in the chiral model with the lagrangian

$$L=i\overline{\psi }\gamma ^\mu _\mu \psi g\overline{\psi }\left(\sigma +i\gamma _5\tau ^a\pi ^a\right)\psi $$

$$+\frac{1}{2}(_\mu \sigma )^2+\frac{1}{2}(_\mu \pi ^a)^2U(\sigma ,\pi )$$
(3)
with the mexican-hat potential $`U`$. Actually most of the arguments of Ref. ABDM-02  and of this paper are valid for a wider class of models with the action

$$S[\psi ,\sigma ,\pi ]=\underset{c=1}{\overset{N_c}{}}d^4x$$

$$\times \left[i\overline{\psi }_c\gamma ^\mu _\mu \psi _cg\overline{\psi }_c\left(\sigma +i\gamma _5\tau ^a\pi ^a\right)\psi _c\right]+N_cS_1(\sigma ,\pi )$$
(4)
containing a chiral invariant mesonic piece $`S_1(\sigma ,\pi )`$. We have explicitly written the summation over the quark color indices $`c`$ and have inserted the number of colors $`N_c`$ in front of the bosonic part of the action in the form which allows a systematic $`1/N_c`$ expansion in the large-$`N_c`$ limit. It is assumed that the parameters of the chiral model are chosen so that the chiral $`SU(2)_LSU(2)_RU(1)_V`$ symmetry is spontaneously broken to the vector symmetry group $`U(2)_V`$ and that the model contains states with the nucleon quantum numbers. We also assume that the action $`S_1(\sigma ,\pi )`$ is $`C`$, $`P`$ and $`T`$ invariant.

The main conclusion of Ref. ABDM-02  is that the $`T`$-odd quark distribution functions computed in the model (3) might be nonvanishing. This result is rather strange. Indeed, the lagrangian (3) is explicitly $`T`$-invariant, so that the theory is invariant under the *standard* time reflection (STR) transformations of the quantum fields:

$`T\psi (x^0,𝐱)T^1`$ $`=`$ $`C\gamma _5\psi (x^0,𝐱),`$
$`T\sigma (x^0,𝐱)T^1`$ $`=`$ $`\sigma (x^0,𝐱),`$ (5)
$`T\pi ^a(x^0,𝐱)T^1`$ $`=`$ $`\epsilon ^a\pi ^a(x^0,𝐱),`$
$$\epsilon ^1=\epsilon ^3=1,\epsilon ^2=1.$$
(6)
Under the following conditions:

* the chiral model makes sense (the ground state is stable and the states with the quantum numbers of nucleon exist),
* no exotic phenomena happen like spontaneous breakdown of the $`T`$-invariance,
* one uses the naive definition of the $`T`$-odd functions (1), (2),

the STR symmetry $`T`$ (5) of the model (3) automatically leads to the conclusion that the $`T`$-odd functions (1), (2) vanish in this model in contrast to the argument of Ref. ABDM-02 .

Let us try to understand what stands behind this discrepancy. First one notices that the authors of Ref. ABDM-02  actually do not compute the $`T`$-odd parton distributions. Their conclusion about nonvanishing $`T`$-odd distributions is based only on certain symmetry arguments.

Instead of the STR transformation $`T`$ (5) the authors of Ref. ABDM-02  suggest to consider what they call “nonstandard time reversal” (NSTR) transformation $`T_{\mathrm{NSTR}}`$ which differs from STR $`T`$ simply by the isotopic 180 rotation:

$$T_{\mathrm{NSTR}}=e^{i\pi I_2}T$$
(7)
so that

$`T_{\mathrm{NSTR}}\psi (x^0,𝐱)T_{\mathrm{NSTR}}^1`$ $`=`$ $`i\tau _2C\gamma _5\psi (x^0,𝐱),`$
$`T_{\mathrm{NSTR}}\sigma (x^0,𝐱)T_{\mathrm{NSTR}}^1`$ $`=`$ $`\sigma (x^0,𝐱),`$ (8)
$`T_{\mathrm{NSTR}}\pi ^a(x^0,𝐱)T_{\mathrm{NSTR}}^1`$ $`=`$ $`\pi ^a(x^0,𝐱).`$
Since the isotopic invariance is a symmetry of the model (3), it does not matter whether one works with STR or NSTR. But one should not forget that NSTR is not the only symmetry of the theory, there is also the isotopic symmetry. However, the authors of Ref. ABDM-02  concentrate on NSTR and ignore other symmetries. They make the correct observation that NSTR taken alone allows nonzero $`T`$-odd functions and they declare without any calculation that the sigma model may lead to nonzero $`T`$-odd distributions. The fault of this argument is that NSTR is not the only symmetry of the model — if one combines NSTR with the isotopic invariance to recover STR, then one can use the standard argument to show that the $`T`$-odd functions vanish.

One could wonder, why in Ref. ABDM-02  such a significance is attributed to NSTR, whereas the isotopic invariance and STR are ignored, although they are exact symmetries of the model. It seems that the authors mix two different issues:

1) the symmetries of the action which are not spontaneously broken (these symmetries include both STR and NSTR as well as the isotopic invariance),

2) the symmetries of the mean field describing the nucleon states.

As it is well known, the chiral quark-soliton models of type (4) allow a systematic $`1/N_c`$ expansion in the limit of the large number of colors $`N_c`$ Witten-79 ; GS-84 . The leading order of this $`1/N_c`$ expansion is nothing else but the mean field approximation. In terms of the path integral approach, the $`1/N_c`$ expansion can be constructed as a systematic expansion about the corresponding saddle point of the effective action, which can be obtained by the Gaussian integration over the quark fields. The properties of the vacuum, in particular the spontaneous breakdown of the chiral symmetry, are controlled by the constant (space and time independent) saddle point $`\sigma _0,\pi _0^a`$. As for the nucleon states, they are described by the $`𝐱`$ dependent (but static) saddle point usually referred to as “soliton”. This soliton field has the form

$`\sigma _{\mathrm{sol}}(𝐱)`$ $`=f_1(|𝐱|),`$
$`\pi _{\mathrm{sol}}^a(𝐱)`$ $`=x^af_2(|𝐱|),`$ (9)
which violates a lot of symmetries: spatial translations and $`O(3)`$ rotations, $`SU(2)`$ isotopic rotations, STR. On the other hand, the soliton field is invariant with respect to combined space-isotopic rotations and with respect to NSTR. Probably due to this privileged status of NSTR, the authors of Ref. ABDM-02  have concentrated on this symmetry ignoring the others. But the symmetries of the soliton field have nothing to do with the symmetries of the Hamiltonian and of the ground state. For example, the soliton field is not translationally invariant but this does not mean that the translational invariance of the model is broken.

Thus both STR and NSTR are exact symmetries of the model. Another way to understand that STR is an exact symmetry is to remember that according to Eq. (7) NSTR differs from STR by an isotopic rotation. If STR were broken (spontaneously or in some other way) with NSTR remaining unbroken, then this would automatically lead to the broken isotopic symmetry. The third argument in favor of STR comes from the *CPT* theorem: the breakdown of STR would mean that the $`CP`$ symmetry is violated in the model (3), which is probably not the aim of the authors of Ref. ABDM-02 .

The illusion of the violation of STR by the soliton field actually disappears if one properly uses the standard methods of the treatment of collective coordinates in the semiclassical quantization GS-84 ; GS-7x . Using this technique, one can in principle check the STR symmetry explicitly in any order of the $`1/N_c`$ expansion.

To summarize, as long as the isotopic $`SU(2)`$ invariance is not broken, there is no difference whether one uses STR or NSTR — anyway the $`T`$-odd distribution functions (1), (2) vanish in the model considered in Ref. ABDM-02  in contrast to the conclusion of the authors.

From the above analysis it is obvious that the vanishing of $`T`$-odd distributions (1), (2) in the model (3) is an internal feature of this model and has nothing to do with the nonvanishing $`T`$-odd distributions in QCD.

The reason allowing the existence of the $`T`$-odd parton distributions in QCD is the nontrivial role played by the direction of the Wilson lines accompanying the quark fields Collins-02 . If one really wants to model the QCD effects, then one has to find a way to incorporate the gauge links into the effective model.

On general grounds it is expected that the large-$`N_c`$ QCD in the color singlet sector should be equivalent to some effective theory with (generally bilocal) meson fields Witten-79 . However, the corresponding action is not known even in the leading order of the $`1/N_c`$ expansion. Instead one usually deals with ad hoc models like the sigma model (3). In the best case these models can reproduce some features of QCD like spontaneous breakdown of chiral symmetry, correct large-$`N_c`$ counting, certain pieces of the hadronic spectrum etc. Since these models are not derived from QCD, the question how to implement the Wilson lines in these models can hardly be treated seriously.

An interesting attempt to describe the gluonic degrees of freedom within a low energy effective model was made in the framework of the instanton vacuum model DP-86 . Although the “derivation” of this model from QCD involves quite a number of approximations and simplifications, in principle one can trace the relation between the gluonic degrees of freedom of QCD and the quark-pion fields of the effective theory DPW . It would be interesting whether within this approach one could say something about the $`T`$-odd quark distribution functions.

Acknowledgement. The author appreciates discussions with Ya. I. Azimov, D. Boer, A. Metz, P. J. Mulders, M. V. Polyakov and P. Schweitzer.


text:
# ITEP-LAT/2002-24 KANAZAWA 02-33 The flux distribution of the three quark system in SU(3)

## 1 Introduction

The SU(3) lattice study on the three quark (3Q) system is important for the clarification of the baryonic structure. For more than 20 years, there is a crucial question about the long range force of the 3Q system : whether there is a genuine three body force, or it can be described by the sum of the two body forces. In the former case, the flux structure is expected to be of “$`Y`$-shape”, which has a junction at the point where the total length of strings from a quark to the junction is minimal. On the other hand, in the latter case, the flux structure is expected to be of “$`\mathrm{\Delta }`$-shape”, which consists of the 3 sets of two body interactions.

In order to distinguish between the two $`Ans\ddot{a}tze`$, the 3Q potential was investigated using lattice simulations by several groups more than 15 years ago and also recently. The question, however, hasn’t been settled yet. Some authors support the $`\mathrm{\Delta }Ansatz`$, others find evidence for the $`YAnsatz`$  (see also ). This discrepancy is due to the difficulty of measuring the 3Q potential accurately for a large quark separation and moreover the difference between the two $`Ans\ddot{a}tze`$ in the 3Q potential is not large, at most 20$`\%`$. Therefore, it is useful to measure the flux directly using lattice simulations.

A couple of groups measured the flux in multi-quark systems. However, for example, the Helsinki and the UK groups  investigated them only in SU(2), where mesons and baryons have no differences, and no one has succeeded in the direct measurement of flux in the 3Q system in SU(3), since it needs huge statistics and large computer memory.

We shall study the flux of the three quark system in SU(3) after fixing the gauge into the maximally abelian (MA) gauge based on the dual superconductor scenario, because it is easy to interpret the flux structure and the validity of this scenario for the flux tube has been confirmed by Bali et al.. In addition, one can get better signals on the abelian flux in the MA gauge than on the full flux in SU(3) and its computation needs only less statistics and memory than the full one, which means that there is the possibility to succeed in the measurement of the flux without waiting for a new bigger supercomputer.

## 2 Operators and simulation details

The local quantities describing the abelian flux of the 3Q system are obtained from the correlation between an appropriate operator and the 3Q abelian Wilson loop:

$`O(s)={\displaystyle \frac{O(s)W_{3\mathrm{Q}}}{W_{3\mathrm{Q}}}}O,`$
where the 3Q abelian Wilson loop is defined by

$`W_{3\mathrm{Q}}={\displaystyle \frac{1}{3!}}|ϵ_{abc}|u_\mathrm{R}^au_\mathrm{G}^bu_\mathrm{B}^c`$
with a path product

$`u_𝒞^a={\displaystyle \underset{s\mathrm{\Gamma }_𝒞}{}}u_\mu ^a(s)\mathrm{\Gamma }_𝒞=\mathrm{\Gamma }_\mathrm{R},\mathrm{\Gamma }_\mathrm{G},\mathrm{\Gamma }_\mathrm{B}`$
of abelian link variables along the path $`\mathrm{\Gamma }_𝒞`$ in fig. 1. Differently from the nonabelian case, the color of quarks does not change during the propagation, because the off-diagonal components of the gauge field are frozen in the abelian gauge.

The simulation is performed in quenched QCD at $`\beta `$=6.0 on $`16^332`$ lattices. The link variables are fixed into the maximally abelian gauge with a simulated annealing algorithm. The 3 static charges, red, green and blue quarks, are located at $`(17,14)`$ $`(22,6)`$ and $`(12,6)`$ in the $`x`$-$`y`$ plane, respectively. For noise reduction, we use the smearing method for spatial links of the 3Q abelian Wilson loop.

## 3 Flux of the 3 quark system

Figs. 2 and 3 show the color electric field and the magnetic super current in the 3Q system. For example, the green component of the electric field flows from the green quark into the blue and red quarks. For each flux, we observe a solenoidal monopole current in the plane perpendicular to the flux. Corresponding to the amount of flux of the electric field, there is a magnetic super current around each flux tube. Our results are in a qualitative agreement with dual Ginzburg-Landau theory .

In order to clarify the flux structure, we decompose the abelian gauge field into the monopole part and photon part. Fig. 4 shows the action density in the abelian, monopole and photon parts, respectively. The 3 peaks originating from the sources appear in the photon part, while the plateau, which is responsible for a confinement potential, appears in the monopole part. There does not appear a dent, which is expected in the case of $`\mathrm{\Delta }Ansatz`$, in the center of the plateau. Then, the abelian action density has a clear junction in the center of the 3 quarks, which suggests $`YAnsatz`$.

As shown in fig. 5, in the monopole part, there appear two rotating induced electric fields, which cancel the electric field by original sources outside of the flux. Thus, we can expect the following scenario: When the three quark sources are put into the vacuum, the solenoidal monopole current appears for each flux, and then they induce the two rotating electric fields such that the flux is squeezed into the “$`Y`$-shape”.

Moreover, we investigated the flux of the 3Q system also in the full QCD generated with the non-perturbatively $`O(a)`$ improved Wilson fermion action at $`\beta `$=5.29 and $`\kappa `$=0.1355. For $`m_\pi /m_\rho 0.7`$, the flux structure is almost the same and we did not observe any big effect of dynamical quarks, like in the case of the quark-antiquark flux tube.

This work is partially supported by grant INTAS-00-00111. H. I. thanks Humboldt University for hospitality. H. I. and V. B. are supported by JSPS. Computations have been done on COMPAQ AlphaServer ES40 at Humboldt University.


text:
# Methyl group dynamics in a confined glass

## I Introduction

Porous materials have attracted great interest in the last years, since they provide a framework for the study of the physical propierties of solids under restricted geometry . The dynamics of systems confined in these materials is expected to differ considerably from that in the bulk state, first of all because the length of the intermolecular interactions and correlations is limited by the finite size of the pore if the latter is sufficiently small. Moreover, the interactions between the pore wall and the confined molecules -which strongly depend on the porous system- would affect the dynamics of these molecules in a non-well defined surface layer on the confined system. Due to their simple dynamics , small molecular rotors as ammonia, methane, ammine ions or methyl groups can be considered as adequate dynamical probes of the effects of confinement by comparing the rotational barriers in the bulk and in confined geometry .

A thorough neutron scattering study on methyl groups dynamics in bulk glassy toluene has been recently presented by some of us . The experimental results were successfully analysed in terms of a well-established barrier distribution model, such a distribution resulting from the structural disorder present in the glass. The average barrier of the distribution was found to be notably higher than that of the $`\beta `$-crystalline reference phase, which shows a short-range structure similar to that of the glassy phase . This fact could indicate that the relevant intermolecular interactions contributing to the rotational potential of the methyl groups extend beyond nearest-neighbours.

In order to make a comparison with the bulk state and to shed new light about this question and about the interaction lengths controlling the parameters of the distribution, we have carried out neutron scattering measurements on glassy toluene confined in mesoporous silicates with different pore sizes.

## II Theoretical aspects

The usual single-particle model for methyl group dynamics in crystalline systems at very low temperatures ($`T`$ 2 K) is that of a rigid rotor tunnelling through a one-dimensional rotational potential $`V(\mathrm{\Phi })`$, which is restricted to keep the rotational symmetry of the methyl group . In most of cases only the leading threefold term of the Fourier expansion is retained, higher order corrections being small: $`V(\mathrm{\Phi })=V_3(1\mathrm{cos}3\mathrm{\Phi })/2`$. The torsional levels are tunnel-split by the overlapping of the neighbouring single-well wavefunctions, and the energy splitting of the ground torsional level, or rotational frequency $`\omega _\mathrm{t}`$, can be observed in neutron scattering spectra -for barrier heights $`V_3`$ below $``$ 700 K- in the $`\mu `$eV-range as two resolution-width inelastic peaks centered at $`\pm \omega _\mathrm{t}`$. The corresponding incoherent scattering function for rotational tunnelling, normalized to scattering from one hydrogen is :

$`S_{\mathrm{MG}}^{\mathrm{inc}}(Q,\omega )=[A(Q)+B(Q)]\delta (\omega )`$ (1)
$`+B(Q)[\delta (\omega +\omega _\mathrm{t})+\delta (\omega \omega _\mathrm{t})]`$ (2)
with $`\mathrm{}Q`$, $`\mathrm{}\omega `$ respectively the momentum and energy transfer of the neutron, $`A(Q)=[1+2j_0(Qr)]/3`$, $`B(Q)=2[1j_0(Qr)]/9`$, $`j_0(Qr)=\mathrm{sin}(Qr)/(Qr)`$ and $`r`$ the H-H distance in the methyl group.

At temperatures typically above $``$ 50 K, a picture of thermally activated classical hopping over the barrier is valid. The corresponding hopping rate follows an Arrhenius-like temperature dependence $`\mathrm{\Gamma }=\mathrm{\Gamma }_{\mathrm{}}\mathrm{exp}(E_\mathrm{A}/kT)`$, with $`E_\mathrm{A}`$ the classical activation energy, defined as the difference between the top of the barrier and the ground torsional state. $`\mathrm{\Gamma }_{\mathrm{}}`$ is a temperature-independent preexponential factor. The incoherent scattering function for the classical hopping regime is :

$$S_{\mathrm{MG}}^{\mathrm{inc}}(Q,\omega )=A(Q)\delta (\omega )+3B(Q)L(\omega ;\mathrm{\Gamma })$$
(3)
with $`L(\omega ;\mathrm{\Gamma })`$ a normalized Lorentzian of half-width at half-maximum (HWHM) equal to $`\mathrm{\Gamma }`$.

A series of investigations in structural glasses has evidenced the need of introducing a barrier distribution $`g(V_3)`$ to give account for the different features observed in the spectra of these systems (rotation-rate-distribution-model, RRDM) . Such a distribution has its origin in the different local environments for the methyl groups resulting from the structural disorder in the glassy state. According to this approximation, the spectra for the glassy system are obtained as a superposition of crystal-like spectra (tunnel-like as (1) at $`T`$ 2 K, or classical-like as (2) at high temperature) weighted by $`g(V_3)`$:

$$S_{\mathrm{MG}}^{\mathrm{inc}}(Q,\omega )=S_{\mathrm{MG}}^{\mathrm{inc}}(Q,\omega ,V_3)𝑑V_3$$
(4)
Numerical relationships  between $`V_3`$, $`\mathrm{}\omega _\mathrm{t}`$ and $`E_\mathrm{A}`$ allow the transformation between the corresponding distributions of these quantities $`g(V_3)`$, $`h(\mathrm{}\omega _\mathrm{t})`$ and $`f(E_\mathrm{A})`$. By assuming that the preexponential factor $`\mathrm{\Gamma }_{\mathrm{}}`$ is $`V_3`$-independent , $`f(E_\mathrm{A})`$ can also be transformed, through the Arrhenius law for the temperature dependence of the hopping rate, into the corresponding distribution of hopping rates $`H(\mathrm{log}\mathrm{\Gamma })`$. Thus, the general equation (3) for the RRDM is reduced to the tunnelling and hopping cases by doing the substitutions $`\delta (\omega \pm \omega _\mathrm{t})h(\pm \mathrm{}\omega _\mathrm{t})d(\mathrm{}\omega _\mathrm{t})`$ in equation (1) and $`L(\omega ,\mathrm{\Gamma })H(\mathrm{log}\mathrm{\Gamma })L(\omega ,\mathrm{\Gamma })d(\mathrm{log}\mathrm{\Gamma })`$ in equation (2). The consistency of the model is reached when the distributions of tunnelling frequencies and classical hopping rates, obtained respectively at low and high temperature, follow from the same barrier distribution $`g(V_3)`$.

## III Experimental details

Neutron scattering measurements were carried out at the backscattering spectrometer IN16 of the Institute Laue Langevin (ILL, Grenoble, France). A wavelength of 6.27 Å was selected providing an energy resolution of 0.5 $`\mu `$eV (HWHM). The energy window available by Doppler shift covered from -15 to 15 $`\mu `$eV. The scattering angle covered a range from 11 to 149o, resulting in a $`Q`$-window from 0.2 to 1.9 Å-1. The instrumental resolution was calibrated by a vanadium sample, which shows purely elastic scattering. Raw data were corrected for detector efficiency, scattering from the sample holder and absorption by means of ILL standard programs.

Spectra were taken on three samples of ring-deuterated toluene confined in mesoporous silicates of the series MCM-41 and SBA-15, which geometry consists of a honeycomb-type lattice of parallel cylindrical pores. They were synthesized with pore diameters $`D`$ = 24, 35 and 47 Å, according to an hydrothermal procedure described in Ref. . The porous geometry parameters were confirmed by neutron diffraction and nitrogen absorption experiments. Deuteration of the silanol groups was achieved by chemical H/D exchange with no further chemical treatment of the surface. Complete filling of the outgased matrices was achieved with an appropriate mass of toluene .

Flat samples of thickness 2 mm were used to get transmissions above 90 %, allowing to neglect multiple scattering effects in the $`Q`$-range investigated, 1.3 - 1.9 $`\mathrm{\AA }^1`$, where the ratio $`B(Q)/A(Q)`$ is most favourable in the $`Q`$-window of IN16. The cooling rate from room temperature to 2 K achieved in the cryostat was sufficient to get the glassy state of confined toluene. Spectra were also taken for the empty matrixes, showing purely elastic scattering in the investigated temperature interval.

## IV Results

Figs. 1,2 show the incoherent scattering functions, normalized to maximum unity, for confined and bulk glassy toluene -the latter taken from the previuos investigation of Ref. -, obtained after substraction of all the elastic contributions due to coherent scattering and to incoherent scattering from other atoms different from the methyl group protons. In this way the comparison between the dynamics in the confined and in the bulk state is not affected by the different ratios from the methyl group incoherent cross section to the total one. From a direct visual inspection, it is evidenced that confinement results in a strong change in the tunnelling dynamics at 2 K (Fig. 1), and also at the hopping dynamics at high temperature (Fig. 2). As previously reported for the bulk state , the experimental data for the confined samples at 2 K could be reproduced in terms of a superposition of tunnelling frequencies, as exposed in Section 2, such a distribution following from an asymmetric Gamma-distribution of rotational barriers $`g(V_3)`$:

$$g(V_3)=\frac{(p/\mathrm{e})^p}{V_{3_0}\mathrm{\Gamma }(p)}\left(\frac{V_3}{V_{3_0}}\right)^{p1}\mathrm{exp}\left[p\frac{V_3V_{3_0}}{V_{3_0}}\right]$$
(5)
with $`\mathrm{\Gamma }(p)`$ the Euler Gamma function, $`V_{3_0}`$ the average barrier and $`p`$ an adimensional shape parameter measuring the asymmetry of the function. The corresponding standard deviation $`\sigma _V`$ is obtained as $`\sigma _V=V_{3_0}/\sqrt{p}`$. The thick lines in Fig.1 are the RRDM theoretical functions for rotational tunnelling after convolution with the instrumental resolution. The obtained parameters for the three pore diameters investigated are given in Table 1. The error bars are $``$ 5%. The resulting distributions are of similar width, but much more asymmetric and with a notably lower average barrier than that of the bulk. A systematic increase of the average barrier is obtained with growing pore size, though far from the bulk limit in all cases, suggesting an important contribution of surface effects to the barrier distribution.

In order to check the consistency of the RRDM for the obtained parameters of $`g(V_3)`$, we derived the distribution of classical hopping rates through Arrhenius relationship from the distribution of classical activation energies derived from $`g(V_3)`$. The corresponding RRDM theoretical function for classical hopping was constructed for different values of $`\mathrm{\Gamma }_{\mathrm{}}`$, convoluted with the instrumental resolution and compared with the experimental data. Thick lines in Fig. 2 correspond to the values $`\mathrm{\Gamma }_{\mathrm{}}=`$ 9 and 8 meV for confined and for bulk glassy toluene  respectively. The error bar for $`\mathrm{\Gamma }_{\mathrm{}}`$ can be estimated in the order of $`\pm `$ 3 and $`\pm `$ 2 meV for the confined and the bulk cases respectively.

## V Discussion

The single-particle potential for methyl group rotation in a given system is built up by intra- and intermolecular interactions extending until a certain radius $`r_\mathrm{c}`$. The results above described show that the distribution of potential barriers for methyl group rotation in glassy toluene is strongly modified by confinement. This has to be certainly related with the presence of the pore-wall and the interactions between the wall and the methyl groups. In a rough approximation we can expect that the methyl groups which are at a distance from the pore wall $`r<r_\mathrm{c}`$ will feel the wall and their single-particle potential will become modified. On the other hand, methyl groups at $`r>r_\mathrm{c}`$ from the pore wall will in principle behave as bulk-like methyl groups. In this framework, $`g(V_3)`$ can be separated into a “bulk-like” and a “surface-like” contribution, corresponding, respectively, to molecules at a distance of the pore surface beyond and below $`r_\mathrm{c}`$. Previous results corresponding to dynamics of small rotors in confined crystalline systems seem to support this kind of separation . From the cylindrical geometry of the pores, and assuming an homogeneus density of methyl groups within the pores, it is straightforward to calculate the fraction of surface-like molecules as $`f_\mathrm{S}=1(12r_\mathrm{c}/D)^2`$. The corresponding surface-like contribution to the total barrier distribution can be calculated as $`g_\mathrm{S}(V_3)=[g(V_3)(1f_\mathrm{S})g_\mathrm{B}(V_3)]/f_\mathrm{S}`$, with $`g(V_3)`$ and $`g_\mathrm{B}(V_3)`$ the distributions for confined and bulk toluene defined by the parameters given in Table 1. For $`r_\mathrm{c}`$ below a certain value, this procedure yields to unphysical negative values of $`g_\mathrm{S}(V_3)`$, so that value can be understood as a lower limit of $`r_\mathrm{c}`$. Unphysical distributions were obtained for $`r_\mathrm{c}<`$ 5, 6 and 7 $`\mathrm{\AA }`$ for $`D`$ = 24, 35 and 47 $`\mathrm{\AA }`$ respectively. Thereby these values correspond to the maximum possible bulk-like contribution to $`g(V_3)`$ in each case. Fig. 3 shows the distributions for one of the pore diameters and for the bulk case, and also the surface contributions for different values of $`r_\mathrm{c}`$, including some of those with unphysical negative values due to a underestimated value of $`r_\mathrm{c}`$ (see inset). As the interaction radius $`r_\mathrm{c}`$ should not depend on the pore size and $`r_\mathrm{c}=`$7 Å is also compatible with the $`g(V_3)`$ obtained for the pore diameters $`D=`$ 24 and 35 Å (see Table II), we can take 7 Å as a first estimation for $`r_\mathrm{c}`$ in glassy toluene. It is noteworhthy that this value has to be considered only as a lower limit, because larger values could also be compatible with the different $`g(V_3)`$ obtained. Measurements at larger pore sizes would be required to give a more accurate estimation of $`r_\mathrm{c}`$. In any case our results seem to support the idea that the relevant interactions on the methyl group extend beyond the typical nearest-neighbour distance between centers-of-mass (4.7 Å) , as was suggested by the comparison of the results in the bulk glassy and $`\beta `$-crystalline phases .