text:
# Static three quark potential in the quenched lattice QCD

## 1 Introduction

Determination of the static quark-antiquark ($`Q`$-$`\overline{Q}`$) potential has been one of most vital subjects of the lattice QCD simulation . On the other hand, the status of the three quark (3$`Q`$) potential is much less settled, in spite of its importance in the baryon spectroscopy. In this paper, we present our recent result of lattice calculation of the static three quark potential .

It has been widely used the static $`Q`$-$`\overline{Q}`$ potential expressed as the sum of the short range Coulomb term and long range confining linear term. In the case of the 3$`Q`$ potential, two distinct treatments of the linear term are frequently applied. In this paper they are denoted as “$`\mathrm{\Delta }`$-type” and “$`Y`$-type” forms. The former expresses the linear potential as the sum of two-body potentials, which are proportional to the distance between each pair of quarks. The latter treats the potential with the string linking three quarks at the point with which the total length of the string takes the minimum length, $`L_{min}`$. $`Y`$-type ansatz has been used in the phenomenological studies with the string tension $`\sigma _{3Q}\sigma _{Q\overline{Q}}`$, and successfully describes the gross feature of the baryon spectrum . $`\mathrm{\Delta }`$-type ansatz is, however, also used in some model calculations with the string tension $`\sigma _\mathrm{\Delta }`$ effectively around half of $`\sigma _{Q\overline{Q}}`$. The lattice simulations have been done mainly in 1980’s , and there is one recent preliminary work . Among them, several authors interpreted their results as the support of the $`\mathrm{\Delta }`$-type picture . This embarrassing situation is mainly caused by the insufficient extraction of the ground state contribution from the baryonic Wilson loop. In nineties, calculation of $`Q`$-$`\overline{Q}`$ potential became much reliable by making use of the smearing technique which enhances the ground state contribution to the Wilson loop . This procedure is directly applicable to the three quark system, and is expected to enable us reliable extraction of the ground state contribution.

Our first goal of this work is an extraction of reliable static 3$`Q`$ potential using the smearing technique. With the lattice result, we discuss preferable description of the potential. These subjects are described in the next section.

Another important subject on the 3$`Q`$ potential is the thermal effect on it at finite temperature. Again the lattice computation is much retarded compared with the $`Q`$-$`\overline{Q}`$ potential. We present our preliminary result at $`T>T_c`$ in Section 3.

## 2 Zero temperature result

The three quark potential is extracted from the baryonic Wilson loop,

$$W_{3Q}(\stackrel{}{r}_1,\stackrel{}{r}_2,\stackrel{}{r}_3,T)\frac{1}{3!}ϵ_{abc}ϵ_{a^{}b^{}c^{}}U_1^{aa^{}}U_2^{bb^{}}U_3^{cc^{}},$$
(1)
$$U_kP\mathrm{exp}\{ig_{\mathrm{\Gamma }_k}𝑑x_\mu A_\mu (x)\},$$
(2)
where $`P`$ denotes the path-ordered product along the path $`\mathrm{\Gamma }_k`$ ($`k=`$1, 2, 3) as shown in Figure 1. The position of the contraction of three lines $`U_k`$ can be chosen arbitrarily, and the result should not depend on this point. We choose the contracting point as the spatial origin, and locate three static quarks on each three spatial axis. The same constellation is averaged making use of the translational invariance.

The numerical simulation is done on the lattice of size $`12^3\times 24`$, with the standard Wilson gauge action at $`\beta =5.7`$ in the quenched approximation. The lattice spacing is determined as $`a0.19`$ fm by setting the string tension of $`Q`$-$`\overline{Q}`$ potential to be $`\sigma =0.89`$ GeV/fm. The (baryonic) Wilson loop is measured on 210 configurations with applying the smearing technique to enhance the ground state contribution. The potentials are extracted from the (baryonic) Wilson loop by the single exponential fit in appropriate fit regions, where the Wilson loop is dominated by the ground state contribution .

Let us start with the $`Q`$-$`\overline{Q}`$ potential. It is widely accepted that the lattice result of the static $`Q`$-$`\overline{Q}`$ potential is well described by the form

$$V(r)=C\frac{A}{r}+\sigma r.$$
(3)
In Table 1, we list the results of the fit in three kinds of fit ranges: on-axis fit, fits with off-axis data in the whole $`r`$ region and at $`r2`$. Since our lattice spacing is rather coarse, the fit with the off-axis data may suffer from the rotational symmetry breaking effect. In this sense, the fit with the on-axis data is preferable to determine the parameters. For the 3$`Q`$ potential, however, we inevitably treat the off-axis constellations, and then it would be meaningful to examine the typical size of the rotational symmetry breaking effect in the $`Q`$-$`\overline{Q}`$ system. For this purpose, we carry out the fit including the off-axis data with two fit ranges. The fit using all data (all $`r`$ in Table 1) results in large $`\chi ^2`$ value, and this is mainly caused by the result at the short distance region where the statistical error is smallest and the cutoff effect is most severe. We alternatively fit at $`r2`$, which is more appropriate for comparison with the result of the 3$`Q`$ potential. In this case the fit results in $`\chi ^2/N_{df}=3.3`$, and we regard this value as a guide of the rotational symmetry breaking effect on this lattice. In both cases, the fit result of the Coulomb coefficient much differs from the value of on-axis fit than the case of the string tension.

Figure 2 shows the lattice result of the 3$`Q`$ potential from the viewpoint of the $`Y`$-type ansatz. The horizontal axis, $`L_{min}`$, is the minimum length of the flux lines linking three quarks. Three flux connecting the junction with quarks take the angles $`2\pi /3`$, which balances three forces with equal strength, in the case that all three angles of the “quark triangle” are less than $`2\pi /3`$ . The data are fitted to the form

$`V(\stackrel{}{r}_1,\stackrel{}{r}_2,\stackrel{}{r}_3)`$ $`=`$ $`C_YA_Y\left({\displaystyle \frac{1}{r_{12}}}+{\displaystyle \frac{1}{r_{23}}}+{\displaystyle \frac{1}{r_{31}}}\right)`$ (4)
$`+\sigma _YL_{min}.`$
The result is listed in Table 1. This form of the potential well describes the observed result, and the string tension is almost same as in the case of $`Q`$-$`\overline{Q}`$ potential. The $`\chi ^2`$ of the fit is in the same order as in $`Q`$-$`\overline{Q}`$ with fit range $`r2`$. Considering the rotational symmetry violation effect at this $`\beta `$, this seems to be a consistent result.

Alternative ansatz, $`\mathrm{\Delta }`$-type, is represented with the fitting function

$`V(\stackrel{}{r}_1,\stackrel{}{r}_2,\stackrel{}{r}_3)`$ $`=`$ $`C_\mathrm{\Delta }A_\mathrm{\Delta }\left({\displaystyle \frac{1}{r_{12}}}+{\displaystyle \frac{1}{r_{23}}}+{\displaystyle \frac{1}{r_{31}}}\right)`$ (5)
$`+\sigma _\mathrm{\Delta }(r_{12}+r_{23}+r_{31}),`$
and the obtained parameters are also listed in Table 1. As the result of fit, we find that the string tension is much smaller than the $`Y`$-type fit, although the constant and the coefficient of the Coulomb term are almost the same as the $`Y`$-type results. In the fit, $`\sigma _\mathrm{\Delta }`$ is slightly larger than the half of $`\sigma _Y`$. Geometrically, $`2L_{min}L_\mathrm{\Delta }=r_{12}+r_{23}+r_{31}`$ roughly holds in most cases. This explains that $`\sigma _\mathrm{\Delta }\frac{1}{2}\sigma _{Q\overline{Q}}`$. There is another possible physical explanation for $`\sigma _\mathrm{\Delta }\frac{1}{2}\sigma _{Q\overline{Q}}`$ using the $`SU(3)_c`$ Casimir scaling. For the Coulomb coefficient, one gluon exchange brings the color factor $`1/2`$ in the perturbation theory as the ratio of the $`SU(3)_c`$ Casimir factor. The linear term is however nonperturbative and then there is no clear explanation to keep the $`SU(3)_c`$ Casimir scaling beyond the perturbation. In addition, the $`\chi ^2`$ of the $`\mathrm{\Delta }`$-type fit is much larger than the case of the $`Y`$-type fit, where the goodness of the fit is comparable with $`Q`$-$`\overline{Q}`$ potential case with $`r2`$. Although the present simulation may suffer from rather large finite $`a`$ effect, our result shows that the $`Y`$-type ansatz is more preferable with the lattice data.

Our conclusion at this stage is the following. $`Y`$-type ansatz is more preferable for the lattice result although $`\mathrm{\Delta }`$-type ansatz is not completely excluded with present systematic uncertainty. Adopting $`Y`$-type picture, we have found a universal feature of the string tension as $`\sigma _{3Q}\sigma _{Q\overline{Q}}`$. To confirm this result and investigate the 3$`Q`$ potential in more detail, simulations at larger $`\beta `$’s are in progress. Preliminary result on the $`16^3\times 32`$ lattice at $`\beta =6.0`$ is consistent with the result at $`\beta =5.7`$.

## 3 Finite temperature result

Now we turn on the temperature and observe how the potential changes with the thermal effect . In this paper, we focus on the static 3$`Q`$ potential in the deconfined phase. Above the critical temperature $`T_c`$, the quarks and the gluons are expected to be liberated and form a plasma state, in which the potential is expected to be screened. Then one expects that the screened potential takes the Yukawa potential like as the $`Q`$-$`\overline{Q}`$ potential , and genuine three-body potential (in $`Y`$-type picture) would disappear.

On the finite temperature lattice, the static potential is extracted from the Polyakov loop correlation,

$`P_3(\stackrel{}{r}_1,\stackrel{}{r}_2,\stackrel{}{r}_3)=P(\stackrel{}{r}_1)P(\stackrel{}{r}_2)P(\stackrel{}{r}_3)`$
$`=c\mathrm{exp}(V_{3Q}(\stackrel{}{r}_1,\stackrel{}{r}_2,\stackrel{}{r}_3)N_\tau ),`$ (6)
where the Polyakov loop $`P(\stackrel{}{x})`$ is defined as $`P(\stackrel{}{x})=\text{Tr}_{t=0}^{N_\tau 1}U_4(\stackrel{}{x},t)`$. The $`Q`$-$`\overline{Q}`$ potential is also measured in similar way.

To simulate the finite temperature system, we adopt the anisotropic lattice . An advantage of the anisotropic lattice for our study is that it enables us to work out several values of temperature without varying the lattice cutoff. As the gauge field action, we use the Symanzik improved action at the tree level .

The numerical calculation is done on the lattices of sizes $`16^2\times 24\times N_t`$ at $`\beta =4.56`$ and $`\gamma _G=3.45`$, in the quenched approximation. At $`N_t=96`$, which is realized as the zero temperature lattice, we determine the renormalized anisotropy as $`\xi =3.95(2)`$ from the Wilson loops in the spatial-spatial and temporal-spatial planes. The static $`Q`$-$`\overline{Q}`$ potential at $`T=0`$ gives the cutoffs $`a_\sigma ^1=1.61(1)`$ GeV and $`a_\tau ^1=6.36(5)`$ GeV. The critical temperature is estimated by the rapid change of the Polyakov loop and its susceptibility, and is found between $`N_t=24`$ and $`25`$, and slightly near to $`N_t=24`$. For the investigation of the static potentials in the deconfined phase, we generate 60 configurations at $`N_t=20`$, which roughly corresponds to $`T1.22T_c`$.

Figure 3 shows the $`Q`$-$`\overline{Q}`$ potential and the 3$`Q`$ potential in the deconfined phase. At the long distance, The Polyakov loop correlators approach to the multiples of the expectation value of single Polyakov loop: $`P_2PP^{}`$ and $`P_3P^3`$. These contributions are subtracted from the potentials. The $`Q`$-$`\overline{Q}`$ potential above $`T_c`$ is known to behave as the Yukawa-type potential . In present calculation, the Debye screening mass of the $`Q`$-$`\overline{Q}`$ potential strongly depends on the fit range and is estimated around 0.5 – 1 GeV. At this $`N_t`$, we have measured the 3$`Q`$ potential with only equilateral constellations. In Figure 3, we show the 3$`Q`$ potential divided by 3, and find good agreement with the $`Q`$-$`\overline{Q}`$ potential,

$$V_{3Q}(r_{12}=r_{23}=r_{31}=r)3V_{Q\overline{Q}}(r)(T>T_c).$$
(7)
Above $`T_c`$, since the quarks are deconfined, there is no reason to assume that the combination of quarks are in the color-singlet states. The question is whether the 3$`Q`$ potential can be represented as the sum of two body potential. Our present result shows that the $`V_{3Q}/3`$ behaves as almost same as the $`Q`$-$`\overline{Q}`$ potential at least for the equilateral 3$`Q`$ constellations. For more definite conclusion, we need precise results including other 3$`Q`$ constellations and a systematic comparison between $`V_{3Q}`$ and $`V_{Q\overline{Q}}`$ at various $`T`$.

The calculations have been done on NEC SX-4 at Research Center for Nuclear Physics, Osaka University and NEC HSP at INSAM (Institute for Nonlinear Science and Applied Mathematics), Hiroshima University.


format: text

text:
# Preliminary Spectral Analysis of the Type II Supernova 1999em

## 1 Introduction

SN 1999em was discovered by the Lick Observatory Supernova Search (IAUC 7294) on Oct 29.44 UT and confirmed by the BAO supernova group on Oct 29.7 UT. It was followed spectroscopically and photometrically by ground-based observatories and was the subject of a Supernova Intensive Studies (SINS) *Hubble Space Telescope* (*HST*) observation on Nov 5, 1999. SN 1999em is in the SABc galaxy NGC 1637, with a heliocentric recession velocity of 717 km s-1. The first spectrum taken on Oct 29, 1999 is extremely blue and displays strong Balmer features typical of Type II supernovae with a characteristic velocity of around 10,000 km s-1.

SN 1999em appears to have been a normal Type IIP, likely the result of a core-collapse in a star with a massive hydrogen envelope (Leonard et al., 2000). The early and weak radio emission (IAU Circulars 7318 and 7336) hints that the star did not experience significant mass loss shortly before the explosion, which is also thought to be the normal situation for SNe IIP (see Weiler et al., 2000, and references therein). The lack of circumstellar interaction makes SNe IIP “clean” cases to model, resulting in more reliable atmospheric modeling which in turn makes SNe IIP ideal for cosmological use. Because SN 1999em was observed earlier, and over a wider wavelength range than any other normal bright SN IIP, it is a perfect test bed for our understanding of SNe IIP. These supernovae are of particular interest as distance indicators.

## 2 Observations

The observed spectra we use in this paper were obtained on Nov 4, 1999 at the Fred L. Whipple Observatory (FLWO), and on Oct. 29, 1999 at CTIO. Leonard et al. (2000) find that on Oct. 29 the *V* light curve is rising. The point at Nov. 4-5 is on a broad maximium prior to the plateau. The light curve settles down to a plateau phase on about Nov. 20.

## 3 SYNOW Models

We have used the fast, parameterized supernova synthetic spectra code SYNOW to make an initial investigation of line identifications and expansion velocities. The code is discussed in detail by Fisher (2000) and recent applications include Fisher et al. (1999), Millard et al. (1999), and Hatano et al. (1999a). SYNOW spectra consist of resonant–scattering profiles superimposed on a blackbody continuum. For the two synthetic spectra presented here, the radial dependence of all line optical depths is a power law of index 8. The observed spectra have been corrected for a redshift of 717 km s-1 and for reddening of $`E(BV)=0.05`$ mag (see below).

In Figure 1 the CTIO observed spectrum of October 29 is compared with a synthetic spectrum for which the continuum blackbody temperature $`T_{bb}`$ and the excitation temperature $`T_{exc}`$ are 13,000 K, the velocity at the photosphere $`V_{phot}`$ is 11,000 km s-1, and only lines of hydrogen, He I, and N II are considered. The relative strengths of the lines of an ion are fixed by local thermodynamic equilibrium (LTE) at the excitation temperature, but the absolute strengths of the lines of each ion are controlled by a free parameter. The strength of the hydrogen lines has been chosen to make H$`\beta `$ and H$`\gamma `$ about the same as in the observed spectrum; the relative optical depths of the hydrogen lines are of course fixed by atomic physics so the poor fit to H$`\alpha `$ is due to a deficiency of the resonant–scattering source function for this transition. The only noticeable line of He I in the synthetic spectrum is that of $`\lambda `$5876. Lines of N II have been introduced so that $`\lambda 4623`$ and $`\lambda 5679`$ can account for the absorptions near 4500 and 5500 Å, in which case N II $`\lambda 5029`$ also affects the synthetic spectrum slightly near 5000 Å. (In SYNOW when a particular species in introduced, all the lines of that species are included with their relative strengths set by assuming LTE populations, but we focus here on the strong optical features which are relevant to the observed spectra.) The hydrogen and He I identifications are definite, and (in anything remotely like LTE) we can offer no alternative to the N II lines. However, see § 4.1.1 regarding the possibility that some of the “N II” lines may actually be due to hydrogen.

In Figure 2 the *HST* observed spectrum of November 5 is compared with a synthetic spectrum that has $`T_{bb}=T_{exc}=9500`$ K, $`v_{phot}=8000`$ km s-1, and lines of hydrogen, Ca II, Mg II, Fe II, and Ni II. The resonance–scattering approximation gives good fits to the features produced by H$`\beta `$, H$`\gamma `$, and H$`\delta `$. Ca II contributes only the blend due to $`\lambda `$3945 (the H and K lines) and Mg II contributes little other than the blend due to $`\lambda 2798`$ (the h and k lines). All other features in the synthetic spectrum are produced by Fe II and, to a lesser extent, by Ni II. The hydrogen, Ca II, and Mg II identifications are definite. The ultraviolet spectrum evidently is mainly a blend of singly–ionized iron–peak lines; Fe II definitely is present and Ni II probably contributes significantly. Other iron–peak ions that have not been introduced here may also affect the observed spectrum.

## 4 PHOENIX Models

### 4.1 Oct 29

The Galactic extinction to NGC 1637 corresponds to a color excess of $`E(BV)=0.03`$ mag, but the observed spectrum of SN 1999em shows a distinct NaD interstellar absorption line that may indicate additional reddening in the parent galaxy. We have calculated a grid of detailed fully line-blanketed PHOENIX models in order to determine the “temperature” of the observed spectrum. With most supernova spectra we can trade off higher temperatures for larger reddening within prescribed limits (such that the atmosphere does not become so hot or cool that strong unobserved lines would predominate). Here, with this very early, very blue spectrum, we found that we were not able to freely exchange temperature for reddening, so that the reddening was well determined to be $`E(BV)0.05`$ and the “temperature” was about 11,000 K. Figure 3 displays a fit with $`T_{\text{model}}=11,000`$ K and solar compositions; the observed spectrum has been dereddened with $`E(BV)=0.05`$ mag. The abundances were taken to be solar throughout the model atmosphere and the density was assumed to follow a power-law with $`\rho r^n`$ and $`n=7`$, somewhat shallower than in the SYNOW fits, though we don’t regard the difference as being too significant, since the true density structure is unlikely to follow an exact power-law. We also present models with $`n=9`$ (see Figure 8).

The fit is rather good, and the $`T_{\text{model}}`$ cannot be significantly reduced since the Ca II H+K line becomes very strong in the synthetic spectrum with $`T_{\text{model}}=10,000`$ K, but it is weak in the observed spectrum. Also, the observed feature due to He I $`\lambda 5876`$ is not reproduced in the synthetic spectrum. Hotter models are not ruled out but at $`T_{\text{model}}=12,000`$ K the Ca II H+K feature becomes weaker and there is a clear flux deficiency in the red. Normally one expects that one can trade off reddening for temperature in the red part of the spectrum, and this is indeed the case within a certain range of reddening. We assume a standard (Cardelli et al., 1989) reddening law with $`R=3.1`$. However, since the spectrum depends strongly on the strength of the Ca II H+K feature, we find that $`T_{\text{model}}=11,000\pm 500`$ K; the $`T_{\text{model}}=12,000`$ K model requires an extinction $`E(BV)0.15`$ to reproduce the red flux, but it then does a poor job in the blue. Thus, we find that $`E(BV)0.05`$ and $`E(BV)<0.15`$. We show a model with $`E(BV)=0.10`$ below.

We also found that we could not reproduce the He I $`\lambda 5876`$ feature by enhancing the amount of gamma-ray deposition. Therefore, guided by the indications of enhanced N in the SYNOW fits, we calculated a series of spectra where we enhanced He by a factor of 2.5 (reducing H appropriately), and we enhanced N by a factor of 10 (reducing C+O). This fit is shown in Figure 4, where now the He I $`\lambda 5876`$ fits well, but $`T_{\text{model}}=11,500`$ K in order to keep the Ca H+K feature from becoming too strong in the synthetic spectrum. Figure 5 with $`T_{\text{model}}=12,000`$ K, the same compositions and $`E(BV)=0.10`$ shows that the effects of the Ca H+K feature can be reduced by increasing the model temperature and increasing the reddening modestly. Therefore we conclude that there is strong evidence for enhanced helium in SN 1999em \[a similar result was found for SN 1987A (Arnett et al., 1989; Eastman & Kirshner, 1989; Lundqvist & Fransson, 1996; Blinnikov et al., 2000)\]. We note at these early times that the material is hot enough for the optical helium lines to be excited even in LTE (see Hatano et al., 1999b), so while gamma-ray deposition was needed in e.g. SN 1993J (see Baron et al., 1995a, and references therein), it is not required at these early times (nor is it precluded at later times, although we did not find it necessary for the spectra a week later).

#### 4.1.1 Complicated Balmer and He I Lines

We attempted to confirm the evidence in the SYNOW fits that nitrogen is also enhanced and that N II contributes to the observed features at $`\lambda 5679,\lambda 5029,\text{and}\lambda 4623`$. To do this we reduced the abundance of nitrogen from its solar value by a factor of ten. This synthetic spectrum showed that the $`\lambda 5679`$ and $`\lambda 5029`$ features do decrease, but the feature that is attributed to $`\lambda 4623`$ in the SYNOW fits is still present (this is not surprising, since the feature also appears in Figure 3 with solar compositions). Guided by the LTE line spectra presented by Hatano et al. (1999b), we reduced carbon by a factor of ten from its solar value and the 4600 Å feature remained. Finally, reducing both carbon and oxygen by a factor of ten from their solar values (leaving N fixed) still produces a feature with if anything a better shape than in previous calculations. Thus, with the “usual suspects” the source of the 4600 Å feature remains unidentified, and the evidence for enhanced nitrogen is weakened, but it does a good job on the $`\lambda 5679,\lambda 5029,`$ features.

The lack of a strong candidate for the 4600 Å feature leads us to consider more exotic line identifications. Careful examination of Figure 3 reveals that even with solar compositions features appear in the synthetic spectrum at the approximately correct wavelength positions. This led us to test the hypothesis that the features are actually produced by “complicated P-Cygni” profiles, that is the usual wide P-Cygni profile with the peak centered near zero velocity and a second P-Cygni profile with a second absorption minimum around 20,000 km s-1. Figure 6 displays the synthetic calculation with the metalicity reduced to $`Z=Z_{}/100`$. One can be confident that all the lines in this synthetic spectrum are produced solely by hydrogen and helium and it is clear that each line is associated with high-velocity absorption, which matches the observed lines very well. (It may be a true P-Cygni feature, with both a peak and a dip, but it is probably premature to decide that at this point.) This is due to the shallow density gradient and high model temperature which keeps the actual electron temperature above 8000 K all the way to the highest velocity in the model, 37,500 km s-1. Complicated NLTE effects vary the Balmer level populations in the mid-velocity range, producing the double P-Cygni lines, i.e. there exist two line forming regions for hydrogen and helium. This is the first time that we have encountered this result and it shows that particularly in differentially expanding flows, one cannot be certain of line identifications until careful synthetic spectral models have been calculated. Thus, we therefore believe that the features that we have heretofore attributed to N II (and which would require enhanced N) are in actuality due to H I and He I, which form secondary features at high velocity. Since this appears in our modeling without any special parameters (i.e. Figure 3), we believe this is a true NLTE effect produced simply by the shallow density gradient.

### 4.2 Nov 4–5

Guided by our results in § 4.1 we have computed full NLTE models for the epoch observed with *HST* using both solar and enhanced He+N compositions (since we have not ruled out enhanced N, we will continue to use it for expediency). Figure 7 shows our best fit for solar compositions and $`n=7`$ with $`T_{\text{model}}=8100`$ K. Figure 8 displays our best fit for the enhanced He+N compositions, where we have also found the fit is improved by setting $`T_{\text{model}}=9000`$ and $`n=9`$. Also displayed is a ground-based optical spectrum obtained on Nov. 4, 1999 at the FLWO, which lends support for enhanced He. All the models show unobserved strong features in the range 2000-3500 Å, which could be due to iron-peak elements. Thus, it is possible that the metallicity of SN 1999em was somewhat less than solar, but we will study this question in future work.

## 5 Conclusions

We have shown that both direct synthetic spectral fits and detailed NLTE models do a good job of reproducing the observed optical+UV spectra of SN 1999em. We have shown that detailed NLTE spectral modeling of very early spectra of Type II supernova can provide an independent estimate to the total reddening as well as abundance information. Specifically we find that $`E(BV)0.050.10`$ and $`E(BV)<0.15`$ mag, and we find strong evidence for helium enhanced by at least a factor of 2 over the solar value. We obtain the striking result that observed features are likely due to “complicated P-Cygni” lines of hydrogen and helium, and thus one needs to be especially careful with line identifications in hot differentially expanding flows. That the observed complicated profiles fall right out of the PHOENIX calculations is remarkable and lends strong support to the reliablity of the modeling.

The fact that planned supernova search programs will preferentially find supernovae very early, combined with our results that the extinction can be determined by modeling the very early spectra, lends promise to the use of Type II supernovae as distance indicators, through the use of phenomenological applications like the “Expanding Photosphere Method” (Baade, 1926; Kirshner & Kwan, 1974; Branch et al., 1981; Schmidt et al., 1992; Eastman et al., 1996) or through sophisticated modeling of the “Spectral-fitting Expanding Atmosphere Method” (Baron et al., 1993, 1994, 1995a, 1995b, 1996).

We thank Doug Leonard and Weidong Li for sharing their results in advance of publication and for helpful discussions. PHH was supported in part by the Pôle Scientifique de Modélisation Numérique at ENS-Lyon. This work was supported in part by NSF grants AST-9731450, AST-9417102, and AST-9417213; by NASA grant NAG5-3505, and an IBM SUR grant to the University of Oklahoma; by NSF grant AST-9720704, NASA ATP grant NAG 5-8425 and LTSA grant NAG 5-3619 to the University of Georgia; and by NASA GO-8243 and GO-8648 to the SINS group from the Space Telescope Science Institute, which is operated by AURA, Inc. under NASA contract NAS 5–26555. Some of the calculations presented in this paper were performed at the San Diego Supercomputer Center (SDSC), supported by the NSF, and at the National Energy Research Supercomputer Center (NERSC), supported by the U.S. DOE. We thank both these institutions for a generous allocation of computer time. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration.


text:
# Absolute frequency measurement of the In+ clock transition with a mode-locked laser

## Abstract

The absolute frequency of the In+ $`5s^{\mathrm{2\hspace{0.17em}1}}S_0`$ \- $`5s5p^{\mathrm{\hspace{0.17em}3}}P_0`$ clock transition at 237 nm was measured with an accuracy of 1.8 parts in $`10^{13}`$. Using a phase-coherent frequency chain, we compared the $`{}_{}{}^{1}S_{0}^{}`$ \- $`{}_{}{}^{3}P_{0}^{}`$ transition with a methane-stabilized He-Ne laser at 3.39 $`\mu `$m which was calibrated against an atomic cesium fountain clock. A frequency gap of 37 THz at the fourth harmonic of the He-Ne standard was bridged by a frequency comb generated by a mode-locked femtosecond laser. The frequency of the In+ clock transition was found to be $`\mathrm{1\hspace{0.33em}267\hspace{0.33em}402\hspace{0.33em}452\hspace{0.33em}899.92}(0.23)`$ kHz, the accuracy being limited by the uncertainty of the He-Ne laser reference. This represents an improvement in accuracy of more than 2 orders of magnitude on previous measurements of the line and now stands as the most accurate measurement of an optical transition in a single ion.

Improving the accuracy and stability of frequency standards using narrow transitions in atomic systems at optical or higher frequencies as a reference for laser oscillators has been a long-standing objective in metrology . Impressive progress has been achieved in short-term stability of oscillators operating in these spectral regions  as well as in long-term stabilization by directly locking lasers onto narrow atomic transitions . Major advances have recently been reported in the ability to determine absolute optical frequencies using harmonic frequency chains  or frequency chains based on precise measurement of large frequency intervals .

Highly accurate and stable optical frequency standards enable precise measurement of fundamental constants  or investigation of their possible variation in time . Better frequency standards permit more accurate determination of atomic transitions in spectroscopy and offer the possibility of stringent testing of QED or general relativity. They are also needed in applications such as navigation or very long baseline interferometry .

We study the $`5s^{\mathrm{2\hspace{0.17em}1}}S_0`$ \- $`5s5p^{\mathrm{\hspace{0.17em}3}}P_0`$ transition at 237 nm in a single trapped In+ ion with a view to applying it as an optical frequency standard . This transition with a natural linewidth of only 0.8 Hz  has a line quality factor of $`Q=\nu /\mathrm{\Delta }\nu =1.6\times 10^{15}`$. The two states participating in the hyperfine-induced $`J=0J=0`$ transition couple only weakly to perturbing external electric or magnetic fields. In particular, they are insensitive to electric quadrupole shifts caused by the field gradient of the trap. The particle can be kept in the trap almost free of collisions for an infinite time and brought to near-rest by laser cooling. Thus, transit-time broadening and first-order Doppler effects can be eliminated and the second-order Doppler effect can be reduced to negligible values . With these characteristics taken into account, a residual uncertainty of $`10^{18}`$ for a frequency standard based on a single stored In+ ion is expected .

In a recent publication , we reported the measurement of the absolute frequency of the In+ clock transition using two optical reference frequencies, a methane-stabilized He-Ne laser at 3.39 $`\mu `$m and a Nd:YAG laser at 1064 nm whose second harmonic was locked to a hyperfine component in molecular iodine. The measurement inaccuracy of 3.3 parts in 1011 was limited by the degree of uncertainty to which the iodine reference was known. Here we report on a new precise measurement of the absolute frequency of the $`{}_{}{}^{1}S_{0}^{}`$ \- $`{}_{}{}^{3}P_{0}^{}`$ transition using a phase-coherent frequency chain which links the 237 nm radiation (1 267 THz) of the In+ clock transition to the He-Ne laser at 3.39 $`\mu `$m (88 THz) alone. In this case the accuracy of the measurement is only limited by the much smaller uncertainty of the CH4/He-Ne reference.

The transportable He-Ne standard was built at the Institute of Laser Physics in Novosibirsk, Russia . It was previously (1996) calibrated for a measurement of the hydrogen 1S - 2S absolute frequency . Here we use the result of a more recent calibration  that was obtained together with a new cesium-based measurement of the hydrogen 1S - 2S interval , with $`f_{HeNe}=\mathrm{88\hspace{0.33em}376\hspace{0.33em}182\hspace{0.33em}599\hspace{0.33em}976}(10)`$ Hz. This value deviates from the previous one by 39 Hz (1.6 combined standard deviations), most likely because the operating parameters were not exactly maintained for several years. The new calibration was done 4 months before the measurement of the In+ clock transition. Unlike in the previous calibration, the laser was not moved between its calibration and the measurement.

The experimental setup of the frequency chain is illustrated in Fig. 1. A NaCl:OH- color center laser at 1.70 $`\mu `$m is phase-locked to the second harmonic of the methane-stabilized He-Ne laser at 3.39 $`\mu `$m. A 848 nm laser diode is then locked to the second harmonic of the color center laser. This is accomplished by first stabilizing it to a selected mode of the frequency comb of a Kerr-lens mode-locked Ti:sapphire femtosecond laser (coherent model Mira 900), frequency-broadened by a standard single-mode silica fiber (Newport FS-F), and then controlling the position of the comb in frequency space . A radio-frequency divider, dividing by 128, helps to overcome the limited servo bandwidth that controls the comb position. Recently, we demonstrated that this frequency comb meets the exceptional requirements of high-precision spectroscopy . The 76 MHz pulse repetition rate, which sets the mode separation, is phase-locked  to a commercial cesium clock (Hewlett-Packard model 5071 A) which is constantly calibrated using the time disseminated by the Global Positioning System. Each mode of the comb is therefore known with the same fractional precision as the He-Ne standard, i.e. within 1.1 parts in $`10^{13}`$ . A diode laser at 946 nm is phase-locked to another selected mode of the comb, positioned $`n=\mathrm{482\hspace{0.33em}285}`$ modes or 37 THz to lower frequencies from the initial mode at 848 nm. The mode number $`n`$ is determined from our previous measurement . The beat frequency with the 946 nm Nd:YAG laser, whose 4th harmonic excites the In+ clock transition, is counted with a commercial radio-frequency counter.

For the absolute frequency measurement, the in-lock beat signals of the chain are continuously monitored with additional counters for possible lost cycles . Points that are off by more than 0.5 Hz from the phase-locked beat signal are discarded. With the frequency chain in lock, the unknown 115In+ $`{}_{}{}^{1}S_{0}^{}`$ \- $`{}_{}{}^{3}P_{0}^{}`$ clock transition frequency $`f_{In^+}`$ at 237 nm is related to the known frequency of the He-Ne standard $`f_{HeNe}`$ through

$$f_{In^+}=16f_{HeNe}4(f_B+nf_r)f_{LO}.$$
(1)
Here $`f_B`$ is the frequency of the beat signal at 946 nm detected by a photodiode and recorded by a radio-frequency counter in a 1 sec interval, $`f_r`$ denotes the repetition rate of the mode-locked femtosecond laser, $`n`$ is the number of modes separating the two selected modes of the comb, and $`f_{LO}`$ = 1 632 MHz contains all contributions from the local oscillator frequencies employed for the phase-locks.

Within several days of measurement we recorded 6 214 excitations of the $`{}_{}{}^{3}P_{0}^{}`$ level by scanning the 946 nm Nd:YAG clock laser over the In+ $`{}_{}{}^{1}S_{0}^{}`$ \- $`{}_{}{}^{3}P_{0}^{}`$ resonance. Excitation of the $`{}_{}{}^{3}P_{0}^{}`$ state is detected by optical-optical double resonance : after applying the clock laser for 15 ms the laser exciting the $`{}_{}{}^{1}S_{0}^{}`$ \- $`{}_{}{}^{3}P_{1}^{}`$ cooling transition is turned on for 40 ms to probe the population of the ground state. If no fluorescence photons are counted on the $`{}_{}{}^{1}S_{0}^{}`$ \- $`{}_{}{}^{3}P_{1}^{}`$ transition, the cooling laser is kept on for up to ten further 40 ms intervals to wait for the decay of the metastable state and an excitation event is recorded. The frequency of the clock laser is typically changed in steps of 80 Hz and 16 excitation attempts are made at each frequency. Scanning of the Nd:YAG master laser is synchronized with detection of the beat signal $`f_B`$ at 946 nm and the counters used for detection of lost cycles. All counters as well as the AOM driving synthesizers are referenced to a local cesium atomic clock.

Figure 2a shows the $`{}_{}{}^{3}P_{0}^{}`$ excitation probability as a function of the beat frequency $`f_B`$ for a typical measurement session, collecting 674 quantum jumps to the $`{}_{}{}^{3}P_{0}^{}`$ state. During the measurement session 21 excitation spectra were recorded where the investigated scanning range varied slightly from spectrum to spectrum. The width of the distribution is due to short-term and mid-term frequency instabilities of the Nd:YAG laser, the methane standard, and the frequency chain. A weighted fit to a Gaussian function is used to determine the line center. By averaging the beat frequency at line center for the eleven measurement sessions performed we obtain $`f_B=\mathrm{49\hspace{0.33em}174\hspace{0.33em}925}(42)`$ Hz (see Fig. 2b). From this value we determine the absolute frequency $`f_{In^+}`$ of the 115In+ $`{}_{}{}^{1}S_{0}^{}`$ \- $`{}_{}{}^{3}P_{0}^{}`$ clock transition according to Eq. (1):

$$f_{In^+}=\mathrm{1\hspace{0.33em}267\hspace{0.33em}402\hspace{0.33em}452\hspace{0.33em}899.\hspace{0.33em}92}(0.23)\text{kHz}.$$

The uncertainty of 1.8 parts in 1013 is derived from a quadratic addition of the uncertainty of the measurement (4 $`\mathrm{\Delta }f_B`$) and the uncertainty in the He-Ne standard calibration (16 $`\mathrm{\Delta }f_{HeNe}`$). Both contributions are due to the limited reproducibility of the He-Ne standard.

Systematic frequency shifts of the clock transition are negligible at the present level of accuracy. The magnetic field dependence is -636 $`\pm `$ 27 Hz/G for the $`F=9/2;m_F=9/2>F=9/2;m_F=7/2>`$ component that we excite , where we adjust the magnetic field to zero with a precision of a few times 10 mG. Other systematic frequency shifts such as the quadratic Stark or second-order Doppler shift are orders of magnitude smaller than the Zeeman shift at the temperatures to which the ion is cooled in our trap (T $``$ 150 $`\mu `$K) .

The new result is well within the error bars of our previous measurement, where we obtained a value of $`\mathrm{1\hspace{0.33em}267\hspace{0.33em}402\hspace{0.33em}452\hspace{0.33em}914}(41)`$ kHz for the clock transition . Compared with that value, the new measurement represents an improvement in accuracy of more than 2 orders of magnitude. It now stands as the most accurate measurement of an optical transition frequency in a single ion. For the future it is planned to use a femtosecond laser frequency comb to compare the indium clock transition directly with a cesium clock. In this case the cesium reference could be provided by a transportable fountain clock, as was done in a recent measurement of the 1S - 2S transition in hydrogen . The stability and reproducibility of the $`{}_{}{}^{1}S_{0}^{}`$ \- $`{}_{}{}^{3}P_{0}^{}`$ clock transition makes indium an attractive candidate for such ultra-precise frequency measurement.

A further interesting feature of our frequency chain is the fact that it can be used to compare two narrow optical transition frequencies with each other, e.g. the clock transition in indium with the hydrogen 1S - 2S transition. In particular, this may allow investigation of possible variations of fundamental constants in time as recently proposed in .

The collaboration with ILP, Novosibirsk, was supported by the Volkswagen Foundation within Project No. I/72 607.


text:
# Role of entropy barriers for diffusion in the periodic potential

## I Introduction

A variety phenomena in physics and other fields can be modeled as Brownian motion in an external periodic potential . One of particular example, the surface diffusion of atoms or small clusters, is of great fundamental and technological interest . During crystal growth the deposited atoms diffuse over the surface until they become incorporated in the lattice. On the semiconductor Si(100) or Ge(100) surface, most of the deposited Si or Ge atoms combine to form dimers, and the diffusion of such dimers has recently been studied experimentally by a scanning tunneling microscope . Moreover, atoms adsorbed on metal surfaces in some cases form closely packed islands which diffuse as a whole .

Theoretically the problem of diffusion can be described by a Langevin equation for the atom(s) or, equivalently, a Fokker-Planck-Kramers equation for the distribution function in the phase space . In the trivial case of Brownian motion without the external potential, the diffusion coefficient is equal to $`D=D_fk_BT/m\eta `$, where $`k_B`$ is the Boltzmann constant, $`T`$ is the temperature, $`m`$ is the particle mass, and $`\eta `$ is the viscous damping coefficient which models the energy exchange with the substrate (thermostat). For a single atom in the one-dimensional (1D) sinusoidal external (substrate) potential the value of diffusion coefficient is well known as summarized in the Risken monograph . Exact results exist for the overdamped (Smoluchowski) case, $`\eta \mathrm{}`$, when

$$D=D_f\mathrm{I}_0^2\left(h\right),$$
(1)
where $`h\epsilon /2k_BT`$, $`\epsilon `$ is the (total) height of the substrate potential and $`\mathrm{I}_0(h)`$ is the modified Bessel function, and in the underdamped limit, $`\eta 0`$, when

$$D=D_fG(h),$$
(2)
where $`G(h)=\left(h/2\pi \right)^{1/2}e^h\mathrm{I}_0^1\left(h\right)J(h)`$, $`J(h)=_0^1𝑑uu^{3/2}e^{2h/u}𝐄^1(u)`$, and $`𝐄(u)`$ is the complete elliptic integral of second kind. At low temperatures, $`k_BT\epsilon `$, both expressions (1) and (2) take the Arrhenius form, $`D=\stackrel{~}{D}A`$ with $`A=\mathrm{exp}\left(\epsilon /k_BT\right)`$ and $`\stackrel{~}{D}\omega _0^2a^2/2\pi \eta `$ in the high-friction case \[here $`a`$ is the period of the substrate potential and $`\omega _0=(2\pi /a)(\epsilon /2m)^{1/2}`$ is the frequency of oscillation at its minimum\], and $`\stackrel{~}{D}\pi D_f/2`$ in the low-friction limit. In a general case the diffusion coefficient can be found numerically with practically any desired accuracy by the matrix continued-fraction-expansion method .

At low temperatures, $`k_BT\epsilon `$, when the diffusion proceeds by uncorrelated thermally activated jumps over the barrier from one minimum of the external potential to another, the diffusion coefficient may be presented as $`D=RA\lambda ^2`$, where $`RA`$ is the rate of escape from a minimum of the external potential and $`\lambda ^2`$ is the mean-square jump length. For a moderate or large damping, $`\eta \omega _0`$, when only jumps for one period $`a`$ of the external potential are possible, one should take $`\lambda =a`$ and $`R=R_{\mathrm{TST}}B(\eta )`$, where $`R_{\mathrm{TST}}=\omega _0/2\pi `$ is the escape rate given by the transition state theory (TST) , and the factor $`B(\eta )=(z^2+1)^{1/2}z`$ with $`z=\eta /2\omega _s`$ provides an interpolation between the TST and overdamped limits as was found by Kramers  \[here $`\omega _s`$ is the “saddle” frequency at the saddle point $`x=x_s`$, near which the external potential has the form $`V(x)\epsilon \frac{1}{2}m\omega _s^2(xx_s)^2`$; for the sinusoidal potential $`\omega _s=\omega _0`$\]. The underdamped limit, $`\eta \omega _0`$, is qualitatively different: in this case $`R2\eta \epsilon /\pi k_BT\eta `$ as was found firstly by Kramers , but the average jump length diverges as $`\lambda \eta ^1`$, thus this again leads to the dependence $`R\lambda ^2\eta ^1`$ similarly to the overdamped case. The occurrence of long jumps, $`\lambda >a`$, has been observed in a number of experiments on surface diffusion . The interval from low to moderate friction is covered by the Mel’nikov-Meshkov formula

$$R_{\mathrm{MM}}R_{\mathrm{TST}}\mathrm{exp}\left(\frac{1}{\pi }_0^{\mathrm{}}𝑑u\frac{\mathrm{ln}\left[1e^{\mathrm{\Delta }\left(u^2+\frac{1}{4}\right)}\right]}{u^2+\frac{1}{4}}\right),$$
(3)
where $`\mathrm{\Delta }=8h\eta /\omega _0`$. Thus, the whole interval of frictions may be described by the interpolation formula $`RR_{\mathrm{MM}}B(\eta )`$, which was checked numerically in . Combining this expression for $`R`$ with the numerically calculated values of $`D`$, one can find the distribution of jump lengths . Note that the widely used TST expression $`DR_{\mathrm{TST}}Aa^2`$, where the diffusion coefficient does not depend on the damping coefficient, operates in fact for a narrow interval of frictions close to the point $`\eta \omega _0`$ only (which, fortunately, often corresponds to experimental situations).

Although the described above results for one-dimensional diffusion are very important and often lead to reasonable estimations for experimentally measured diffusion coefficients, in real systems the motion always takes place in a $`N>1`$ configurational space. Indeed, even for diffusion of a single atom adsorbed on a crystal surface $`N=2`$ at least. Besides, the diffusing object may have internal degrees of freedom. Multi-dimensional effects modify both the escape rate $`R`$ and the jump length $`\lambda `$. The escape rate can be presented as $`R=R_{1\mathrm{D}}F`$, where the coefficient $`F`$ is known as the “entropy factor” . The value of $`F`$ can be found with the help of transition state theory  which yields $`F\left(\mathrm{\Pi }_i\omega _{0,i}\right)/\left(\mathrm{\Pi }_i\omega _{s,i}\right)`$, where $`\omega _{0,i}`$ are the frequencies at the minimum and $`\omega _{s,i}`$ are the “saddle” frequencies for all degrees of freedom $`i`$ except the given one along the diffusion path. In this approach $`F`$ can be interpreted as $`F=\mathrm{exp}(\mathrm{\Delta }S/k_B)`$, where $`\mathrm{\Delta }S`$ is the difference in entropy of the saddle and minimum-energy configurations. The entropy factor is often used to explain the “compensation effect” , when at experiment one observes that a decrease of the activation energy (calculated as a slope of the Arrhenius plot of $`\mathrm{ln}D`$ versus $`T^1`$) is compensated by decreasing of the prefactor. As for the jump length, while for $`\eta \omega _0`$ it still is given by $`\lambda =a`$, in the underdamped limit it is modified qualitatively comparing with the one-dimensional case. In the multi-dimensional space, the path connecting adjoining minima of the external potential may not coincide with the direction of easy crossing at the saddle point. Therefore, the probability of deactivation during long jumps is enhanced, leading to the reduction of jump length, $`\lambda <\lambda _{1\mathrm{D}}`$ . In particular, for the 2D-periodic substrate potential with the square symmetry it was found numerically  that $`D\eta ^{0.5}`$ which gives $`\lambda \eta ^{0.75}`$.

The multi-dimensional effects are also important in diffusion of molecules or small clusters: even for diffusion in the 1D periodic potential (e.g., along “channels” on furrowed or stepped surfaces) one has for the dimer diffusion $`N=2`$ at least. Diffusion of the dimer was studied numerically by Vollmer  with the help of matrix continued-fraction-expansion technique. The adiabatically slow motion of a linear molecule in the 1D sinusoidal potential was analyzed in , where the adiabatic trajectory was found for a general case. This allowed to find the activation barriers and the minimum-energy and saddle-state frequencies and then to estimate the diffusion coefficient.

The aim of the present paper is to study the multi-dimensional effects in diffusion processes. We consider two typical examples: motion of a single atom in a “channel” which is periodic in one dimension and parabolic in others, and diffusion of a dimer (two-atomic molecule) in the 1D sinusoidal potential. We find an analytical solution for the overdamped case and analyze numerically the dependence of diffusion coefficient on the damping constant $`\eta `$. The numerical results were obtained with the Verlet algorithm by calculating the trajectory $`x(t)`$ and then splitting it into $`N_{\mathrm{tr}}`$ pieces, each of the time duration $`\tau `$. The diffusion coefficient was then calculated as $`D=\mathrm{\Delta }x^2/2\tau `$.

The paper is organized as follows. In Sec. II we obtain the analytical expression for the diffusion coefficient in the overdamped limit. In Sec. III A we analyze the case of pure entropic barriers. In Sec. III B the activated diffusion of a single atom is studied. In Sec. III C the diffusion of a dimer is described. Finally, Sec. IV concludes the paper.

## II Overdamped limit

Consider a particle moving in the $`N`$-dimensional external potential $`V_N(x,y_1,\mathrm{},y_{N1})`$ which is periodic in the $`x`$ direction,

$$V_N(x+a,\mathrm{})=V_N(x,\mathrm{}),$$
(4)
and grows unboundently in the other $`N1`$ dimensions,

$$V_N(x,y_1,\mathrm{},y_{N1})\frac{1}{2}\omega _i^2(x)y_i^2\text{if}y_i\pm \mathrm{},$$
(5)
where $`\omega _i(x)>0`$ for all $`x`$.

With presence of a viscous friction, the particle motion should be diffusive at long-time scale. The diffusion coefficient $`D`$ can be found with the Einstein relation $`D=T\mu `$, where the mobility $`\mu `$ describes the proportionality between the linear current $`j`$ and the infinitesimal external dc force $`f`$ which causes this current, $`j=\mu f`$. Therefore, we have to consider the particle motion in the external potential

$$V_f(x,y_1,\mathrm{},y_{N1})=V_Nfx,$$
(6)
and then take the limit $`f0`$.

In the overdamped case, when the friction coefficient $`\eta `$ is much larger than the characteristic system frequencies, the motion of the particle is described by the Smoluchowski equation

$$\frac{W}{t}+\stackrel{}{}\stackrel{}{J}=0,\stackrel{}{J}=\eta ^1(W\stackrel{}{}V_f+T\stackrel{}{}W),$$
(7)
where $`W(x,y_1,\mathrm{},y_{N1};t)`$ is the distribution function, $`\stackrel{}{J}(x,y_1,\mathrm{},y_{N1};t)`$ is the density of particle’s current, and the particle mass and Boltzmann constant are put to unity, $`m=1`$ and $`k_B=1`$.

For a steady state, Eq. (7) takes the form

$$T\frac{W}{x}+W\frac{V_f}{x}=\eta J_x$$
(8)
for the $`x`$ component, and a similar form for other degrees of freedom. The density $`\stackrel{}{J}`$ of the current should satisfy the equation

$$\frac{J_x}{x}+\underset{i=1}{\overset{N}{}}\frac{J_{y_i}}{y_i}=0.$$
(9)
To reduce notations, below we consider the case of $`N=2`$ only; generalization to the $`N>2`$ case is trivial. Let us introduce the one-dimensional density and current as

$$\rho (x)=_{\mathrm{}}^+\mathrm{}𝑑yW(x,y),$$
(10)
$$j(x)=_{\mathrm{}}^+\mathrm{}𝑑yJ_x(x,y).$$
(11)
Owing to the condition (5), the current $`j(x)`$ does not depend on $`x`$,

$$\frac{dj(x)}{dx}=J_y(x,+\mathrm{})+J_y(x,\mathrm{})=0,$$
(12)
where we have used Eq. (9). Thus, integrating both parts of Eq. (8) over $`y`$, we obtain the one-dimensional equation

$$T\frac{d\rho (x)}{dx}+\rho (x)\frac{dV_F(x)}{dx}=\eta j,$$
(13)
where we introduced the potential $`V_F(x)`$ defined by the equation

$$\frac{dV_F(x)}{dx}=\left[\rho (x)\right]^1_{\mathrm{}}^+\mathrm{}𝑑yW(x,y)\frac{V_f(x,y)}{x}.$$
(14)
Now, if $`V_F(x)`$ may be presented in the form

$$V_F(x)=V_N(x;f)fx,$$
(15)
where $`V_N(x;f)`$ is a periodic function on $`x`$, Eq. (13) takes the form studied in , and the diffusion coefficient can be calculated as

$$D=D_f(I_+I_{})^1,I_\pm (T)=(2\pi )^1_0^{2\pi }𝑑xe^{\pm V_{\mathrm{eff}}(x)/T},$$
(16)
where $`D_f=T/\eta `$ and $`V_{\mathrm{eff}}(x)=lim_{f0}V_N(x;f)`$. Thus, the diffusion coefficient $`D`$ is determined by the one-dimensional function $`V_N(x;0)`$. In the limit $`f0`$ we may substitute the equilibrium distribution function $`W=W_{\mathrm{eq}}\mathrm{exp}(V_N/T)`$ into Eq. (14), thus obtaining

$$\frac{dV_{\mathrm{eff}}(x)}{dx}=\frac{_{\mathrm{}}^+\mathrm{}𝑑ye^{V_N(x,y)/T}V_N(x,y)/x}{_{\mathrm{}}^+\mathrm{}𝑑ye^{V_N(x,y)/T}}.$$
(17)
Emphasize that this is the key approximation which is rigorous in the overdamped limit only. For the underdamped case, $`\eta 0`$, a similar multiplicative separation in the Fokker-Planck-Kramers equation, $`W(x,y,v_x,v_y,f)W(x,v_x,f)W_{\mathrm{eq}}(y,v_y)`$, does not work even in the $`f0`$ limit.

Let $`V_N(x,y)`$ has the form

$$V_N(x,y)=V(x)+U(y)+v(x,y),$$
(18)
where the function $`v(x,y)`$ describes the coupling between the two degrees of freedom. Then the effective potential $`V_{\mathrm{eff}}(x)`$ can be presented in the following form,

$$V_{\mathrm{eff}}(x)=V(x)TS(x,T),$$
(19)
where the “entropy potential” $`S(x,T)`$ is defined by the expression

$$S(x,T)=\mathrm{ln}_{\mathrm{}}^+\mathrm{}𝑑y\mathrm{exp}\{[U(y)+v(x,y)]/T\}.$$
(20)
Notice that $`S(x)`$ does not depend on $`V(x)`$.

## III Applications

### A Pure entropy barriers

Let $`V(x)=0`$ in Eq. (18),

$$U(y)=\frac{1}{2}m\omega _1^2y^2,$$
(21)
and

$$v(x,y)=\frac{1}{4}m(\omega _2^2\omega _1^2)(1\mathrm{cos}x)y^2,$$
(22)
so that the atomic motion is inactivated in the $`x`$ direction, but the frequency of transverse oscillation depends on $`x`$, $`\omega =\omega _1`$ at $`x=0`$ and $`\omega =\omega _2`$ at $`x=\pi `$. Then the integral in Eq. (20) can be easily evaluated analytically, and the entropy potential is given by the expression

$$S(x)=\frac{1}{2}\mathrm{ln}\left\{1+\frac{1}{2}\left[\left(\frac{\omega _2}{\omega _1}\right)^21\right](1\mathrm{cos}x)\right\}.$$
(23)
Notice that the entropy potential (23) does not depend on temperature, because both potentials (21) and (22) depend on $`y`$ in the same way ($`y^2`$). The function $`S(x)`$ is shown in Fig. 1. It is periodic with the period $`a=2\pi `$ and the height $`\epsilon _S=|\mathrm{ln}(\omega _2/\omega _1)|`$. The diffusion coefficient is given by $`D=D_fF`$, where the entropy factor $`F`$ depends on the ratio of frequencies $`z=\omega _2/\omega _1`$ only,

$$F(z)=[I_+(z)I_{}(z)]^1,$$
(24)
where

$$I_\pm (z)=\pi ^1_0^\pi 𝑑x\left[1+\frac{1}{2}(z^21)(1\mathrm{cos}x)\right]^{\pm 1/2}.$$
(25)
Equations (24,25) yield $`F(z)=(\pi /2)^2𝐊^1(\sqrt{1z^2})𝐄^1(\sqrt{1z^2})`$, where $`𝐊`$ is the complete elliptic integral of first kind. Near $`z1`$ the function $`F(z)`$ has the expansion $`F(z)1\frac{1}{8}(1z)^2`$, while at $`z0`$ it behaves as $`F(z)(\pi /2)^2\mathrm{ln}^1(4/z)`$. The function $`F(z)`$ is presented in Fig. 2. One can see that in the overdamped limit, the effect of entropy barriers is not too strong, in particular, even for $`\omega _1/\omega _2=0.1`$ the diffusion coefficient reduces comparing with the free-diffusion value by a factor of $`F(0.1)0.66`$ only. Indeed, although the height $`\epsilon _S`$ tends to infinity at $`z0`$, the width of barriers becomes very narrow and thus cannot strongly modify the diffusion coefficient.

In the underdamped case, on the contrary, the role of entropy barriers is essential. The dependence of the diffusion coefficient on the damping constant $`\eta `$ was obtained numerically and shown in Fig. 3. One can see that the function $`D(\eta )`$ exhibits a typical behavior of activated diffusion ($`D\eta ^1`$ at small and large frictions with a crossover between the limits) as might be expected from the shape of the entropy potential $`S(x)`$ of Fig. 1. In the overdamped limit the average jump length is equal to the period of the potential $`S(x)`$, $`\lambda 2\pi `$, while in the underdamped limit long jumps with $`\lambda /2\pi 1`$ play the dominant role as shown in Fig. 4 (in these simulations we assumed that the atom is trapped in a given well if it has sojourned in this well for a time lapse not shorter than $`(2\eta )^1`$ ). The effect of entropy barriers is even stronger than might be expected from the analogy with the energy barriers of the same height. For example, for the frequencies $`\omega _1/\omega _2=0.1`$ used in the simulation, the height of the barrier is $`\epsilon _S=S(\pi )2.3`$, that would give the ratio $`D(\eta \mathrm{})/D(\eta 0)2\epsilon _S/k_BT4.6`$ for the $`T=1`$ case, while the simulation leads to the ratio $`D(\eta \mathrm{})/D(\eta 0)>33`$. From Fig. 4b one can see that $`\lambda /2\pi 10^2`$ for the case of $`\eta =10^3`$, while for the one-dimensional diffusion it should be $`\lambda /2\pi \eta ^1=10^3`$. Thus, multi-dimensional effects result in a strong reduction of jump’s length in the underdamped limit which leads to a decrease of the diffusion coefficient comparing with the 1D motion. Note also that the dependence on temperature (shown in inset of Fig. 3) is almost negligible as has to be expected for the entropy potential.

### B Atom in a corrugated channel

Let now the dependence of the external potential $`V_N(x,y)`$ on $`y`$ is still given by Eqs. (21) and (22), but the motion in the $`x`$ direction is activated,

$$V(x)=\frac{1}{2}\epsilon (1\mathrm{cos}x),$$
(26)
where $`\epsilon `$ is the height of the external potential. At the minima of the potential (26) the transverse vibrations are characterized by the frequency $`\omega _1`$, while at the saddle points, by the frequency $`\omega _2`$. In the one-dimensional case, as well as for the 2D case with $`\omega _1=\omega _2`$, in the overdamped limit we have, according to Eq. (1), $`D_{\mathrm{Smoluchowski}}=D_f\mathrm{I}_0^2(\epsilon /2T)`$. Because the entropy potential $`S(x)`$ does not depend on the function $`V(x)`$, it is still given by Eq. (23), and the integral (16) can be easily evaluated. The results for the overdamped limit are shown in Fig. 5, which can be compared with the simulation results for different frictions presented in Figures 6 and 7.

From the $`D(T)`$ dependence of Figures 5 and 6 one can see that at high temperatures, when the motion is inactivated, the $`\omega _1=\omega _2`$ case leads to the maximum of the diffusion coefficient similarly to the case with pure entropic barriers. With temperature decreasing, the energy barriers and the entropy barriers play “in phase” for the “narrow-barriers” case of $`\omega _1<\omega _2`$, and “in antiphase” for the “wide-barriers” case of $`\omega _1>\omega _2`$. At low temperatures $`D>D_{1\mathrm{D}}`$ for the case of $`\omega _1>\omega _2`$ at high and moderate frictions in agreement with predictions of the TST approach. The effect, however, is smaller than the TST predicts: in simulation we found that the diffusion coefficient changes only in three times when the ratio of frequencies is equal to ten. At very low frictions (e.g., $`\eta <10^2`$ in Fig. 7), the entropy barriers become more important than the energy barriers, and the diffusion coefficient again becomes smaller than the 1D one for all cases of $`\omega _1\omega _2`$ as it was for the case of pure entropic barriers.

For moderate and low frictions the simulation results of Fig. 7 can be fitted by a dependence $`D(\eta )\eta ^{1/3}`$. Because the escape rate $`R`$ is still proportional to $`\eta `$ in the multi-dimensional case , we may conclude that in the present case, the average jump length scales as $`\lambda \eta ^{2/3}`$, which is in agreement with the results of pure entropic barriers presented in the previous subsection, and also may be compared with the 1D law $`\lambda \eta ^1`$ and the 2D simulation result  $`\lambda \eta ^{3/4}`$. Thus, in the underdamped limit multi-dimensional effects lead to decreasing of diffusivity (comparing with the 1D case) due to reduction of jump length which scales as $`\lambda \eta ^{2/3}`$ instead of the 1D scaling law $`\lambda \eta ^1`$.

### C Diffusion of the dimer

Now we can study diffusion of a dimer in the 1D sinusoidal potential. Let $`x_1`$ and $`x_2`$ are the coordinates of two atoms coupled by the elastic spring with the constant $`g`$, and $`a_0`$ is the equilibrium distance ($`0a_0\pi `$). Then the Hamiltonian of the system takes the form

$`H={\displaystyle \frac{1}{2}}m_a\dot{x}_1^2+{\displaystyle \frac{1}{2}}m_a\dot{x}_2^2+{\displaystyle \frac{1}{2}}\epsilon _s(1\mathrm{cos}2\pi x_1/a_s)+`$ (27)
$`{\displaystyle \frac{1}{2}}\epsilon _s(1\mathrm{cos}2\pi x_2/a_s)+{\displaystyle \frac{1}{2}}g(x_2x_1a_0)^2.`$ (28)
In what follows we put $`\epsilon _s=2`$, $`m_a=1`$, $`a_s=2\pi `$, and in the present paper we consider the case of $`a_0=0`$ only. Introducing the coordinates $`x=x_1+x_2`$ and $`y=x_2x_1`$, the Hamiltonian (28) can be rewritten as

$$H=\frac{1}{2}m\left(\dot{x}^2+\dot{y}^2\right)+V_N(x,y),$$
(29)
$$V_N(x,y)=\frac{1}{2}\epsilon \left(1\mathrm{cos}\frac{x}{2}\mathrm{cos}\frac{y}{2}\right)+\frac{1}{2}gy^2,$$
(30)
which describes the motion of one particle of mass $`m=m_a/2=1/2`$ in the $`x`$-periodic potential of height $`\epsilon =2\epsilon _a=4`$ and period $`a=2a_s=4\pi `$.

The adiabatic trajectory for this system was studied in . Its shape depends on a value of the elastic constant $`g`$. The points $`(x,y)=(4\pi n,0)`$, where $`n`$ is an integer, always correspond to the absolute minimum of the potential energy. Near the minimum, the potential energy has the expansion $`V_N(x,y)\frac{1}{2}m\left(\omega _{0x}^2x^2+\omega _{0y}^2y^2\right)`$ with $`\omega _{0x}=1`$ and $`\omega _{0y}=(2g+1)^{1/2}`$. For a strong spring, $`g1/2`$, there is only one saddle point at $`(x_s,y_s)=(2\pi ,0)`$ between two adjacent minima $`(0,0)`$ and $`(4\pi ,0)`$. Near the saddle, the potential energy has the expansion

$$V_N(x,y)\epsilon _s+\frac{1}{2}m\left[\omega _{sx}^2(xx_s)^2+\omega _{sy}^2(yy_s)^2\right]$$
(31)
with $`\omega _{sx}=1`$ and $`\omega _{sy}=(2g1)^{1/2}`$, so that the activation energy for dimer motion is equal to $`\epsilon _s=\epsilon =4`$ (see Fig. 8). Therefore, dimer diffusion can be approximately described as motion of one atom in the corrugated periodic potential with the transverse frequencies $`\omega _{1,2}=(2g\pm 1)^{1/2}`$, i.e. it corresponds to the case of “wide” barriers studied in the previous subsection. Thus, although the shape of adiabatic trajectory does not depend on the elastic constant for the case of strong coupling, the diffusion coefficient does depend on $`g`$, it increases when $`g1/2`$ due to decreasing of the transverse curvature at the saddle point. The simulation results of Fig. 9 show that the harmonic approximation describes the $`D(g)`$ dependence with a good accuracy. From Fig. 10, where the ratio $`D(g)/D(0)`$ is presented for different temperatures, one can see also that close to the critical point $`g=1/2`$, when anharmonicity of transverse vibrations at the saddle point is large, the entropy factor strongly depends on $`T`$, especially at low temperatures.

For intermediate values of the elastic constant, $`1/\pi g<1/2`$, the adiabatic trajectory still has only one saddle point $`(2\pi ,y_s)`$ between the adjacent minima, where $`y_s`$ is now a solution of the transcendental equation $`\mathrm{sin}(y_s/2)=gy_s`$. Near the saddle, the potential energy has the expansion (31) with the frequencies $`\omega _{sx}=\left[1(gy_s)^2\right]^{1/4}`$ and $`\omega _{sy}=(2g\omega _{sx}^2)^{1/2}`$. The saddle is characterized by the energy $`\epsilon _s(g)=\frac{1}{2}\epsilon \left[1+\mathrm{cos}(y_s/2)\right]+\frac{1}{2}gy_s^2`$, so that $`2+\pi /2<\epsilon _s<4`$.

Finally, for a weak coupling between dimer’s atoms, $`g<1/\pi `$, there are two saddle points between the adjacent minima $`(0,0)`$ and $`(4\pi ,0)`$, with a local minimum of the potential energy between these saddle points. The coordinates of the saddle points are $`(2\pi x^{},\pi )`$ and $`(2\pi +x^{},\pi )`$, where $`x^{}=2\mathrm{cos}^1(g\pi )`$. These saddle points are characterized by the energy $`\epsilon _s(g)=\frac{1}{2}(\epsilon +g\pi ^2)`$, so that $`2<\epsilon _s<2+\pi /2`$. Near the saddle, the potential energy has the expansion (31) with coefficients $`\omega _{sx}=(gG)^{1/2}`$ and $`\omega _{sy}=(g+G)^{1/2}`$, where $`G=\left[1(\pi ^21)g^2\right]^{1/2}`$.

The whole dependence $`\epsilon _s(g)`$ is shown in Fig. 8a. The activation energy monotonically increases from the single-atom value $`\epsilon _s=2`$ at $`g=0`$ to the rigid-dimer value $`\epsilon _s=4`$ at $`g=1/2`$ and then remains constant. Thus, one could expect that the diffusion coefficient should monotonically decrease with $`g`$ increasing. However, the simulation results of Fig. 10 show that often this is not true. The peculiarity in the transverse frequencies at the point $`g=1/2`$, where the saddle transverse frequency reaches zero, leads to a maximum of the function $`D(g)`$ close to this point, if the damping is small, $`\eta 0.5`$, and the temperature is not too low, $`T1`$ (recall $`\epsilon =4`$). Thus, multi-dimensional effects may strongly affect dimer’s diffusivity.

## IV Conclusion

In the present paper we studied in details the diffusion of a particle in two-dimensional space which is periodic along $`x`$ and unbounded in the transverse direction. We calculated the entropy factor which emerges due to transverse degree of freedom, both in the overdamped limit (analytically) and in the underdamped case (numerically), and compared it with the prediction of the transition-state theory. We showed that in the underdamped limit, the multi-dimensional effects lead to reduction (comparing with the one-dimensional motion) of jump lengths between subsequent trapping of the atom in bottoms of the external periodic potential. The simulation predicts that jump lengths scale as $`\lambda /\lambda _{1\mathrm{D}}\eta ^{1/3}`$. This leads to a decrease of diffusivity which now scales as $`D\eta ^{1/3}`$ instead of the 1D dependence $`D_{1\mathrm{D}}\eta ^1`$.

In the overdamped limit, the entropy factor (and, therefore, the prefactor in the Arrhenius formula for activated diffusion) does not depend on temperature as long as the transverse motion near the adiabatic trajectory may be described by the harmonic approximation. Simulation shows that this remains true, at least approximately, for low damping as well. Thus, in most cases the experimentally observed dependence of the prefactor on temperature has to be attributed to collective effects due to interaction between diffusing particles or/and between the atom and (deformable) substrate. However, in the case of dimer diffusion at some value of the interaction between the atoms, when the saddle transverse frequency is equal zero, the anharmonicity of the transverse potential begins to play the important role and the entropy factor strongly depends on $`T`$.

###### Acknowledgements.

Helpful discussions with T. P. Valkering are gratefully acknowledged. This work was partially supported by the INTAS Grant 97-31061.


text:
# An Efficient Modified “Walk On Spheres” Algorithm for the Linearized Poisson-Boltzmann Equation

## Abstract

A discrete random walk method on grids was proposed and used to solve the linearized Poisson-Boltzmann equation (LPBE) . Here, we present a new and efficient grid-free random walk method. Based on a modified “Walk On Spheres” (WOS) algorithm  for the LPBE, this Monte Carlo algorithm uses a survival probability distribution function for the random walker in a continuous and free diffusion region. The new simulation method is illustrated by computing four analytically solvable problems. In all cases, excellent agreement is observed.

Random walk methods have been used to solve a wide variety of parabolic and elliptic partial differential equations (PDEs) . Generally, there are two broad classes of random walk methods; one uses discrete random walks on grids , and the other continuous random walks in free space . One of the widely used continuous random walk methods, the “Walk On Spheres” (WOS) method , uses the first-passage probability distribution on a sphere to facilitate large steps in random walks. (The first-passage probability, $`w(𝐱;𝐱_\mathrm{𝟎})`$, is the probability of hitting the vicinity of $`𝐱`$ on the bounding surface for the first time when the random walker starts from $`𝐱_\mathrm{𝟎}`$, a point inside the bounding surface.) This continuous random walk method needs to discretize neither space nor time, nor the diffusing trajectory, and so it is particularly advantageous when the geometry of the region of interest is very complex or if the solution of the PDE is required at only a relatively small number of points.

We are interested in solutions to the Dirichlet problem for the linearized Poisson-Boltzmann equation (LPBE) in the domain $`\mathrm{\Omega }`$:

$$^2\psi (𝐱)=\kappa ^2\psi (𝐱),𝐱\mathrm{\Omega },$$
(1)
$$\psi (𝐱)=\psi _0(𝐱),𝐱\mathrm{\Omega },$$
(2)
where $`\kappa `$ is called the inverse Debye length . Notice that when $`\kappa ^2`$ is zero, the above problem becomes a Dirichlet problem for the Laplace equation. The WOS method for the Dirichlet problem for the Laplace equation has been widely used . We will combine this WOS method with a survival probability density function which incorporates the term involving $`\kappa ^2`$ in the LPBE.

In a discrete random walk method  for the LPBE, the corresponding Master equation relates $`\kappa ^2`$ to the removal probability of the random walker on the grid. During each step of the discrete random walk, the walker either moves to one of the neighboring sites, or stays fixed, or is removed. This probabilistic interpretation of $`\kappa ^2`$ can be also extended to continuous random walk methods once we know the survival probability distribution function of a random walker in continuous space.

In this letter, we obtain this survival probability distribution function of a random walker in continuous space by reinterpreting the weighting function in the previous modified WOS method  for the LPBE (see the Appendix for more details). The survival probability of a random walker in a continuous and free diffusion region is given by :

$$p(d)=d\kappa /\mathrm{sinh}(d\kappa ),$$
(3)
where $`d`$ is the distance from the starting point in the diffusion region. Figure 1 shows this probability density function. We modify the WOS method to incorporate the survival probability to solve the LPBE via a continuous random walk method. This probability density combined with the WOS method is used to remove a random walker during the random walk by the acceptance-rejection method . We generate a random number, $`\eta `$ in $`[0,1)`$ when we perform a WOS step, and we compare $`\eta `$ with $`p(d)`$, the survival probability at $`d`$, the radius of WOS. If $`\eta >p(d)`$, the random walker is removed at this WOS step.

An estimate for the solution of the LPBE at $`𝐱_\mathrm{𝟎}`$, where random walkers start, is given by $`S_N`$:

$$S_N=\frac{1}{N}\underset{i=1}{\overset{N_s}{}}\psi _0(X_{n_i}).$$
(4)
Here, $`N`$ is the total number of random walkers, $`N_s`$ is the number of survived-and-absorbed random walkers, and $`X_{n_i}`$ is the final position of the walker on the boundary when it is absorbed after $`n_i`$ WOS steps.

In this method, like the WOS method, errors are due to both statistical sampling and the $`\delta `$-absorption layer which captures random walkers near the boundary to terminate their random walk. However, the error from the $`\delta `$-absorption layer can always be made smaller than the statistical error . For the same random walk, the estimate difference between using $`\delta `$ and $`\delta /10`$ gives a measure of the error due to the finite width of the $`\delta `$-absorption layer. By adjusting $`\delta `$ we can make the error from the absorption layer less than the statistical error. This means that if we increase the number of random walkers to decrease the statistical error, consequently we must reduce $`\delta `$ and so increase the running time.

In the following, we compare our simulation results with the analytic results for four problems, which were used as examples for the discrete random walk method . In all cases, the results are given as those normalized by the boundary condition $`\psi _0`$, which is assumed sufficiently small for the LPBE to be valid. The number of random walks used for the solution at a point is $`10^5`$, and the absorption layer thickness is $`\delta =10^4`$. The analytic results  are shown with solid lines in Figs. 234 and 5 and our simulation results with circles. For the all four cases, our simulation results show excellent agreement.

Our method has several features. First, it is easier to implement and will be faster than the other discretized methods, such as the discrete random walk method , the finite difference method  and the boundary-element method , especially with complicated geometries. For a desired point, it takes only a few seconds to compute a solution with $`10^5`$ random walkers and $`\delta =10^4`$ on a 550 Mhz PC. However, it is hard to compare to other methods because they compute solutions at all grid points. We can safely say that continuous Monte Carlo methods are more efficient when the solution is required only at relatively small number of points. Secondly, the accuracy and the running time of our method depends primarily on the number of statistical samples, and so it is naturally parallel. Thirdly, it is certain that our new method is faster than the old modified WOS method , because while some of our random walkers are removed during their random walk, in the old method all random walkers must complete their random walks to contribute to the solution according to their weightings. Also, in open boundary cases, like the three examples except the parallel plates, it is necessary to use a certain cut-off in the old modified WOS  to kill random walkers, which will bias the results. As an example, in Table 1 in the case of the parallel plates, we compare our new method with the old modified WOS method, . We use the customary comparison method for Monte Carlo methods, the time consumption (or laboriousness) : $`t\times \underset{¯}{\text{D}}\xi `$, where $`t`$ is the CPU time expended in calculating a single estimate and $`\underset{¯}{\text{D}}\xi `$ is the variance of the estimates. The less laborious the algorithm, the more efficient it is. In Table 1, the time consumption (or laboriousness) of our algorithm is better than that of the old modified WOS method. Finally, our method is easy to extend to solve the LPBE with source terms . That will be the subject of our upcoming research with biochemical applications.

APPENDIX

In this appendix, we show how the weighting function in the old modified WOS method  can be interpreted as the survival probability distribution function. For simplicity, consider the LPBE in the old modified WOS method . The solution at $`𝐱_\mathrm{𝟎}`$ in the domain can be expressed as follows :

$$u(𝐱_\mathrm{𝟎})=\frac{1}{N}\underset{i=1}{\overset{N}{}}Q_i^{n_i}\psi _0(X_{n_i}),$$
(5)
where

$$Q_i^0=1,Q_i^{n_i}=Q_i^{n_i1}\frac{d_i^{n_i1}\kappa }{\mathrm{sinh}(d_i^{n_i1}\kappa )},d_i^{n_i}=d(P_i^{n_i}).$$
(6)
Here, $`N`$ is the total number of diffusing random walkers, $`i`$ refers to $`i`$th random walker, $`X_{n_i}`$ is the position where the $`i`$th random walker is absorbed in the $`\delta `$-absorption layer after $`n_i`$ WOS steps, and $`d_i^{n_i}`$ the radius of $`n_i`$th WOS of the $`i`$th random walker.

If we interpret $`Q_i^{n_i}`$ as a survival probability of $`i`$th random walker, $`_{i=1}^NQ_i^{n_i}`$ is the total number of survived-and-absorbed random walkers. Furthermore, due to the property of probabilistic random sampling from the total random walkers, only the survived-and-absorbed random walkers can be regarded as contributors to the solution. This reinterpretation of the weighting function as the survival probability distribution function is a kind of the fractional sampling method, i.e. ‘Russian Roullete’,  which has been used extensively in neutron transport and similar problems.


text:
# 1 RG in large extra dimensions

## 1 RG in large extra dimensions

Recently the extra compactified dimensions have been attracting much attention as possibilities to explain various hierarchy problems; the gauge hierarchy, the Yukawa hierarchy and so on. The effective field theories in extra dimensions contain towers of massive Kaluza-Klein (KK) excitations, whose quantum effects alter the behavior of the running couplings from logarithmic to power. Therefore the traditional picture of gauge coupling unification may be drastically changed . In this talk we examine the running coulings explicitly in several schemes and consider the origin of power law running and their scheme (in)dependence .

In practice the notion of running couplings is not well defined in extra dimensions due to nonrenormalizability. We cannot help but defining the field theories by a certain cutoff. Then the Wilson RG is supposed to offer a natural framework to define the RG flows in such cases. There is another reason why the Wilson RG is suitable for the calculation of the $`\beta `$-functions in extra dimensions. The power law behavior of the running coupling is supposed to be generated by successive threshold corrections by tower of massive KK modes. The Wilson RG is faithful to their decoupling effects.

Therefore we apply the Exact RG, which is a continuum formulation of the Wilson RG, to a scalar field theory in $`M_4\times S^1`$ space-time. In this formulation a certain cutoff is performed to the internal loop momenta, and the running coupling is defined by variation of the cutoff scale. We derive the $`\beta `$-functions for the four scalar coupling in the 1-loop level, assuming that the coupling is weak enough. We also calculate the $`\beta `$-functions in other schemes: the proper time regularization and the momentum subtraction, and compare these results.

In all schemes we obtain the $`\beta `$-functions as

$$\beta _\lambda =\frac{d\lambda }{dt}=bϵ_k(t)\lambda ^2,$$
(1)
where scale parameter $`t=\mathrm{ln}R\mathrm{\Lambda }`$ is introduced in terms of radius of the compact space $`R`$ and the cutoff scale $`\mathrm{\Lambda }`$. $`b`$ is the 1-loop coefficient in four dimensions. The function $`ϵ_k(t)`$ is dependent on each scheme $`k`$, and it’s asymptotic form is given by

$$\begin{array}{ccccc}ϵ_k(t)\hfill & & 1\hfill & \mathrm{for}\hfill & t0,\hfill \\ & & B(k)e^t\hfill & \mathrm{for}\hfill & t0,\hfill \end{array}$$
(2)
where $`B_k`$ is a scheme dependent constant. This function is shown in Fig.1 by truncating the KK modes at $`N`$-th level.

It is seen that the $`\beta `$-function shifts to the five dimensional one smoothly by the KK threshold corrections. We may define the effective dimensions  by

$$D_{\mathrm{eff}}(t)=4+\frac{d\mathrm{ln}ϵ_k}{dt},$$
(3)
which shows smooth transition from 4 to 5.

## 2 Scheme dependence

The $`\beta `$-functions given in Eq. (1) are scheme dependent even in the 1-loop level. However we may redefine the coupling constants in different schemes so that the $`\beta `$-functions match in the asymptotic region. The scheme dependence remained even after this procedure represents the net ambiguity in the RG in extra dimensions. The results obtained by explicit calculations are shown in Fig. 2. It is seen that the unremovable scheme dependence is rather small. In consequence, it is found that the GUT predictions for the low energy gauge couplings  are almost scheme independent .

## 3 RG in the warped extra dimensions

The Randall-Sundrum type of compactification  is now under the most vigorous investigation. The extra dimension is bounded by two three branes, in one of which the Planck scale is reduced to TeV scale. A peculiarity of this compactification is that the mass scale of the KK tower appears at TeV order, even if the bare mass is Planck scale . If we apply the RG scheme mentioned above to this case, the running couplings are found to show power law behavior too. This should be compared with the results obtained by the PV regularization . However it should be said that the running coupling above TeV scale is totally regularization dependent, since the KK spectrum is influenced by the Planck scale (string) physics.


text:
# 1 Overview

## 1 Overview

Anomalies correspond to the fact that sometimes symmetries of a classical theory cannot be implemented in the corresponding quantum theory. The study of field theoretical anomalies has a long history. They were originally observed as unexpected terms in the commutators of quantum currents among each others (Schwinger terms) or with products of field operators (anomalous Ward identities). Much insight was subsequently obtained in the nature of anomalies by topological and cohomological methods.

One line of development lead to a description of gauge anomalies as local cocycles on the gauge group expressed in terms of field operators. Relating anomalies to the lack of invariance of the determinant of operators of Dirac type yielded an understanding in terms of characteristic classes of determinant bundles over the space of connections.

These characteristic classes are all obtainable from first quantization of a Dirac type particle in an external classical gauge field, either in a Lagrangian or in a Hamiltonian formulation. Sometimes it is also appropriate to consider bundles of Fock spaces over the space of connections (). In any case, the information on anomalies is related to the bundles of null spaces of families of Dirac operators. The anomalies are traces of products of inverse powers of the Dirac operator $`D`$ and certain vertex operators, just those expressions also obtainable from one-loop Feynman integrals, where $`D^1`$ corresponds to the propagator of the Dirac field.

In yet another language, the anomalies are obtained from the symbol calculus of differential operators on the Dirac bundle. Further insight is gained from non-commutative geometry. Starting from the algebra of differential operators on the Dirac bundle, non commutative characteristic classes can be found, which are just the non commutative counterpart of the anomalies and which, evaluated on suitable non-commutative cycles, just give the characteristic numbers.

These results can also be looked at from the point of view of star products or deformation quantization (). The algebra $`𝔄`$ of symbols of differential operators on the Dirac bundle can be considered as the algebra of observables of a system consisting of a Dirac particle in an external gauge field. The symbol calculus in this algebra is thus a special example of a star product on $`𝔄`$. More precisely, one has to go over to the algebra $`𝔄[[\lambda ]]`$ of formal power series in an indeterminate $`\lambda `$ (eventually to be substituted by $`\mathrm{}`$) with coefficients in $`𝔄`$.

This restriction to formal power series is by no means only a shortcoming, but it also has definite advantages.

It opens the way to a concentration on a better understanding of algebraic aspects of the quantization procedure. In particular, there is the powerful notion of equivalence of star products that are related, for example, by changing ordering prescriptions in the quantization procedure. Anomalies, anyhow, are always of low order in $`\mathrm{}`$ and should be treatable without any loss of information in terms of formal power series in $`\lambda `$. In the language of star product quantization, anomalies are objects in the non-commutative geometry of the algebra $`𝔄[[\lambda ]]`$ endowed with the star product. The trace functional on $`𝔄[[\lambda ]]`$, applied to non-commutative characteristic classes on $`𝔄[[\lambda ]]`$ give characteristic numbers. Uniqueness theorems on the trace functional and non commutative cohomologies can successfully be applied. The theory of star products on finite dimensional phase spaces like the phase space of the Dirac system is well defined. Anomalies are expected also to arise as imperfections in invariance properties of the star product.

## 2 A mechanical example

Let $``$ be a classical phase space of finite dimension. Let a group $`𝒢`$ with Lie algebra $`𝔤`$ act on $``$ in a Hamiltonian way with equivariant momentum map ()

$$𝒥:𝔤^{},$$

such that for every $`\xi 𝔤`$, $`fC^{\mathrm{}}()`$

$$L_\xi f=\{f,𝒥_\xi \}\text{and}\{𝒥_\xi ,𝒥_\eta \}=𝒥_{[\xi ,\eta ]},$$

where $`𝒥_\xi (p)=𝒥(p),\xi `$.

A star product $``$ on $`C^{\mathrm{}}()[[\lambda ]]`$ is called (see, e.g., )

1. covariant, if $`[𝒥_\xi ,𝒥_\eta ]_{}=𝒥_\xi 𝒥_\eta 𝒥_\eta 𝒥_\xi =\text{i}\lambda 𝒥_{[\xi ,\eta ]}`$,
2. invariant, if $`\{𝒥_\xi ,fg\}=\{𝒥_\xi ,f\}g+f\{𝒥_\xi ,g\}`$,
3. strongly invariant, if $`𝒥_\xi ff𝒥_\xi =\text{i}\lambda \{𝒥_\xi ,f\}`$.

Clearly, the last property implies the two others. The second equation is the infinitesimal version of the requirement of invariance under canonical transformations $`\mathrm{\Phi }`$:

$$\mathrm{\Phi }^{}(fg)=\mathrm{\Phi }^{}f\mathrm{\Phi }^{}g.$$

Schwinger terms should correspond to a violation of the covariance condition (1.), whereas anomalous terms in Ward identities are expected to reflect themselves in a lack of strong invariance (3.).

The group of all classical canonical transformations, for example, cannot be implemented in quantum theory without anomalies. Indeed, for this large group, all observables are momentum maps, conditions (1.) and (2.) become identical and the impossibility of a star product to transform Poisson brackets to star commutators directly follows from the theorem of Groenewald and van Hove. These anomalies are the reason for the notorious incompabilities of quantization procedures and canonical transformations. If conditions (1.), (2.) or (3.) are not fulfilled, one may try to save the situation by adding non leading terms to $`𝒥_\xi `$,

$$𝕁_\xi =𝒥_\xi +\lambda (\mathrm{}),$$

such that $`𝕁_\xi `$, the so-called quantum momentum map (), shows the above properties.

For example, quantum covariance means

$$𝕁_\xi 𝕁_\eta 𝕁_\eta 𝕁_\xi =\text{i}\lambda 𝕁_{[\xi ,\eta ]}.$$

One easily convinces oneself, that covariance holds for a star product $``$ if and only if it holds for all equivalent star products. Hence, the transformation

$$T:𝒥_\xi 𝕁_\xi $$

rendering a non-covariant $`𝒥_\xi `$ covariant cannot be an equivalence transformation but must depend on $`\xi `$. The redefinition of $`𝒥_\xi `$ is ad hoc rather than global.

As an illustration we can consider a mechanical system in one space dimension. Galilei invariant systems are not appropriate, because there is no equivariant momentum map. The generators of translations and boosts $`p`$ and $`q`$ commute as Lie algebra elements, but as generators of canonical transformations we have $`\{p,q\}=1`$. The Galilei group can be realized by canonical transformations only after central extension.

The Poincaré group does not have this defect. The Ruijsenaars-Schneider model () is a mechanical system in one space dimension with Poincaré invariance and non trivial interaction.

Let us denote the generators of time and space translations by $`𝔥`$ and $`𝔭`$ and the boost generator by $`𝔟`$. In the Lie algebra of the Poincaré group we have

$$[𝔥,𝔭]=0,[𝔟,𝔥]=𝔭,[𝔟,𝔭]=\frac{1}{c^2}𝔥.$$

In the Ruijsenaars-Schneider model, the same relations hold,

$$\{𝔥,𝔭\}=0,\{𝔟,𝔥\}=𝔭,\{𝔟,𝔭\}=\frac{1}{c^2}𝔥,$$

where $`𝔥`$, $`𝔭`$, and $`𝔟`$ now are certain functions over phase space. This is not in contradiction with the no-interaction theorem, because the translations are realised in the conventional way only to order $`\frac{1}{c^2}`$. Moreover, in a canonical quantization one has

$$[,]=0,[𝔹,]=\frac{\text{i}}{\mathrm{}},[𝔹,]=\frac{\text{i}}{\mathrm{}c^2},$$

or (the factor $`\text{i}`$ can be absorbed in the formal parameter $`\lambda `$)

$$[\widehat{𝔥},\widehat{𝔭}]_{}=0,[\widehat{𝔟},\widehat{𝔥}]_{}=\widehat{𝔭},[\widehat{𝔟},\widehat{𝔭}]_{}=\frac{1}{c^2}\widehat{𝔥}$$

for the standard ordered star product $``$. Here, $`\widehat{𝔥}`$, $`\widehat{𝔭}`$, and $`\widehat{𝔟}`$ differ from the corresponding classical functions only in higher orders of $`\lambda `$. In addition the classical Ruijsenaars-Schneider system is completely integrable. It has a series of conserved quantities $`𝒥_k`$ in involution,

$$\{𝒥_k,𝒥_l\}=0.$$

In the quantized theory one has, as shown by Ruijsenaars (), that a modification $`𝒥_k𝕁_k`$ is possible, such that quantum integrability

$$[𝕁_k,𝕁_l]_{}=0$$

holds, where $``$ again is the standard ordered star product.

## 3 Field theory

One would, of course, like to apply the star product formalism directly to field theory, not only to first quantization. Unfortunately, for this case, the theory of star products is not yet fully developed.

One reason lies in the fact that even for the classical algebra of fields products like $`\phi _1(x_1)\phi _2(x_2)`$ are in general undefined for $`x_1=x_2`$, this deficiency beeing even more serious for quantum fields. The normal way out is to consider (quantum) fields as (operator valued) distributions. For objects

$$\phi _f=d^nx_1\mathrm{}d^nx_r\phi (x_1)\mathrm{}\phi (x_r)f(x_1,\mathrm{},x_r)$$

with suitable test functions $`f`$, products are now defined ( and references therein). For the same reason, one expects star products $`\phi _f\phi _g`$ to be well defined.

In quantum field theory, one normally concentrates on expectation values like

$$w(x_1,\mathrm{},x_n)=0|\phi (x_1)\mathrm{}\phi (x_n)|0.$$

This should also be a promising strategy for field theory in star product formulation. One can interprete the distribution

$`w(x_1,\mathrm{},x_n)`$ $`=0|\phi (x_1)\mathrm{}\phi (x_n)|0`$
$`=\mu \left(\phi (x_1)\mathrm{}\phi (x_n)\right)`$
as the value of a state on the algebra $`𝒜[[\lambda ]]`$ of classical fields endowed with a star product. In fact, $`w`$ even resembles a trace functional on $`𝒜[[\lambda ]]`$. $`w`$ has to be interpreted as a formal series in $`\lambda `$, whose coefficients are just the quasi-classical expansion of the quantum correlation functions. This will result in considerable technical simplifications as compared to quantum field theory. The star product algebra can now be reconstructed by a generalized GNS construction in the sense of the method introduced by M. Bordemann and S. Waldmann (). Recently, part of this program has been implemented by K. Fredenhagen and M. Dütsch in .

For the discussion of anomalies, the momentum map $`𝒥_\xi `$ has to be replaced by the conserved current $`𝒥_\xi ^\mu `$. Adding higher terms to find a quantum momentum map corresponds to the need of regularization. Indeed, by using the canonical (anti)commutation relations the classical expression for normal ordered currents is found formally only to order $`\mathrm{}`$.

Conditions similar to (1.), (2.), and (3.) can be formulated, and their violation will be largely analogous to anomalous Ward identities in quantum field theory but in the sense of formal power series. Violations of the above conditions should be characterized cohomologically.

There should be a formulation of anomalies in terms of the non-commutative geometry in the field algebra $`𝒜[[\lambda ]]`$ and not just in the algebra $`𝒜_D[[\lambda ]]`$ of a Dirac particle.

The Dirac Operator plays a fundamental rôle in the non-commutative geometry of $`𝒜_D[[\lambda ]]`$. $`D`$ can also be interpreted as a derivation in $`𝒜`$, which can also be written in the form

$$𝒜[[\lambda ]]a[Q,a]_{}𝒜[[\lambda ]],$$

with

$$Q=d^nx\overline{\phi }D\phi .$$

A more ambitious task would be to formulate anomalies as algebraic properties on the algebra $`[[\lambda ]]`$ of observables of a field theory rather than on the larger field algebra $`𝒜[[\lambda ]]`$. The Dirac operator will still be a derivation in $`[[\lambda ]]`$. Gauge anomalies primarily concern non observable quantities, their effect should be hidden more deeply in $`[[\lambda ]]`$.

## Acknowledgements

This note is dedicated to H. D. Doebner.
We thank M. Bordemann for many elucidative discussions. H. R. thanks J. Mickelsson for hospitality at the Royal Technical Highschool Stockholm and for useful discussions.

The revised version was induced by learning about a paper by Dütsch and Fredenhagen (), in which important elements of the construction sketched here were implemented in considerable detail. We also thank M. Dütsch and K. Fredenhagen for illuminating discussions.