text:
# Determination of Primordial Metallicity and Mixing in the Type IIP Supernova 1993W1footnote 11footnote 1Partially based on observations performed at ESO La Silla

## 1 Introduction

Supernova 1987A (SN 1987A) confirmed our basic theoretical understanding of a Type II supernova (SN II) as the core collapse of a massive star, which leaves behind a compact object (neutron star or black hole) and expels the outer mantle and envelope into the interstellar medium (see Arnett et al., 1989, and references therein). While the explosion mechanism itself remains a subject of active research, the propagation of the shock wave through the mantle and envelope is reasonably well understood. Theoretical models for the light curve do a good job of reproducing the observations (Blinnikov, 1999; Blinnikov et al., 2000) and detailed radiation transport calculations verify that these models reproduce the observed spectra from the UV to the IR (Mitchell et al., 2001, 2002). Unlike Type Ia supernovae (SNe Ia) where the actual compositions as a function of velocity are an important subject of current research (Fisher et al., 1999), the compositions of the envelopes of red or blue supergiant stars is primarily hydrogen and helium and thus we show that the primordial compositions, and the degree of mixing of hydrogen and 56Ni, can be determined to high precision by the detailed spectral modeling of observed SN II spectra.

SNe II have a very large spread in their intrinsic brightness, from the very dim SN 1987A, to the exceedingly bright SNe 1979C and 1997cy. The observed spread in intrinsic luminosity is greater than a factor of 500. This is not surprising given the fact that the progenitors span a wide range of initial stellar masses, possible binary membership, and prior star formation histories. Clearly, SNe II do not meet the astronomers requirement of being a *standard candle*, however we believe that the fact that our models can determine the stellar compositions, degree of mixing, and kinetic energy of the explosion shows that their atmospheres can be well understood. Although we do not present distances in this paper, we believe that these results increase the attractiveness of SNe II as cosmological probes and in the future SNe II will become complementary with SNe Ia as distance indicators. Both the SEAM (see Baron et al., 1995b, 1996; Mitchell et al., 2002, and references therein) and the EPM (see Baade, 1926; Schmidt et al., 1994; Hamuy et al., 2001; Leonard et al., 2002) methods for determining distances to SNe II (and other supernova types for SEAM) depend on the ability to model the spectral energy distribution (SED) of the supernova atmosphere accurately. Clearly, quantities like the composition and degree of mixing play a role in the output SED and thus the work presented here helps to place both these methods on firmer footing.

While the average SN II is several times dimmer than a SN Ia, current ground-based searches and proposed future space-based searches for supernovae will easily detect these objects at cosmologically interesting distances. In the very deep high-redshift supernova searches by the Hi-Z Team (Riess et al., 1998) and the Supernova Cosmology Project (SCP) (Perlmutter et al., 1999) several SNe II have been discovered even though the searches are geared to finding SNe Ia. The SCP has positively identified via spectroscopy 5 SNe II out to a redshift $`z0.45`$. Gilliland et al. (1999) detected SN 1997fg through differencing two epochs of Hubble Deep Field North $`I`$-band photometric images. While the supernova was not spectroscopically classified, the combination of the host galaxy type (an irregular), its redshift $`z=0.952`$, and the supernova’s apparent brightness makes it a strong candidate for a SN II. A. Clocchiatti et al. (IAUC 7549) found SN 2000fp at $`z=0.30`$, G. Altavilla et al. (IAUC 7762) found SN 2001gg at $`z=0.61`$, SN 2001gh at $`z0.16`$, and SN 2001gj at $`z=0.27`$, and G. Altavilla et al. found SN 2002co at $`z=0.318`$. The proposed *SNAP* satellite, a wide-field optical imaging telescope, (see http://snap.lbl.gov) will easily be able to find, follow, and spectroscopically identify SNe II beyond redshifts of $`z1`$. *NGST*, a large, small field IR telescope, will be able to detect SNe II to the initial star formation period in the universe, $`z>5`$.

The explosion mechanism is now thought to result in a non-spherically symmetric shock-wave (LeBlanc & Wilson, 1970; Khokhlov et al., 1999) and this has been confirmed by the detection of significant polarization in the spectra of SNe II (Wang et al., 1996; Leonard et al., 2000). Nevertheless, the large expansion of the ejecta before the supernova becomes visible leads to significant “sphericalization” (Chevalier, 1984), and asphericity effects may only be important at late times, or in SNe II that are strong circumstellar interacters (SNe IIn). These supernovae are clearly distinguishable observationally from “normal” Type II supernovae.

Calculations (Miralda-Escudé & Rees, 1997; Dahlén & Fransson, 1999; Yungelson & Livio, 2000) indicate that SNe II may in fact be the first stellar objects visible in the universe and hence they serve as important probes for the star formation rate, the rate of chemical evolution by measuring their primordial abundances, and directly for the cosmological parameters. We show that the primordial abundances of the progenitor star can be reliably determined by detailed synthetic spectral modeling.

## 2 Calculations

We have chosen to model SN 1993W in detail because there are several high signal to noise spectra very near the time of explosion as well as at later times and because we found that this object is clearly a low metallicity supernova although both the photometry (see Table 1) proves that it is indeed an SN IIP, closely resembling the light curve of SN 1969L and SN 1988A.

The calculations were performed using the multi-purpose stellar atmospheres program PHOENIX version 11.7 (Hauschildt & Baron, 1999; Baron & Hauschildt, 1998; Hauschildt et al., 1997a, b, 1996). PHOENIX solves the radiative transfer equation along characteristic rays in spherical symmetry including all special relativistic effects. The non-LTE (NLTE) rate equations for many ionization states are solved including the effects of ionization due to non-thermal electrons from the $`\gamma `$-rays produced by the radiative decay of 56Ni, which is produced in the supernova explosion. The atoms and ions calculated in NLTE are: H I, He I–II, C I-III, N I-III, O I-III Na I-III, Mg II, Ca II, Si I–III, S I–III, Fe I–III, Ni I-III, and Co I-III. These are all the elements whose features make important contributions to the observed spectral features in SNe II.

Each model atom includes primary NLTE transitions, which are used to calculate the level populations and opacity, and weaker secondary LTE transitions which are are included in the opacity and implicitly affect the rate equations via their effect on the solution to the transport equation (Hauschildt & Baron, 1999). In addition to the NLTE transitions, all other LTE line opacities for atomic species not treated in NLTE are treated with the equivalent two-level atom source function, using a thermalization parameter, $`\alpha =0.05`$. The atmospheres are iterated to energy balance in the co-moving frame; while we neglect the explicit effects of time dependence in the radiation transport equation, we do implicitly include these effects, via explicitly including the rate of gamma-ray deposition in the generalized equation of radiative equilibrium and in the rate equations for the NLTE populations.

The models are parameterized by the time since explosion and the velocity where the continuum optical depth in extinction at 5000 Å ($`\tau _{\text{std}}`$) is unity, which along with the density profile determines the radii. This follows since the explosion becomes homologous ($`vr`$) quickly after the shock wave traverses the entire star. The density profile is taken to be a power-law in radius:

$$\rho r^n$$

where $`n`$ typically is in the range $`610`$. Since we are only modeling the outer atmosphere of the supernova, this simple parameterization agrees well with detailed simulations of the light curve (Blinnikov et al., 2000) for the relatively small regions of the ejecta that our models probe.

Further fitting parameters are the model temperature $`T_{\text{model}}`$, which is a convenient way of parameterizing the total luminosity in the observer’s frame. We treat the $`\gamma `$-ray deposition in a simple parameterized way, which allows us to include the effects of nickel mixing which is seen in nearly all SNe II. Detailed fitting of a time series of the observed spectra determines all the parameters, i.e., $`n`$, $`T_{\text{model}}`$, the amount of nickel mixing, the amount of helium mixing, and the metallicity.

## 3 Results

SN 1993W was discovered on Aug. 19, 1993 by C. Pollas (IAUC 5848). Spectroscopy was obtained on Aug. 20 and Aug. 24 at ESO; and on Sept. 10, Sept. 11, Oct. 22, and Nov. 8 at Lick Observatory.

Figures 1 and 2 display the observed and best fit synthetic spectra for the two early epochs that we will study here. The observed spectra were de-redshifted by 5400 km s-1 and de-reddened using the reddening law of Cardelli et al. (1989) with $`E(BV)=0.047`$, i.e., accounting for galactic foreground extinction. There is no evidence of interstellar reddening in the parent galaxy since no narrow Na ID lines are visible in the observed spectra. Table 2 displays the model parameters used for each date. The fit quality is superb, no line deviates from the observed feature by more than 5% in the Aug 20 spectrum and the deviation is below 50% (except in Ca II which we discuss below) in the Aug 24 spectrum, which is excellent. Table 3 displays the synthetic colors for the best fit models and Table 4 lists the observed values. For Aug. 24, the synthetic $`BV`$, agrees with the observed colors to within 0.04 mag which is within the observational uncertainties. The $`UB`$ colors differ somewhat more, but $`U`$ photometry is often uncertain and when the Ca H+K line is well fit the synthetic value of $`UB=0.20`$ (see Figure 16). Clearly, our ability to fit observed spectra in such detail gives us confidence in the physical parameters that we can extract from our models. Figure 3 displays the Balmer lines as a function of velocity for Aug 20, 1993 and Figure 4 displays both the Balmer lines and the weak Fe II $`\lambda `$5169 line which is often used in the expanding photosphere method (Schmidt et al., 1992; Eastman et al., 1996; Schmidt et al., 1992; Hamuy et al., 2001; Leonard et al., 2002). This result shows that the Spectral-fitting Expanding Atmosphere method (SEAM, Baron et al., 1993, 1994, 1995a, 1995b, 1996; Mitchell et al., 2002) fits not only the weak lines such as Fe II $`\lambda `$5169 but also the prominent SNe II Balmer lines. The velocities determined by detailed spectral analysis are extremely accurate (close to the instrumental resolution) and indicates that early spectra are quite important for SEAM analyses.

Besides velocity, the other important fitting parameter is the total observed bolometric luminosity which provides the outer boundary condition for the requirment of global energy conservation. It is convenient to parameterize the luminosity in terms of the model temperature $`T_{\text{model}}`$. Figures 5 and 6 show the effect of varying the the model temperature on the standard reference models. From Figure 5 it is clear that the 7000 K model is too cool, producing both too little flux in the blue and far too much in the red. The 9000 K model produces too much flux in the blue, but with the spectral coverage of the observation, it is difficult to distinguish it in the red from our standard 8000 K model (see Fig. 1). Thus we can obtain the model temperature at this single epoch to $`\pm 500`$ K. With larger spectral coverage we should be able to determine the model temperature to $`\pm 100`$ K. Fig. 6 shows that on Aug 24, the cool model produces poorer fits to the Balmer lines, produces a very strong Ca II IR triplet (which also depends somewhat on the metallicity, discussed below) and also does less well at reproducing the Ca II H+K line. The hotter model also does a poorer job of reproducing the Balmer lines, and although the Ca II IR triplet is weaker, as needed by the observations, the Ca II H+K line is still far too strong.

### 3.1 Aug 20

Because the Aug 20 spectrum is so hot and so early, it contains (in the optical) only lines of hydrogen and helium. Thus we may use this spectrum to determine the amount of hydrogen and helium mixing, although we must simultaneously determine the amount of non-thermal excitation due to the radioactive decay of 56Ni. We have calculated a very large number of models ($`500`$) and we find that no reasonable value of $`\gamma `$-ray deposition can produce the observed He I $`\lambda 5876`$ line without enhancing the helium abundance over the solar value. We find that in order to fit this feature there must significant mixing of hydrogen and helium in the outer envelope and the helium must also be excited by non-thermal $`\gamma `$-ray deposition which comes from the radioactive decay of 56Ni. A similar result was found for SN 1987A (Arnett et al., 1989; Eastman & Kirshner, 1989; Lundqvist & Fransson, 1996; Blinnikov et al., 2000; Mitchell et al., 2001, 2002) and for SN 1999em (Baron et al., 2000). Figures Figure 7 and 8 show the effects of varying the helium mass fraction. Figure 8 shows that with a fixed amount of $`\gamma `$-ray deposition, solar abundances of helium do not produce a strong enough He I $`\lambda 5876`$ feature. Figures 9 and 10 show the effect of varying the $`\gamma `$-ray deposition, with a constant $`Y=0.76`$. Even with enhanced helium abundances in the outer layers, the He I $`\lambda 5876`$ feature does not appear with no non-thermal ionization, therefore mixing of nickel (in order to produce $`\gamma `$-rays to excite the helium) into the outer layers is required. The mixing of nickel, hydrogen, and helium in SNe II is not well understood and the situation is unlikely to improve until the explosion mechanism is better understood. Therefore we use a very simple prescription for $`\gamma `$-ray deposition, the $`\gamma `$-rays are deposited locally by a constant mass fraction of nickel $`\text{X}_{\text{Ni}}`$ — the deposition function follows the density profile. We have used this prescription successfully in modeling other SNe II (Lentz et al., 2001; Mitchell et al., 2001, 2002). Both enhanced surface helium and the excitation of the helium by $`\gamma `$-ray deposition appear to be a common feature of SNe IIP and support the notion that fluid instabilities or MHD jets are intimately involved in the explosion mechanism in order to produce large amounts of mixing.

Finally, Figure 11 displays the effects of varying the density profile $`n`$. The Balmer profiles are too narrow in the $`n=10`$ case and the $`n=6`$ model does nearly as well as our standard $`n=8`$ case. There is some degeneracy between the value of $`n`$ and the velocity, since the narrow lines in the $`n=10`$ case could be made broader by increasing the velocity of the ejecta, however in our large grid of models we found that lower values of $`n`$ are favored, largely due to the fact that changing $`n`$ and increasing the velocity alters the slope of the P-Cygni feature as it rises from the bottom of the absorption trough to the emission peak.

### 3.2 Aug 24

Figure 2 shows that the Fe II $`\lambda 5169`$ line has now appeared in the observed spectrum and a weak feature is also visible that is likely due to Fe II $`\lambda \lambda 4924,5018`$. However, there is not significant line blanketing in the blue even though the model temperature $`T_{\text{model}}`$ has dropped to the hydrogen recombination temperature (see Fig. 12). Figure 12 shows the effect of varying the model metallicity with gamma-ray deposition turned off. The results are expected, higher metallicity leads to large absorption due to iron line blanketing in the region 4000–4500 Å. In fact even at very low metallicity the line blanketing is much larger than observed, gamma-ray deposition significantly reduces the line blanketing as shown in Figure 13. In fact, the comparison of the two figures strengthen the need for significant gamma-ray deposition at early times. While Fig. 13 clearly shows that there is more line blanketing as the metallicity increases, it also shows that metallicity and $`\gamma `$-ray deposition radically alter the flux in the blue part of the spectrum.

Figure 14 shows that even at very low metallicity, the early spectra are also quite sensitive to the primordial metallicity, when gamma-ray deposition is included. The higher metallicity $`Z=0.10Z_{}`$ is favored to fit the two Fe II features in the blue (see Fig. 15), but the Ca II IR triplet is too strong with $`Z=0.10Z_{}`$, and is even a bit too strong for $`Z=0.03Z_{}`$. The Ca II H+K feature is too strong in both cases and is likely affected by the choice of gamma-ray deposition as can be seen from Figure 12. Since we have not tried to alter individual abundances, but only scaled total overall metallicity our preferred value for metallicity is set by the iron features in the blue. Clearly the progenitor star for SN 1993W had a significantly lower metallicity than the sun.

Fig. 16 shows a model with the same parameters as Figure 2, but with the Ca abundance reduced by a factor of 10. Now both the H+K and IR triplet lines are reasonably well reproduced in the synthetic spectrum; however since we have held everything else fixed the Balmer lines are not as well fit as they are in Figure 2. It is interesting to note that nucleosynthetic models for Type Ia and Type II supernovae from low mass progentitors produce Fe/Ca ratios that are about a factor of four higher than that produced by more massive SNe II (Tsujimoto et al., 1995). It is possible that the progenitor of SN 1993W just happened to occur in a region that was recently enriched by the explosion of a Type Ia supernova, resulting in enhanced Fe/Ca. Although low S/Fe and Si/Fe would also be expected in this scenario these are not ruled out by the spectra presented here. Si and S abundances can best be determined in later (nebular) spectra of SNe IIP. Since calcium lines are so strongly dependent on the temperature, small changes in the non-thermal ionization could possibly produce the suppression in the calcium strengths that we obtain here. A detailed study of a full time series of a similar SN IIP would thus be illustrative and we intend to carry out such a program in future work.

We find that $`0.03<Z_{\text{SN 1993W}}/Z_{}<0.15`$, and we strongly prefer a value of $`Z_{\text{SN 1993W}}/Z_{}0.10`$ (particularly if Fe/Ca is enhanced). The metallicity likely could be even better pinned down with UV observations that can be obtained by *HST* and by *SNAP*. By studying the progenitors of SNe II, Smartt et al. (2002), found that metallicities near solar with on SN 1987A and SN 1980K having metallicities of about $`0.5Z_{}`$. Thus the value of $`0.10Z_{}`$ for SN 1993W is interestingly low. Therefore early observations of SNe II are crucial if we wish to use these objects as probes of galactic chemical evolution.

## 4 Explosion Time

We can use simple arguments to obtain a limit on the explosion time. The total bolometric luminosity of the supernova is approximately

$`L`$ $`=`$ $`4\pi R^2\sigma T^4`$
$`=`$ $`4\pi v_{\mathrm{ph}}^2t^2T^4`$
where $`v_{\mathrm{ph}}`$ is the velocity of the photosphere. The ratio of the luminosities between time $`t_1`$ (Aug 20) and time $`t_2`$ is

$`{\displaystyle \frac{L_1}{L_2}}`$ $`=`$ $`{\displaystyle \frac{v_{\mathrm{ph}_1}^2t_1^2T_1^4}{v_{\mathrm{ph}_2}^2t_2^2T_2^4}}`$
$`=`$ $`{\displaystyle \frac{v_{\mathrm{ph}_1}^2t_1^2T_1^4}{v_{\mathrm{ph}_2}^2(t_1+3.90)^2T_2^4}}`$
where we have used the fact that $`t_2=t_1+3.90`$ days. From our synthetic spectra we can determine that the bolometric correction was $`BC=0.40`$ on Aug 20, and $`BC=0.04`$ on Aug 24. The $`B`$ magnitude is estimated to be 18.3 on both days (see Tab. 1) and using the $`BV`$ color from Table 3 we find that the ratio of $`L_1`$ to $`L_2`$ is 1.17. Plugging in our values for the velocities and temperatures we find that Aug 20 is $`4.7\pm 0.7`$ days after the explosion. The largest contributor to the uncertainty in the explosion time is the error in the observed photometry.

## 5 Conclusions

We have shown that synthetic spectra of SNe IIP can be fit to high accuracy, and that the primordial metallicity, degree of mixing of H and He, as well as the amount of mixing of 56Ni can also be determined by early spectra. Baron et al. (2000) showed that the total amount of extinction due to dust either in the parent galaxy or in our own galaxy can also be estimated with UV+optical spectra at very early times. While the reduced calcium abundance that was needed to reduce the strength of the calcium lines is intriguing, we do not as yet believe that our results require a reduced value of Ca/Fe. In a future paper we will explore a grid of models and determine a distance to SN 1993W, along with an estimate of the systematic error. We have shown that primordial metallicity, H/He mixing, and 56Nimixing are derivable from the detailed analysis of early spectra of SNe IIP. The analysis of abundances of individual elements is in reach, but will require a larger grid of models intercompared with a set of well studied SNe IIP. The data for this will be obtained as part of searches for nearby SNe Ia and will just require followup spectrophotometry. We believe that we have demonstrated that SNe IIP atmospheres can be well understood theoretically and hence are likely to make excellent independent cosmological probes, both as distance indicators and as probes of cosmic nucleosynthesis.

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 NASA grant NAG5-12127, NSF grant AST-0204771, and an IBM SUR grant to the University of Oklahoma; and by NSF grant AST-9720704, NASA ATP grant NAG 5-8425 and LTSA grant NAG 5-3619 to the University of Georgia. MT aknowledges the support of the Italian Ministry of University and Research through grant COFIN MM02905817. 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.


format: text

text:
# Wigner Centennial: His Function, and Its Environmental Decoherence

## 1 Introduction

Wigner, in his paper On the Quantum Correction For Thermodynamic Equilibrium , constructs a map from the quantum state of a particle into a certain classical phase space distribution:

$$\widehat{\rho }W(x,p).$$
(1)
(Wigner adds in footnote that it was found ‘by. L.Szilard and the present author some years ago for another purpose.’) The Wigner-function $`W(x,p)`$ became the prototype of all subsequent trials to simulate a quantum state through a statistical distribution over the classical phase space. As it was already obvious to Wigner himself, the sign of $`W(x,p)`$ is indefinite. By now, it is widely accepted that the existence of domains where $`W`$ is negative should mean that the state $`\widehat{\rho }`$ is essentially quantum in a sense that it can not be simulated by a statistical distribution $`W(x,p)`$ over the classical phase space. In the contrary case, when $`\widehat{\rho }`$’s Wigner-function is positive, one can say that $`\widehat{\rho }`$ exhibits classicality in the above particular sense.

After half century , Wigner was certainly the first among the most influential to support the concept later called environmental decoherence . Wigner, according to his Review of the quantum mechanical measurement problem, was impressed  by the work of Zeh  claiming that a macroscopic body can actually not be a closed system of its microscopic degrees of freedom. Wigner adopts such a reality and emphasizes the need of a new equation for apparently non-isolated objects. He modifies the Schrödinger equation by adding a second term to the r.h.s.:

$$\frac{d\widehat{\rho }}{dt}=\frac{i}{\mathrm{}}[\widehat{H},\widehat{\rho }]+\mathrm{irreversible}\mathrm{term}.$$
(2)
The irreversible term models the inevitable interactions of the macroscopic body with its environment. Perhaps, this was the first closed equation with the explicite intention of modeling environmental decoherence.

Decoherence  is responsible for the emergent classicality in quantized systems. In particular, the negative valued domains of the Wigner-function $`W(x,p)`$ might become washed out and grow positive under the influence of a modified Schrödinger equation like Eq. (2). Kiefer and the present author have recently proved an even stronger theorem . The positivity of $`W(x,p)`$ is achieved in finite time, contrary to our intuition that continuous variables would only asymptotically decohere.

Sec. 2 outlines the function and the equation proposed by Wigner in 1932 and 1983, respectively. In Sec. 3 the theorem of exact decoherence is recapitulated.

## 2 His Function (1932) and Equation (1983)

According to Wigner , the quantum state $`\widehat{\rho }`$ can be mapped into a normalized phase space distribution:

$$W(x,p)=\frac{1}{2\pi }xr/2|\widehat{\rho }|x+r/2\mathrm{e}^{ipr/\mathrm{}}dr.$$
(3)
There is an equivalent form of the above map, see e.g. :

$$W(x,p)=\mathrm{tr}\left[\widehat{\rho }𝒮\delta (x\widehat{x})\delta (p\widehat{p})\right]$$
(4)
where $`𝒮`$ stands for total symmetrization defined by successive application of the rules

$$𝒮\widehat{x}\widehat{f}=\frac{1}{2}\left(\widehat{x}\widehat{f}+\widehat{f}\widehat{x}\right),𝒮\widehat{p}\widehat{f}=\frac{1}{2}\left(\widehat{p}\widehat{f}+\widehat{f}\widehat{p}\right),$$
(5)
to the already $`𝒮`$-ordered function $`\widehat{f}=f(\widehat{x},\widehat{p})`$. The form (4) shows directly that the Wigner-function has the correct quantum mechanical marginal distributions for both canonical variables separately:

$`{\displaystyle W(x,p)dp}`$ $`=`$ $`\mathrm{tr}[\widehat{\rho }\delta (x\widehat{x})],`$ (6)
$`{\displaystyle W(x,p)dx}`$ $`=`$ $`\mathrm{tr}[\widehat{\rho }\delta (p\widehat{p})].`$ (7)
The form (4) of the Wigner-function is covariant against linear canonical transformations:

$$\left(\genfrac{}{}{0pt}{}{\widehat{x}}{\widehat{p}}\right)\left(\genfrac{}{}{0pt}{}{\widehat{x}^{}}{\widehat{p}^{}}\right)W(x,p)W^{}(x^{},p^{})=W(x,p).$$
(8)
Wigner proposes the following modification of Schrödinger’s equation in case of a macroscopic body :

$$\frac{d\widehat{\rho }}{dt}=\frac{i}{\mathrm{}}[\widehat{H},\widehat{\rho }]\underset{\mathrm{}m}{}\epsilon _{\mathrm{}}[\widehat{L}_\mathrm{}m,[\widehat{L}_\mathrm{}m,\widehat{\rho }]],$$
(9)
where $`\widehat{H}`$ is the total Hamiltonian and $`\widehat{L}_\mathrm{}m`$ are the multipole operators of the angular momentum. The second term is intented to decohere macroscopically different multipole moments. The strengths of their decoherence are being controlled by the rotation invariant parameters $`\epsilon _{\mathrm{}}`$. It turns out, however, that more typical is the decoherence between macroscopically different center of mass positions $`\widehat{x}`$. The modified Schrödinger equation retains the mathematical structure of Wigner’s one (9). Indeed, for a free object of mass $`m`$ it is this simple :

$$\frac{d\widehat{\rho }}{dt}=\frac{\mathrm{i}}{\mathrm{}}[\frac{\widehat{p}^2}{2m},\widehat{\rho }]\frac{D}{2\mathrm{}^2}[\widehat{x},[\widehat{x},\widehat{\rho }]].$$
(10)
This equation simplifies the environment as if it corresponded to a random external force $`F(t)=\sqrt{D}w(t)`$ where $`w(t)`$ is standard white-noise. The solution of the standard Schrödinger equation with the randomly fluctuating potential $`F(t)\widehat{x}`$, averaged over the noise $`w`$, yields the above modified Schrödinger equation .

The strength of decoherence is governed by a single parameter $`D`$. Assume a certain coherent width $`\sigma `$ for the initial quantum state. Dimensional analysis of the modified Schrödinger equation shows that the Hamiltonian and the irreversible terms have opposite tendencies. The unitary term increases $`\sigma `$ at a ‘coherent broadening’ characteristic time $`m\sigma ^2\mathrm{}`$. The irreversible term is tending to decrease $`\sigma `$ at a time scale $`\mathrm{}^2/D\sigma ^2`$ of ‘incoherent localization’. The two contrary effects may have a balance at the ‘stationary coherence width’

$$\sigma _0=\left(\frac{\mathrm{}^3}{Dm}\right)^{1/4},$$
(11)
achieved typically after a characteristic ‘decoherence time’

$$t_0=\sqrt{\mathrm{}m/D}.$$
(12)
Then the physical picture suggests that the quantum state has become a random (incoherent) mixture of wave packets whose characteristic width is the stationary value (11). The corresponding Wigner-function must obviously be positive since each contributing wave packet is assumed to have a positive Wigner-function.

## 3 Exact decoherence

Let us consider the evolution of the quantum state $`\widehat{\rho }(t)`$ under the influence of environmental decoherence as described by Eq. (10). The following theorem holds . Independently of the initial state $`\widehat{\rho }(0)`$, the solution $`\widehat{\rho }(t)`$ will exhibit exact classicality for $`tt_D`$ in a sense that the corresponding Wigner-function becomes non-negative:

$$W(x,p;t)0\mathrm{iff}tt_D.$$
(13)
The exact value of the decoherence time is calculable: $`t_D=3^{1/4}t_0`$.

The proof is the following. Substituting Eq. (3) into Eq. (10), we obtain the standard classical Fokker-Planck equation for the Wigner-function (3):

$$\frac{dW}{dt}=\frac{p}{m}\frac{W}{x}+\frac{D}{2}\frac{^2W}{p^2}.$$
(14)
Consider the Gaussian of width $`\sqrt{C}`$ in both $`x`$ and $`p`$:

$$g(x,p;C)=\frac{1}{2\pi C}\mathrm{exp}\left[\frac{x^2+p^2}{2C}\right],$$
(15)
and use it to coarse-grain an arbitrarily chosen Wigner-function. The result is always non-negative if and only if the coarse-graining scale is greater than one-half:

$$g(x,p;C)W(x,p)0\mathrm{iff}C1/2,$$
(16)
where $``$ denotes convolution. Indeed, for $`C=1/2`$ the coarse-graining yields the positive Husimi-function ; greater values of $`C`$ yield further coarse-graining which preserves the positivity of the Husimi-function. Now we can generalize the lemma (16). Let us generalize the Gaussian profile (15) first:

$$g(x,p;𝐂)=\frac{1}{2\pi }|𝐂|^{\frac{1}{2}}\mathrm{exp}\left[(x,p)\frac{1}{2𝐂}\left(\genfrac{}{}{0pt}{}{x}{p}\right)\right].$$
(17)
The following lemma holds, see also Ref.  by Khalfin and Tsirelson (1992):

$$g(x,p;𝐂)W(x,p)0iff|𝐂|1/4.$$
(18)
The correlation matrix $`𝐂`$ can always be transformed into the form $`C𝐈`$ with $`C=\sqrt{|𝐂|}`$. Hence, due to the covariance (8) of the Wigner-function, the lemma (18) follows from the special case (16).

Coming back to the Fokker-Planck equation (14), its solution can be written as the following time-dependent Gaussian coarse-graining:

$$W(x,p;t)=g(x,p;𝐂_W(t))W(xpt/m,p;0).$$
(19)
The correlation matrix of the coarse-graining profile is time-dependent:

$$𝐂_W(t)=Dt\left(\begin{array}{cc}t^2/3m^2& t/2m\\ t/2m& 1\end{array}\right),$$
(20)
as it can be inspected if we insert it into Eq. (19) which we substitute into (14). The determinant yields:

$$|𝐂_W(t)|=\frac{D^2t^4}{12m^2}.$$
(21)
The proof of the theorem (13) culminates in the observation that the determinant is monotone function, i.e., the Wigner-function $`W(x,p;t)`$ results from progressive coarse-graining. Applying the lemma (18) to the determinant, we conclude that the Wigner-function $`W(x,p;t)`$ is non-negative if $`|𝐂_W(t)|1/4`$. This condition is equivalent with $`t3^{1/4}t_0t_D`$. This completes the exact proof of the theorem whose earlier versions can be found in Refs. .

## 4 Conclusions

I discussed two items among Wigner’s remarkable contributions to the foundations of quantum mechanics: the first phase-space quasi-distribution (1932) and the first explicite equation of decoherence (1983), respectively. The simplified version of his decoherence equation leads to a theorem when applied to his function $`W(x,p;t)`$. I recapitulated the exact proof, given recently by Claus Kiefer and myself, guaranteeing the positivity of Wigner’s function after a finite time of decoherence. The purpose of my talk was to demonstrate Wigner’s continued impact even on the second century of quantum mechanics.

## Acknowledgement(s)

My research was supported by the Hungarian OTKA Grant No. 032640.


text:
# Robustness and Enhancement of Neural Synchronization by Activity-Dependent Coupling

## Abstract

We study the synchronization of two model neurons coupled through a synapse having an activity-dependent strength. Our synapse follows the rules of Spike-Timing Dependent Plasticity (STDP). We show that this plasticity of the coupling between neurons produces enlarged frequency locking zones and results in synchronization that is more rapid and much more robust against noise than classical synchronization arising from connections with constant strength. We also present a simple discrete map model that demonstrates the generality of the phenomenon.

Synchronous activity among individual neurons or their ensembles is a robust phenomenon observed in many regions of the brain, in sensory systems and in other neural networks. With constant synaptic connections the regions of neural synchronization are quite narrow in parameter space and the origin of the observed robustness of synchronization is not clear. It is known that many neurons in the cortex, in the cerebellum and in other neural systems are coupled through excitatory synaptic connections whose strength can be altered through activity-dependent plasticity. Indeed, this plasticity is widely thought to underlie learning processes, and in itself constitutes a broadly interesting phenomenon. Here we discuss its role in the synchronization of neurons in a network.

There have been recent experimental advances in the understanding of such plasticity, and, in particular, of the critical dependence on timing in presynaptic and postsynaptic signaling. Two manifestations of this kind of synaptic plasticity are the Spike-Timing Dependent Plasticity (STDP) Markram et al. (1997); Bi and Poo (1998) seen in excitatory connections between neurons, and its inverse, observed, for example, in the connections between excitatory and inhibitory neurons in the electrosensory lobe of fish Bell et al. (1999). The connections between excitatory neurons through inhibitory interneurons are typical in sensory systems Laurent (1999); Rabinovich et al. (2001) and cerebral cortex Bear and Linden (2000). These also express synaptic plasticity Perez et al. (2001) and play an important role in the control and synchronization of neural ensembles in hippocampus.

We report here on the synchronization of two model neurons coupled through STDP or inverse STDP synapses. We demonstrate that such coupling leads to neural synchronization which is more rapid, more flexible and much more robust against noise than synchronization mediated by constant strength connections. (For reviews, see Glass (2001); Elson et al. (1998); Coombes and Bressloff (1999)). We also build a simple discrete map that illustrates the enhancement of synchronization by activity-dependent coupling. The map allows us to speculate about the general applicability of learning-enhanced synchronization.

We consider here the simplest neural network: two neurons with unidirectional, activity-dependent excitatory synaptic coupling. Each neuron is described by the Hodgkin-Huxley equations with standard Na, K, and ‘leak’ currents Traub and Miles (1991):

$`C{\displaystyle \frac{dV_i(t)}{dt}}`$ $`=`$ $`g_{Na}m_i(t)^3h_i(t)(V_i(t)E_{Na})`$ (1)
$``$ $`g_Kn_i(t)^4(V_i(t)E_K)g_L(V_i(t)E_L)`$
$``$ $`I_{syn}(t)+I_{stim},`$
where $`i=1,2`$.

Each of the activation and inactivation variables $`y_i(t)=\{n_i(t),m_i(t),h_i(t)\}`$ satisfies first-order kinetics

$$\frac{dy_i(t)}{dt}=\alpha _y(V_i(t))(1y_i(t))\beta _y(V_i(t))y_i(t).$$
(2)
The parameters in these equations are given in par .

Each neuron receives a constant input $`I_{stim}`$ forcing it to spike with a constant, $`I_{stim}`$–dependent frequency. The second neuron is synaptically driven by the first via an excitatory current dependent on the postsynaptic $`V_2(t)`$ and presynaptic $`V_1(t)`$ membrane voltages:

$$I_{syn}(t)=g(t)S(t)V_2(t).$$
(3)
$`S(t)`$ is the fraction of open synaptic channels. It satisfies first-order kinetics:

$$\frac{dS(t)}{dt}=\alpha (1S(t))H(V_1(t))\beta S(t),$$
(4)
with $`H(V_1(t))=(1+\mathrm{tanh}(10V_1(t)))/4`$.

The time dependent synaptic coupling strength $`g(t)`$ is conditioned by the dynamics of the pre- and postsynaptic neurons. We consider two types of activity-dependent couplings: (1) an excitatory synapse with STDP, and (2) an excitatory synapse with inverse STDP. Through STDP $`g(t)`$ changes by $`\mathrm{\Delta }g(t)`$ which is a function of the time difference $`\mathrm{\Delta }t=t_{post}t_{pre}`$ between the times of post- and presynaptic spikes. We use the additive update rule

$$\mathrm{\Delta }g(t)=G(\mathrm{\Delta }t)=A\text{ sgn}(\mathrm{\Delta }t)\mathrm{exp}(\gamma |\mathrm{\Delta }t|)$$
(5)
for STDP, and $`\mathrm{\Delta }g(t)=G(\mathrm{\Delta }t)`$ for inverse STDP. We used $`A=0.004\mu `$S and $`\gamma =0.15\text{ms}^1`$.

We studied the synchronization properties of this coupled system by setting the autonomous period of the postsynaptic neuron to 15 ms, then evaluating the actual period of its oscillation $`T_2`$ as a function of the imposed autonomous oscillation period $`T_1`$ of the presynaptic neuron. In Fig. 1 we show $`T_1/T_2`$ as a function of $`T_1`$ in two cases: (a) a synaptic coupling with constant strength $`0.008\mu S`$ and (b) a synaptic coupling with inverse STDP. In the later case the steady-state coupling strength depends on the ratio of neuronal frequencies (c). Its average over all $`T_1`$ values is $`0.002\mu S`$, which is much lower than the strength in the case of constant coupling.

In Fig. 1a we see the familiar ‘Devil’s Staircase’ associated with frequency locking domains of a driven nonlinear oscillator. Only frequency locking with ratios 1:1, 2:1, 3:1, and 4:1 leads to synchronization plateaus with significant width. In Fig. 1b we see that the synchronization domains are substantially broadened due to activity-dependent coupling, especially for $`T_1/T_2=1`$. Some synchronization plateaus exhibit multistability, which we confirmed by observing the associated hysteresis. These results show that even a weak, but adaptive connection with strength that is determined dynamically is able to greatly enhance and enrich synchronization.

We also studied the robustness of this enhanced synchronization in the presence of noise by adding zero mean, Gaussian, white noise to the membrane currents of each neuron. We examined the behaviour of the system with RMS noise amplitudes $`\sigma =0.01,\mathrm{\hspace{0.17em}0.05},\mathrm{\hspace{0.17em}0.1},\text{and }0.5`$ nA.

For $`\sigma =0.01`$ nA no phase-locking plateaus were destroyed. At $`\sigma =0.05`$ nA the 4:1 plateau became distorted. Larger $`\sigma `$ sequentially eliminated synchronization plateaus until only the 1:1 plateau remained. The 1:1 plateau was seen for all $`\sigma `$. In Fig. 2 we illustrate the effect of the noise on synchronization when $`\sigma =0.1`$ nA with (a) constant and (b) inverse STDP coupling. While in (a) most of the plateaus have disappeared, in (b) the 1:1, the 2:1 and even the 3:1 frequency locking regimes remained. In sharp distinction to classical synchronization, frequency locking through activity-dependent coupling is significantly more robust in the presence of noise.

To understand the mechanisms behind such a remarkable robustness we studied the diffusion of oscillation phase caused by noise. For $`\sigma =0.5`$ nA in Fig. 3a we show that in the case of 1:1 synchronization and coupling with constant strength $`0.008\mu S`$ noise-induced phase diffusion results in 2$`\pi `$ phase slips that destroy synchronized state. Quite contrary Fig. 3b shows that in the case of activity-dependent coupling phase slips are absent and the phase difference does not increase. In this particular case the strength of coupling varied around the mean of $`0.0064\mu S`$ with standard deviation of $`0.0026\mu S`$.

In Fig. 4 we plot the average rate of phase slips for different amplitudes of the noise. In line with the above observation we see that in the case of activity-dependent coupling (dashed line) phase slips are suppressed in a wide range of noise amplitudes. We argue here that this suppression of phase slips is the primary mechanism responsible for robustness of synchronization mediated by activity-dependent coupling. After the introduction of a discrete map model we will discuss this mechanism in more detail.

We also considered synchronization through an activity-dependent synapse in the interesting case when the presynaptic neuron produces bursts of spikes and the postsynaptic neuron spikes irregularly. We found that synchronization through an STDP synapse is very fast; even a few spikes are enough for the frequency locking to establish itself. Neurons in the same set up with constant coupling synchronize much more slowly and only if the strength of the connection is appropriate for the given ratio of their frequencies. Hence, activity dependent synapses allow adaptation ‘on the run,’ synching a postsynaptic neuron to the firing properties of its presynaptic partner.

To understand the above results in a general way we have constructed a discrete time map model of periodic generators with STDP-like coupling. This map accounts for the dependence of the coupling strength on the activity of generators. Take $`T_1^0`$ and $`T_2^0`$ as the autonomous periods of the first and second generators. As a result of unidirectional coupling, the period of the second generator will change by some amount $`\mathrm{\Delta }T`$ each time it receives a spike from the first generator. Assuming initial phases to be 0, the time of the $`n+1`$-st spike of the first generator and $`m+1`$-st spike of the second generator are taken to satisfy

$`t_{n+1}^{(1)}`$ $`=`$ $`t_n^{(1)}+T_1^0`$ (6a)
$`t_{m+1}^{(2)}`$ $`=`$ $`t_m^{(2)}+T_2^0\mathrm{\Delta }T_{m,n},`$ (6b)
where $`n`$ and $`m`$ are such that $`t_m^{(2)}t_n^{(1)}t_{m+1}^{(2)}`$. In general, $`\mathrm{\Delta }T_{m,n}`$ would be a function of $`T_1^0`$, $`T_2^0`$, $`t_n^{(1)}`$, $`t_m^{(2)}`$, and the coupling strength $`g_{m,n}`$. We argue that the two main variables here are $`t_n^{(1)}t_m^{(2)}`$, and $`g_{m,n}`$. In the simplest case $`\mathrm{\Delta }T_{m,n}`$ can be approximated by

$$\mathrm{\Delta }T_{m,n}=g_{m,n}F(t_n^{(1)}t_m^{(2)}).$$
(7)
where the function $`F(x)`$ is the analog of a phase response curve Winfree (1980) for our model. To obtain results quantitatively comparable with our neuronal model, we fit it by non-negative quadratic function that describes phase response of our model neurons: $`F(x)=835+63x9x^2`$ for $`0xT_2^0`$ and 0 otherwise. $`g_{m,n}`$ obeys the inverse STDP update rules:

$`g_{m+1,n}`$ $`=`$ $`g_{m,n}G(t_{m+1}^{(2)}t_n^{(1)})`$ (8a)
$`g_{m,n}`$ $`=`$ $`g_{m,n1}G(t_m^{(2)}t_n^{(1)}).`$ (8b)
In Fig. 5 we show the Arnol’d Tongues calculated for the map (6-8) in the cases of (a) constant and (b) inverse STDP coupling. As with the model neurons, we see that activity-dependent coupling greatly enlarges the zones of synchronization.

This discrete map can be further analyzed to find its fixed points corresponding to $`n:m`$ synchronization and to examine their stability. We present here only the case of 1:1 synchronization. Then $`m=n`$, and the system of equations (6-8) can be written in the following simple form:

$`\tau _{n+1}`$ $`=`$ $`\tau _n+T_1^0T_2^0+g_nF(\tau _n)`$ (9a)
$`g_{n+1}`$ $`=`$ $`g_nG(T_1^0\tau _{n+1})G(\tau _{n+1}),`$ (9b)
where $`\tau _n=t_n^{(1)}t_n^{(2)}`$. The fixed points of (9) are given by $`g_n^f=(T_2^0T_1^0)/F(\tau _n^f)`$ and $`\tau _n^f=T_1^0/2`$. Stability calculations show that for such $`F(\tau )`$ and $`G(\tau )`$ these fixed points are stable. The second fixed point illustrates that activity-dependent coupling introduces a new limitation on the relationship between the phases of two oscillators. It is this limitation that causes the suppression of phase slips under the influence of noise. Detailed analysis shows that in the course of noise-affected synchronization the strength of activity-dependent coupling ajusts dynamically to keep this phase relationship close to satisfaction and, hence, suppresses phase slips.

In conclusion, we have analyzed the effects of activity-dependent coupling on synchronization properties of coupled neurons. We showed that such coupling results in a substantial extension of the temporal synchronization zones, leads to more rapid synchronization and makes it much more robust against noise. The enlargement of synchronization zones means that with STDP-like learning rules the number of synchronized neurons in a large heterogeneous population must increase. In fact, this is an aspect of the popular idea due to Hebb Hebb (1949). It is supported by the results in Miltner et al. (1999); Fell et al. (2001) which indicate that the coherence of fast EEG activity in the gamma band increases in a process of associative learning.

Based on our discrete map model results, we argue that the particular details of the signal-generating devices (e.g. neurons) and their connections (e.g. synapses) are not essential and the obtained results have general applicability. In fact, we observed similar phenomena of robust and enhanced synchronization in computer simulations of other types of periodic generators (such as Van-der-Pol and $`\theta `$-oscillators) with STDP-like activity-dependent coupling.

###### Acknowledgements.

This work was partially supported by U.S. Department of Energy Grants No. DE-FG03-90ER14138 and No. DE-FG03-96ER14592, NSF Grant No. PHY0097134, Army Research Office Grant No. DAAD19-01-1-0026, Office of Naval Research Grant No. N00014-00-1-0181, and NIH Grant No. R01 NS40110-01A2. R. Huerta thanks MCyT(Spain) BFI2000-0157.


text:
# Interaction between electronic structure and strain in Bi nanolines on Si(001)

## 1 Introduction

Nanometer-scale electronic technologies require not only the formation of nanoscale devices, but also nanoscale interconnections between the devices. At present, nanoscale structures may be made via at least three different approaches, all of which are likely to used in any scheme: top-down lithographic methods, such as e-beam lithography and AFM lithography; bottom-up methods, which involve positioning of ex situ “prefabricated” structures such as conducting molecules, semiconductor nanowires and carbon nanotubes using, e.g., scanning probes or micro-fluidics; or by in situ self-assembly, as in semiconductor quantum dots, or silicide nanowires. We are pursuing a hybrid self-assembly-based fabrication route, whereby a Bi nanoline is used as a nanowire or as a template for nanowires of other metals. These nanolines are very long—over 400 nm in some cases—and straight: a kink in a nanoline has never been seen. Unlike silicide nanowires, their width is constant, occupying the space of 4 dimers (1.5 nm) in the Si(001) surface. However, variable-bias STM of the nanolines indicates that they have a band gap which is larger than the surrounding surface, and suggests that they are not conductive.

Recently, we used a combination of atomistic structure calculations and experiments to identify the structure of the nanoline (which we have called the Haiku structure). An example of the Bi nanoline on a clean Si(001) surface is shown in Fig. 1, with the Haiku structure shown in the inset. It has a complex reconstruction in the silicon substrate reaching down five layers below the surface, though the surrounding surface is not reconstructed. The depth of the reconstruction suggests that the nanoline may be coupled strongly to its surroundings. Knowledge of the nanoline structure allows us for the first time to make detailed calculations of the electronic structure of the nanoline, which bear out the observations made in STM. The calculations also provide an explanation for the additional features that we have observed in recent atomic-resolution variable-bias STM of the Bi nanoline on the clean Si(001) surface at room temperature.

In this letter, we present high resolution, variable bias images of a Bi nanoline, and DFT calculations which explore the electronic structure of the Bi nanoline. In the next section, we discuss details of the experimental and theoretical techniques used. We then present the atomic-resolution images of the Bi-exposed Si(001) surface and interpret them in light of the various calculations.

## 2 Experimental and Theoretical Methods

The Si(001) substrate was cleaned using a standard process before being transferred into vacuum. The Si surface was prepared by flashing repeatedly to 1100C for a few seconds, until there was only a small pressure rise. The clean surface was checked with STM before Bi deposition began. Bi was evaporated from an effusion cell, a typical dose being Bi at 470C for 10 mins. STM images were taken at the deposition temperature between 570-600C, and at room temperature, using a JEOL 4500 XT UHV STM. The standard recipe for the formation of the Bi nanolines is to anneal the Si(001) surface close to 600C, under a Bi flux of around 1 ML/min.. In this case, the surface was annealed at a slightly lower temperature, around 570C, and quenched to 300C as soon as the presence of nanolines had been confirmed. In this way, we hoped to capture Bi nanolines while still growing, and hence make some observations which would lead to a likely nucleation mechanism. To minimise surface contamination, the sample was held at 300C while the Bi cell was cooled and the chamber pressure returned to its base level, before the sample was cooled to room temperature.

The density functional theory (DFT) calculations were performed in the Generalised Gradient Approximation (GGA) using the VASP code, with ultrasoft pseudopotentials, a plane wave cutoff of 150 eV (sufficient for energy difference convergence) and a Monkhurst-Pack k-point mesh with 4$`\times `$4$`\times `$1 points. The unit cell used contained a single Bi nanoline, and had ten layers of Si, with twenty atoms in each layer (forming a single dimer row ten dimers long with the p($`2\times 2`$) reconstruction) with the bottom two layers constrained to remain fixed and dangling bonds terminated in hydrogen. This unit cell was of sufficient size to for convergence of energy with cell depth, and long enough that the end Si dimers were representative of the clean, undistorted surface.

## 3 Electronic structure of Bi nanoline

Figure 1 shows a typical image from our STM experiments. Note that the Si(001) surface is clean, and that the image was taken at room temperature (in contrast to all previous experiments, where either the substrate was hydrogenated or the temperature was $`570^{}`$C). There are two key features to note: first, the Bi nanoline, which is the double line extending diagonally across the image from top left to bottom right; second, the features (marked ‘A’ in Fig.1) which appear as a double dot, and which are scattered across the surface; note that these features tend to align diagonally with each other, even forming a local c$`(4\times 4)`$ pattern in places (marked ‘B’). In this paper, we will concentrate on the electronic structure of the nanoline. We are actively investigating a possible structure for the double-dot features, which may be a precursor to the nanoline, and we shall present our results in future work.

In order to probe more closely the physical and electronic structure of the nanoline, a series of atomic-resolution variable bias STM images (in the area shown by the dotted box in Fig. 1) were taken and are shown in Fig. 2. The STM bias voltages are given in the figure caption. Filled-states images range from -2.0 V down to -0.3 V. An empty-states image at a bias voltage of +1.5 V is also presented. As has been previously reported, the Bi nanoline appears bright at large bias voltages of either sign, but as the bias voltage is reduced, the contrast of the nanoline reduces so that at around -0.8 V, the relative contrast of the nanoline and the Si substrate is the same, and at lower bias voltages, the nanoline becomes dark relative to the clean Si surface. By contrast, images of the Bi nanoline on the H:Si(001) surface always show the nanoline as a bright feature.

At very low biases, as in Fig. 2(d), substrate dimers either side of the nanoline become enhanced, i.e. appear to be brighter than the rest of the Si dimers, out to a distance of at least two dimers. (This is shown schematically in (d), by the lighter and darker rectangles.) This phenomenon of enhancement of neighbouring dimers at low bias voltages has been seen previously for single missing dimer defects (1DV) on Si(001). In that case, the enhancement was explained by the distortion of these dimers away from their normal structure by the strain field of the 1DV, resulting in a local change of the electronic structure, and the top-most occupied states being moved higher in energy. The tensile strain around the Bi nanoline is expected to be the cause of the enhancement seen here; this enhancement is modelled below. In Fig. 2(f) (which shows extremely high resolution), tensile strain near the Bi nanoline is seen directly: the darkening between the first and second Si dimers is more pronounced than between the second and third Si dimers, indicating an increase in separation between these dimers. In our DFT modelling, the dimer-dimer spacing close to the nanoline is increased by 3%, a significant strain.

Another feature visible in the images is a change of contrast between Bi dimers in the nanoline. While at elevated temperatures the Bi nanoline always appears uniform, at room temperature, higher resolutions may be achieved, and some subtle details resolved. Some dimers in the Bi nanoline in Fig. 2(a) appear different to the rest of the line (they seem to have a dark boundary around them); these separated dimers correspond to brighter dimers on the otherwise dark nanoline in the lower bias images (especially Fig 2(c)) and darker patches (marked as ‘A’) in the empty states image shown in Fig. 2(e). The reason for these features have not been determined. However, they are not seen in elevated-temperature images of the clean Bi nanoline, even at the same bias voltages, while similar dark patches were observed during experiments with adsorbed hydrogen and oxygen. It is therefore quite possible that they are the result of contamination by ambient water or hydrogen while cooling down to room temperature.

We have calculated the electronic structure of the Haiku model for the Bi nanoline using DFT, as described above. In order to understand the voltage contrast seen in STM, we have projected out the charge density associated with various states within certain energies of the Fermi level. We show the partial charge densities associated with all states within 0.2 eV, 0.6 eV and 1.0 eV of the Fermi energy in Figure 3 along with the complete charge density, for three systems: (a) the clean Si(001) surface; (b) the Haiku model of the Bi nanoline; and (c) the Haiku model with a hydrogenated Si(001) surface. For the last, the first image is missing as the first states occur 0.55 eV below the Fermi level (a side-on view of the Haiku structure is shown in its place). The partial densities are shown in a plane passing through a Si dimer atom. As the results are for a statically buckled Si(001) surface, alternate dimers are “up” and “down” (except on the hydrogenated surface, where the buckling does not occur). Only the up atoms appear bright, resulting in an appearance of alternate missing dimers; the difference is most clearly seen by comparison of the electron charge densities with the structural model of the Haiku shown in the top right-hand corner of Fig. 3.

The partial densities are important for the comparison with experiment, since the states close to the Fermi level will contribute strongly to the STM images at low biases. Localisation of a state near the Fermi level in the area around particular atoms or bonds (in comparison to bulk or clean Si) is indicative of a local strain, because strain forces structures away from their equilibrium state, and hence minimum energy configuration. We first identified this connection between strain and STM contrast (or states moved higher in energy) in the case of the single missing dimer defect on Si(001), and we will use it here to elucidate the substrate distortions around the Haiku structure. In the case of the clean Si(001) surface (first column of Fig. 3), the states near the Fermi level (a,i) are mainly localised on the surface dimers, and the density fades away beneath the surface layer. This is as expected, since the dimerisation of the Si(001) surface is associated both with dangling bonds (at the surface) and with strain (sub-surface). By contrast, the states near the Fermi level for the Bi nanoline on clean Si(001) (b,i) are reduced in density below the Si dimers, and are entirely excluded from the Haiku structure and substructure. Instead, the two Si dimers immediately to the side of the Haiku show a greatly enhanced charge density in the states within 0.2 eV of the Fermi level, both by comparison with their neighbours (the left-most dimers in (b,i)) and with the clean Si(001) surface (a,i), which is seen in STM as the brightening of the dimers either side of the Haiku at very low biases (Fig. 2(d) in particular). The localisation of these low-energy states indicates that the substructure of the Haiku, far from being a strained structure is more relaxed than the layers beneath a clean Si(001) surface. Furthermore, it demonstrates that the strain induced by the Bi nanoline is concentrated on the Si dimers either side of the Haiku structure, and on the second-to-third layer bonds beneath the dimers.

The electronic structure of the Bi nanoline with a H-terminated Si(001) surface is shown in Fig. 3(c). The major effect is the removal of the surface states associated with the Si dimer $`\pi `$-bonds, and the removal of the buckling on the surface. The Si atoms immediately beneath the Bi nanoline now show an increased charge density relative to the neighbouring Si atoms, which is clearly visible in the 0.6 eV and 1.0 eV images (c,ii & iii), indicating that the region below the Bi nanoline is more strained than bulk Si. Furthermore, enhancement can be seen in layers as deep as the sixth layer, indicating the deep strain caused by the presence of the Bi nanoline. This difference in the electronic states suggests that the contrast seen with the clean surface is due to the interaction of the buckled dimers and the nanoline (and in particular the “up” atom of a buckled dimer) which is removed with a hydrogenated surface.

The lack of any states on the Bi near the Fermi level is clearly the cause of the darkening of the nanoline relative to the Si(001) surface seen in STM at low voltages. Indeed, the first states seen on the nanoline lie over 0.5 eV away from the Fermi level (which agrees qualitatively with the bias voltage of 0.8 V where parity of appearance between the nanoline and the substrate occurs). In the full charge densities(Fig. 3(iv)), the density of states associated with the Bi dimers, combined with their higher physical height, explains their relative brightness in STM at higher bias voltages. In the case of the hydrogenated surface, the Bi dimers have a similar charge density as in the clean surface, but in this case, the Si $`\pi `$-bonds have been eliminated, so the Bi dimer remains bright in STM at all biases. The large density associated with the nanoline in the 0.6 eV(c,ii) and 1.0 eV(c,iii) images underline the mechanism behind the brightening of the Bi nanoline.

## 4 Conclusions

We have presented high resolution room temperature STM images of Bi nanolines on a clean Si(001) surface for the first time. These images allowed us to probe the electronic structure of the nanoline and the surrounding Si(001) surface, and, by extrapolation, the localised strain associated with the Haiku. We found that the two Si dimers neighbouring the nanoline show enhancement at low bias voltages, and we confirmed that the nanoline is darker than the surrounding surface for biases of less than 0.8 V. We have also presented first principles electronic structure calculations of the Haiku structure with both a clean and a hydrogenated Si(001) surface, and we have explained the voltage contrast of the nanolines using these calculations.

The localisation of the states near the Fermi level seen in the projected charge densities suggests that the Si atoms immediately surrounding the Haiku structure are strained (in part at least due to the interaction with the buckled Si dimers), while the Si substructure immediately beneath the Haiku itself is somewhat relaxed (though the hydrogenated surface plots suggest that it is slightly strained relative to perfect bulk Si). This goes some way to explaining the stability of the subsurface 5-7-5 ring structures, which serve as a highly effective relief mechanism for the epitaxial stress exerted by adsorbed Bi. The absence of states close to the Fermi level suggests that the nanolines are likely to block surface conduction perpendicular to the nanoline, and are not likely to act as a nanowire.

The Bi nanolines are clearly not suitable for conduction on their own. However, as we have shown before, hydrogen will adsorb preferentially on the Si(001) surface and not on the Bi, giving a natural, automatic masking technique. This should allow us to adsorb metals on the nanolines, creating exceptionally high quality nanowires. This research is under way, and will be presented in future work.

DRB thanks the Royal Society for funding through a University Research Fellowship. Calculations were performed at the HiPerSPACE Centre at UCL (JREI grant JR98UCGI). JHGO was supported by the Japanese Science and Technology Agency (JST) as an STA Fellow. This study was performed through Special Coordination Funds of the Ministry of Education, Culture, Sports, Science and Technology of the Japanese Government (Research Project on active atom-wire interconnects). We would like to thank Bill McMahon for bringing to our attention many of the examples of 5-7-5 structures, and for sharing unpublished STM data.


text:
# Equilibrium and Kinetics: Water Confined in Carbon Nanotube as 1D Lattice Gas

## Abstract

A simple 1D lattice gas model is presented, which very well describes the equilibrium and kinetic behaviors of water confined in a thin carbon nanotube found in an atomistic molecular dynamics(MD) simulation \[ Nature 414, 188 (2001) \]. The model parameters are corresponding to various physical interactions and can be calculated or estimated in statistic mechanics. The roles of every interaction in the water filling, emptying and transporting processes are clearly understood. Our results indicate that the physical picture of the single-file kinetics is very simple.

Transport of molecules through pores on nanometer scale is important for biology. Carbon nanotubes have been widely studied as pores filled with gases  and liquids . Recently, using molecular dynamics(MD) simulation, Hummer and coworkers  studied the filling and transporting of water molecules in a very thin carbon nanotube which is a $`(6,6)`$ tube with two open ends, $`4`$ Å radius and $`5.5`$ unit cells (about $`13.5`$ Å in length). The interaction between carbon atoms and water molecules was expressed by Lennard-Jones(LJ) potential. They found that the LJ parameters ($`ϵ`$ and $`\sigma `$) can sensitively affect the filling of water molecules. In their work, two group LJ parameters were used: (1) unmodified parameters: $`ϵ=0.114`$ kcal/mol, $`\sigma =3.2752`$ Å; (2) modified parameters: $`ϵ=0.065`$ kcal/mol, $`\sigma =3.4138`$ Å. For the former, it was found that the water molecules completely filled the nanotube (the average occupancy number $`n=5`$, the corresponding density is about $`5`$ times of that of bulk water, where the radius of the nanotube pore is about $`0.8`$ Å). For the latter, water molecules fluctuated between the empty ($`n=0`$) and filled ($`n=5`$) states, and it was a two-state like behavior, the lifetimes of the two states follow exponential distributions. The high density filling seems be surprising in the unmodified system due to the hydrophobic property of carbon nanotube, though the sharper distribution of the water molecule binding energy inside the nanotube was found in the MD simulation . The sensitive effect of LJ potential in the filling of water, the two-state like kinetics, and other remarkable MD results need to more clearly be understood in physics.

Here, based on a strict statistic mechanics treatment, we present a one-dimensional(1D) lattice gas(LG) model of water molecules confined in a thin tube with a finite length and two open ends. All parameters of this model have explicit physical sources, can be calculated from interactions. They are: (1) the local excess chemical potential $`\mathrm{\Delta }\mu `$ of water in the tube; (2) the free energy contribution $`_c`$ from LJ potential of carbon atoms; (3) the free energy $`_H`$ from hydrogen bonds between water molecules inside the tube (the rotational free degree of the molecules can be included); and (4) the additive free energy $`_a`$ of the water molecules located at ends of the tube, which come from the interaction with the outside water molecules. In this model, LJ potential from the nanotube wall only affect the parameter $`_c`$. The fitted parameters obtained from MD results  are in agreement with the directly calculated or estimated values from the corresponding interactions, that is, the obtained equilibrium properties are well consistent with the results of MD simulation for both the unmodified and modified systems. In the LG model, by analyzing transition states in which one or two water molecules partially occupied the end(s) of the tube, the whole free energy landscape, and the main reaction paths can be obtained, therefore, the two-state like kinetics of filling and emptying of water can be clearly explained.

Considering a nanotube filled with water molecules, we first partition the tube into $`N`$ cells along its long axis, the length of every cell is corresponding to the length of hydrogen bond between water molecules. We form our model step by step: (1) we consider one water molecule moves into a cell from the bulk water, where the cell is only a box without any interaction except limiting the water molecule in volume $`V_0`$. In this process, the changed free energy is $`\mathrm{\Delta }\mu =\mu _w+\mu _g`$, where $`\mu _w`$ is the chemical potential of bulk water, $`\mu _g=t\text{ }\mathrm{ln}(\lambda ^3/V_0)`$, is the chemical potential of ideal gas with density $`1/V_0`$, $`t`$ is a normalized temperature. In this paper, the energy unit is set as the room temperature $`300`$ K, if no explicit unit is written. $`\lambda =(2\pi \mathrm{}^2/mk_BT)^{1/2}0.23`$ Å is the de Broglie thermal wavelength of water. Actually, $`\mathrm{\Delta }\mu `$ is nothing but $`\mu _{nt}^{ex}`$, where $`\mu _{nt}^{ex}`$ is the local excess chemical potential of water inside the nanotube, about $`6.87`$ kcal/mol $`11.6`$ from MD simulation . (2) if there is interaction $`V_e(𝐫)`$ between the wall of the cell and water molecule, except $`\mathrm{\Delta }\mu `$, an additive free energy is

$`_c=t\text{ }\mathrm{ln}({\displaystyle \frac{1}{V_0}}{\displaystyle _{V_0}}𝑑𝐫\text{ }\mathrm{exp}[V_e(𝐫)/t]).`$ (1)
For nanotubes, $`V_e(𝐫)`$ is the Lennard-Jones(LJ) potential between the water molecule and all carbon atom of the nanotube. Calculation shows that $`V_e(𝐫)`$ is almost independent on the azimuth $`\varphi `$ and axis position $`z`$ in infinite nanotubes , For finite nanotubes, the result is also right, except a $`z`$ dependence near the ends of the nanotubes. In $`(6,6)`$ nanotube, due to its very small radius, for both the unmodified and modified LJ parameters, $`V_e(r)`$ can be well approximated as a box potential well, whose widths (about $`0.8`$ Å) are in agreement with the interior radius of the tube found in MD , and the depths of well are about $`89`$ and $`56`$, respectively. It is not very difficult to numerically calculate the more exact value of $`_c`$, but it is not necessary; (3) if two water molecules locate at two near adjacent cells, a hydrogen bond will be formed, the interaction is denoted $`V_H(𝐫)`$. Under the box potential well approximation of $`V_e(r)`$, the free energy can be simply written as,

$`_H=t\text{ }\mathrm{ln}{\displaystyle \frac{1}{V_0^2}}{\displaystyle _{V_1}}𝑑𝐫_1{\displaystyle _{V_2}}𝑑𝐫_2\text{ }\mathrm{exp}[V_H(𝐫_1𝐫_2)/t],`$ (2)
where, due to the limit of the length of hydrogen bond and the nanotube wall, the two water molecules will occupy two cells $`V_1`$ and $`V_2`$ (both volumes are $`V_0`$), respectively. The freedom degree $`𝐫_2`$ is limited, only some rotation can be allowed, we have,

$`_H=E_Ht\text{ }\mathrm{ln}({\displaystyle \frac{a\lambda ^3}{V_0}})=E_H\mu _g\mu _i,`$ (3)
where $`E_H`$ is the energy of the hydrogen bond. All other contributions, including the rotation of the non-spherical water molecules , the rotation of water molecules around their hydrogen bond, and the possible angle changes of the hydrogen bonds, etc., are writen as the factor $`a`$ ($`a>1`$ indicates the existence of these freedom degrees), the corresponding free energy is noted as $`\mu _i`$. While the number ($`n`$) of water molecules formed a continuous hydrogen-bond chain inside the tube is not very large, we suppose that $`\mu _i`$ is independent on $`n`$; (4) if a water molecule locate at ends of the nanotube, besides $`_c`$ (while water is partially out of the nanotube, even $`_c`$ also increase), there is an additive free energy $`_a`$ from the interaction between the water molecule and the outside water. In fact, $`_a`$ is also from hydrogen bonds of water molecules, but due to the lower density of bulk water, it is higher than $`_H`$.

Based on these considerations, the water molecules inside the nanotube can be treated as a simply 1D LG model with $`N`$ sites, including three parameters. Only $`_c`$ is dependent on LJ parameters, and the dependence can be approximately calculated. Hamiltonian of the model is,

$`={\displaystyle \underset{i=1}{\overset{N}{}}}_i\text{ }c_i+_H{\displaystyle \underset{i=1}{\overset{N1}{}}}c_ic_{i+1},`$ (4)
where $`_i=\mathrm{\Delta }\mu +_c+(\delta _{i,1}+\delta _{i,N})\text{ }_a`$, $`\delta _{i,j}`$ is Kronecker delta, $`c_i=1`$(occupied), or $`0`$(unoccupied). The model is equivalent to an Ising model, if we transform $`s_i=2c_i1`$. For a 1D Ising model, if the external field $`H=0`$, the states with $`M`$ and $`M`$ are symmetry, and the energy barrier is very low, where $`M=\frac{{\scriptscriptstyle s_i}}{N}`$. In the LG model, the corresponding external field $`H=(_i+_H)/2`$. If $`_c`$ is very low, $`H>0`$, water molecules can fill the nanotube, while $`_c`$ increases (the well depth decreases), $`H`$ decreases. If $`H0`$, there is an equilibrium between the filling and emptying states. However, some details must be considered as we try to explain the MD results of Ref.  by using the LG model: (1) the finite sites $`N`$; (2) the additive energy $`_a`$ at ends. Since $`N`$ is not very large, we can directly enumerate all $`2^N`$ states. For general case, Hodak and Girifalco  have written the partition function for the LG model with any finite sites and the additive energy.

For unmodified and modified LJ parameters, Hummer et al. have given the probability of finding exactly $`N`$ water molecules inside the nanotube (Figure 2 of Ref. ). In their MD simulation, $`N=5`$. We can obtain the model parameters ($`_H6.27`$, $`_a2.93`$, and $`_c4.96`$) by fitting the MD probability data of the modified system. The regenerated free energy profile is very well in agreement with the MD probability distribution in the modified system as shown in Fig.1. Here, the fitting $`_c`$ is well consistent with the estimated value from LJ potential by Eq.(1). The hydrogen bond of water inside the nanotube is about $`7`$ kcal/mol $`11.8`$ (Figure 3(a) of Ref. ), is stronger than their average value in bulk water ($`5`$ kcal/mol $`8.4`$). From Eq.(3), we know that the free energy contribution $`\mu _i`$ of all rotational freedom degrees is about $`0.7`$. The $`_a`$ is only half of $`_H`$, the physical reason is the lower density of bulk water. Since both $`_a`$ and $`_H`$ are independent of LJ potential, their values can be used in the unmodified system. For the large statistic errors of MD data in the unmodified system, replacing to fit MD results, we directly use that estimated $`_c8`$. The obtained free energy profile is shown in inset of Fig.1, it is approximately in agreement with the MD probability distribution. For equilibrium states, due to the large $`_H`$, water is favor to the states with a single sequentially occupied sites, for example, $`P(10100)`$ is only about $`1/500`$ of $`P(11000)`$, where $`P`$ is the probability of states. Similarly, water is favor to the states whose end sites are occupied. For example, $`P(01100)`$ is only about $`1/20`$ of $`P(11000)`$ due to $`_a`$. For a good approximation, we can calculate the partition function only considering the states with a single sequence of hydrogen-bond chain and at least occupied an end of the tube. We call it as single-sequence end occupied approximation (SSEA). In a very current MD simulation , for short nanotube($`13.5`$ Å) and long nanotube($`27`$ Å), Hummer and coworkers carried to a longer time MD simulation, by changing LJ parameters by varying a control variable $`\eta `$, $`ϵ=ϵ_1\eta ^2`$, and $`\sigma =\sigma _1/\eta ^{1/6}`$, where $`\eta =1.0`$ is corresponding to the unmodified LJ parameters. The model can very well explain these results. For example, (1) very rare non-SSEA events are found; (2) the approximated linear free energy profile $`G(n)`$ for $`1n<N`$ (Using the SSEA of the LG model, $`G(n)=n\mathrm{\Delta }G+const`$, where $`\mathrm{\Delta }G=\mathrm{\Delta }\mu +_c+_H`$.); (3) $`\mathrm{\Delta }G`$ is proportional to $`\eta `$ (Our calculations show that the LJ interaction well depth is proportional to $`\eta `$, therefore $`\mathrm{\Delta }G`$ is proportional to $`\eta `$.). Actually, for explaining the equilibrium property found in their MD, Hummer and coworkers have formed a free energy formula (Eq.(3) of Ref. ), which is equivalent to the SSEA of our model.

The equilibrium property of water filled into nanotube has been very well described in the LG model. Before investigating the kinetics, we discuss that the states with $`6`$ water molecules inside the nanotube, its probability is not very small in MD simulation . It may be corresponding to a state where the ends of the nanotube are partially occupied by two water molecules as shown in Fig.2. The additive free energy of a molecule partially occupied is denoted as $`_a^{}`$, which is different from $`_a`$, due to both differences of carbon-water interaction and the interaction with the outside water. Since water molecules must pass through these partially occupied states as their filling or emptying, the free energy may be comparable with the kinetic energy barrier. For the states which $`n<6`$, the effects of the partially occupied sites can be ignored, then we have $`_a^{}=_a+[G(6)G(5)_H\mathrm{\Delta }\mu _c]/2`$. We obtain that $`_a^{}`$ is $`0.8`$ and $`1.76`$ for unmodified and modified system, respectively. We write that $`_a^{}=(\alpha 1)_c+\delta _w`$, where $`\alpha `$ is a corrected factor of the carbon-water interaction while water molecule partially move out of the nanotube, and $`\delta _w`$ is the contribution from the outside water near the end of nanotube. If we suppose that both $`\alpha `$ and $`\delta _w`$ are independent of LJ parameters, we have $`\alpha 0.67`$ and $`\delta _w3.4`$. The obtained factor $`\alpha `$ is well in agreement with the calculated result of LJ potential while the water molecule is halfly out of the tube, and $`\delta _w`$ is comparable to $`_a`$. It maybe indicate a real physical picture.

The kinetics of filling and emptying is determined by the real reaction paths. The existence of two minimums of the free energy profile $`G(n)`$ is not a sufficient condition for the two-state like behavior of filling and emptying. The understanding of the whole free energy landscape is necessary. Generally, a rate equation can be written as, $`\frac{dP_a}{dt}=_bk(ba)P_b_bk(ab)P_a`$, where $`P_a`$ is the possibility of state $`a`$, $`k(ba)=\gamma \text{ }\mathrm{exp}[\beta (F_{a,b}F_b)]`$ is the rate from state $`b`$ to $`a`$, where $`\gamma `$ is a factor, $`\beta =1/k_BT`$, $`F_{a,b}`$ is the free energy barrier or saddle point (transition state) separating the state $`b`$ and state $`a`$ , $`F_b`$ is the free energy of the initial state $`b`$. Obviously, only two connected states which they have a common boundary can transit each other. For the filling, emptying or transporting of water, there are $`2^N`$ states (here we have supposed that the outside water molecules be always in equilibrium). The basic motion is that a water molecule jump a step inside nanotube or pass through an end of the tube. Since the barrier of broken hydrogen bonds is too high, the collective motion of a hydrogen-bond chain often occure , it means all molecules of the chains move a basic step at the same time. For any cases, the axis positions $`z`$ of the all moving molecules are good reaction coordinates. In the motion of passing an end of nanotube, the free energy barrier locates nearly at the end point of the tube (molecule half occupy the end cell), so we approximately have the barrier $`_b_a^{}`$, but for the motion of water inside of the tube, the barrier is very small ($`0`$) due to the independence of LJ potential on $`z`$. For the collective motion, the free energy barrier is the sum of the barrier of every moving molecule. Besides this, the facter $`\gamma `$ is different from that of single molecule motion. Actually, $`\gamma `$ is dependent on the thermal velocity $`v_T\sqrt{k_BT/m}1/\sqrt{n}`$, where $`n`$ is the number of molecules in the collective moving chain. While the length of a cell $`l`$ is very small, $`\gamma =v_T/l`$ , so $`\gamma 1.5/\sqrt{n}`$ ps-1.

Consequently, based on the simple estimation, we can draw the whole free energy landscape. Obviously, any non-SSEA state can fastly transit to the connected SSEA states, but the inverse transition is very slow, and the barrier between neighbor SSEA states is lower than that between the SSEA states and non-SSEA states. For example, according to our estimation, from $`(10000)`$ to $`(11000)`$, the overcome barrier is about $`\mathrm{\Delta }G+_a^{}_a1.4`$, but from $`(10000)`$ to an non-SSEA state $`(10001)`$, the overcome energy is about $`\mathrm{\Delta }\mu +_c+_a^{}4.8`$. Therefore, the reaction paths almost do not pass through the non-SSEA states in the filling, emptying and transporting processes. It indicates that SSEA describes the kinetics behavior very well. Based on their MD results that the probability of non-SSEA states is very low at equilibrium, Hummer and coworkers have formed a kinetic model only including the reaction between the neighbor SSEA states (Eq.(2) of Ref. ), and found that it worked very well. However, the reason of all reaction paths including non-SSEA states can be neglected is the high free energy barriers, not only the high free energy itself of the states.

Now, we can easily explain the two-state like behavior of filling and emptying: water molecules collectively moving a step into or out from the nanotube, the barriers of the middle states ($`1n<N`$) are small (about $`1.4`$ and $`1.1`$, respectively), and the effect of the pre-exponent factor is not obvious, so the lifetime of these states is very short, the filling and emptying processes have a two-state like behavior. Actually, we can calculate the rate constant of every reaction and the lifetime of every middle state. For example, we have that $`k(1000000000)0.5`$, the rate constant of the inverse process is about $`0.012`$, where unit is ps-1. The order of magnitude of the values is in agreement with the estimated results of MD (about $`0.3`$ and $`0.009`$ , respectively.).

As this work was being finished, we noticed a work , in which Maibaum and Chandler treated whole system consisting of the channel of nanotube and the outside bulk liquid as a lattice gas by a MC simulation. It was shown that there are similar behavior found in MD simulation in generic liquid states. Under some approximations, they also analytical analyzed the equilibrium property of the system. However, for quantitatively explaining the behavior of water confined in nanotube, we should consider the characters of water, for example, the average density of bulk water is only $`1/5`$ of that of the filled state in the tube, then the interaction ($`_a`$) between a water molecule at an end of the tube and the outside molecules is weaker than the interaction between two adjacent inside molecules. For the equilibrium property of system, the whole outside bulk water can be thought as an aqueous reservoir and contribute an additive free energy $`_a`$ to the inside molecules at ends of nanotube, but for the kinetics, we should identify these states, in which different number outside water molecules contact with an inside molecule one by one by hydrogen bonds. If we denoted the state in which $`m`$ outside molecules contact with $`n`$ inside molecules one by one as “$`(m)\mathrm{\_}[n]`$” , then from the state, after a filling or emptying step, the final state is “$`(m1)\mathrm{\_}[n+1]`$” or “$`(m+1)\mathrm{\_}[n1]`$”, respectively. The complete rate chain will consist of all states with any $`m`$ (from $`0`$ to $`\mathrm{}`$). However, since the longer outside molecule chain is more instable, in our interesting time scale, the longer chain will be thought as rapidly equilibrating with the surrounding water. Depending on the considered largest $`m`$, we have different order approximations of the complete rate equation. In this paper, we used the lowest order ($`0`$order) approximation in which we supposed that one outside molecule rapidly contacted with and departed from the end inside molecule, the equilibrium can rapidly reach in our interested time scale. If we consider higher order ($`1`$order) approximation, the details of energy barrier will be changed, but the differences between the two approximations is not very large, so we think that the principal character of the kinetics of the system has been captured by the simple 1D LG model. We wish to give more quantitative agreement on rate constants by using higher order approximation in next work.

X. Z. thanks Professor Hummer for sending him their MD data and a new paper of Ref.  and telling him the work of Maibaum and Chandler (Ref. ). X. Z also thanks Professor Chandler for comments on an earlier draft of this paper, and Mr. Z.-C. Tu for discussions. X. Z is financially supported by the Grants-in-Aid for Scientific Research of JSPS.

(Correspondent address: zhou@pe.titech.ac.jp).


text:
# Reply to ”Comment on One loop renormalization of soliton quantum mass corrections in (1+1)-dimensional scalar field theory models” (hep-th/0211149)

## Abstract

I agree with the authors of hep-th/0211149 that the claim made in Phys.Lett. B542, 282 (2002) is incorrect and that the derivation of its main formula, although correct, contains two compensating errors. In this reply the main formula of Phys.Lett. B542, 282 (2002) is rederived. This new derivation shows that not only the energy momentum cut off regularization method still works in the calculation of the soliton quantum mass corrections, but also that the so called mode number regularization emerges naturally from it.

In Ref.  the authors refuted the claim made in Ref. , that previous works on the one loop soliton quantum mass corrections did not include a surface term from a partial integration. I agree with the authors of Ref. . As the authors of this reference pointed out, in the derivation of the main formula of Ref.  I used incorrectly the asymptotic behaviour of the continuous phase shift and also there was a wrong counting of the density of states. These two mistakes compensate in such a way that the correct formula, Eq. (19) of Ref. , was obtained. The way I derived that equation, ”showed” that the final result is independent of the regularization method. On the other hand, the authors of Ref.  showed that the energy momentum cut off regularization method gave incorrect results when it was used in the calculation of the soliton quantum mas corrections. Therefore, in Ref.  I concluded that the claim above was wrong. In this reply I rederive the formula given by Eq. (19) of Ref. , using the energy momentum cut off regularization method, showing that it gives the correct result for the soliton quantum mass corrections. Then my position about the claim made in Ref.  remains to be the same.

In the notation of Ref. , the bare one loop soliton quantum mass correction is given by

$$\mathrm{\Delta }M_{bare}=\frac{1}{2}\underset{i}{}\omega _i+\frac{1}{2}\underset{q}{}\omega (q)\frac{1}{2}\underset{k}{}\omega ^0(k),$$
(1)
where $`i`$ and $`q`$ label the eigenfrequencies of respectively the discrete and continuous modes in the presence of the soliton, whereas $`k`$ labels those of the vacuum modes. Since $`\omega (q)=\omega (q)=\sqrt{q^2+m^2}`$ and $`\omega ^0(k)=\omega ^0(k)=\sqrt{k^2+m^2}`$ I can write Eq. (1) as

$$\mathrm{\Delta }M_{bare}=\frac{1}{2}\underset{i}{}\omega _i+\underset{q=0}{\overset{\mathrm{}}{}}\omega (q)\underset{k=0}{\overset{\mathrm{}}{}}\omega ^0(k).$$
(2)
The continuous sums in the above equation are quadratically divergent and generally the difference of these is logarithmically divergent. Then these expressions need to be regularized. As an alternative to the energy momentum cut off regularization method, the authors of Ref.  introduced the so called mode number regularization method. Below I will show that not only the energy momentum cut off regularization method gives correct results in the calculation of the soliton quantum mass corrections but also, that the so called mode number regularization method emerges naturally from the former. The idea is to use the same cut off $`\mathrm{\Lambda }`$ in the sums over $`q`$ and $`k`$ in Eq. (2), then

$$\mathrm{\Delta }M_{bare}=\frac{1}{2}\underset{i}{}\omega _i+\underset{q=0}{\overset{\mathrm{\Lambda }}{}}\omega (q)\underset{k=0}{\overset{\mathrm{\Lambda }}{}}\omega ^0(k).$$
(3)
To write the continuous sums in Eq. (3) in an integral form, I enclose the system in a box of finite size $`L`$, impose periodic boundary conditions and finally let $`L\mathrm{}`$. For the free modes I obtain $`k_n=2\pi n/L`$, $`n=0,\pm 1,\pm 2,\mathrm{}`$, from which the density of states for the free modes $`L/(2\pi )`$ is obtained. On the other hand I have

$$q_n=\frac{2\pi n}{L}\frac{\delta (q_n)}{L},n=0,\pm 1,\pm 2,...$$
(4)
The phase shift $`\delta (q)`$, can be chosen in two ways. The customary in scattering theory is to choose it (for physical reasons) in such a way that $`\delta (\pm \mathrm{})=0`$. In this case the phase shift is discontinuous (with discontinuity at the origin) and I denote it as $`\delta _D(q)`$. The other way is to choose it in such a way that $`\delta (0)=0`$ and in this case it is a continuous function and I denote it as $`\delta _C(q)`$. I expect that the result for Eq. (3) be independent of the phase shift being continuous or discontinuous. For large momentum, $`\delta _D(q)`$ behaves like

$$\delta _D(\mathrm{\Lambda })=\frac{V}{2\mathrm{\Lambda }},$$
(5)
where $`V`$ is given by Eq. (15) of Ref. . Also I have that $`\delta _D(0^+)=𝒩\pi `$, where $`𝒩`$ is the number of discrete eigenfrequencies. On the other hand for $`\delta _C(q)`$ I have

$$\delta _C(\mathrm{\Lambda })=𝒩\pi \frac{V}{2\mathrm{\Lambda }}.$$
(6)
The derivatives of the continuous and discrete phase shifts are the same, and I denote them in both cases as $`d\delta (k)/dk`$.

In terms of discretized eigenfrequencies, Eq. (3) can be written as

$$\mathrm{\Delta }M_{bare}=\frac{1}{2}\underset{i=1}{\overset{𝒩}{}}\omega _i+\underset{n=n_0}{\overset{N^{}}{}}\sqrt{q_n^2+m^2}\underset{n=0}{\overset{N}{}}\sqrt{k_n^2+m^2},$$
(7)
where in the second term the sum starts at $`n=n_0`$ not necessarily equal to zero and ends at $`n=N^{}`$ not necessarily equal to $`N`$, since I require that in both sums, $`q_n`$ and $`k_n`$ ranges, according to Eq. (3), from zero to the same cut off $`\mathrm{\Lambda }`$. For the free eigenfrequencies $`N`$ is given in terms of the cut off $`\mathrm{\Lambda }`$ by,

$$\mathrm{\Lambda }=2\pi N/L.$$
(8)
To find $`n_0`$ and $`N^{}`$ I consider separately the continuous and discontinuous phase shifts.

Discontinuous phase shift In this case, setting $`q_{n_0}=0^+`$ in Eq. (4) and using $`\delta _D(0^+)=𝒩\pi `$, I have $`2\pi n_0/L𝒩\pi /L=0`$, from which I obtain $`n_0=𝒩/2`$. On the other hand setting $`\mathrm{\Lambda }=q_N^{}`$ in Eq. (4) and using $`\delta _D(\mathrm{\Lambda })0`$ I obtain $`\mathrm{\Lambda }=2\pi N^{}/L`$ and comparing with Eq. (8) I get $`N^{}=N`$. I have to call attention to the fact that the solution $`n_0=𝒩/2`$, only makes sense when $`𝒩`$ is even. For odd $`𝒩`$ it is obtained a half integer value for $`n_0`$, a solution not present in the spectrum as given by Eq. (4). Of course such a solution does not satisfies the periodic boundary conditions ***It is easy to see that it satisfies anti periodic boundary conditions. Then, in this last case I have to drop out the term $`n_0=𝒩/2`$ from the sum in the second term of Eq. (7), the sum beginning at the next integer larger than $`𝒩/2`$. If I write $`𝒩=2+1`$, this means that the sum must begin in $`n=+1`$. Also in the last sum of Eq. (7) the contribution from the zero momentum, $`n=0`$, must be multiplied by a factor $`1/2`$, since in the original sum given by Eq. (1) the corresponding term is counted only once. To proceed with my derivation I consider first, the case $`𝒩=2`$ even. In this case $`n_0=`$ and Eq. (7) can be written as

$`\mathrm{\Delta }M_{bare}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _i+{\displaystyle \underset{n=}{\overset{N}{}}}\sqrt{q_n^2+m^2}{\displaystyle \underset{n=0}{\overset{N}{}}}\sqrt{k_n^2+m^2}`$ (9)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _i{\displaystyle \underset{n=0}{\overset{1}{}}}\sqrt{k_n^2+m^2}+{\displaystyle \underset{n=0}{\overset{N}{}}}\left(\sqrt{q_n^2+m^2}\sqrt{k_n^2+m^2}\right),`$ (10)
From Eq. (4) I obtain

$$\sqrt{q_n^2+m^2}=\sqrt{k_n^2+m^2}\frac{\delta _D(k_n)}{L\sqrt{k_n^2+m^2}}+𝒪\left(L^2\right).$$
(11)
Using Eq. (11), $`\delta _D(0^+)=2\pi `$ and Eq. (5) in Eq. (10), and taking $`L\mathrm{}`$ I get

$`\mathrm{\Delta }M_{bare}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _im{\displaystyle \frac{\omega (k)}{2\pi }}\delta _D(k)|_{0^+}^\mathrm{\Lambda }+{\displaystyle _0^\mathrm{\Lambda }}{\displaystyle \frac{dk}{2\pi }}\omega (k){\displaystyle \frac{d}{dk}}\delta (k)`$ (12)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _i+{\displaystyle \frac{V}{4\pi }}+{\displaystyle \frac{1}{2}}{\displaystyle _\mathrm{\Lambda }^\mathrm{\Lambda }}{\displaystyle \frac{dk}{2\pi }}\omega (k){\displaystyle \frac{d}{dk}}\delta (k)+𝒪\left(\mathrm{\Lambda }^1\right).`$ (13)
Then, subtracting the tadpole graph contribution, as given by Eq. (14) of Ref. ,

$$\frac{V}{4}_\mathrm{\Lambda }^\mathrm{\Lambda }\frac{dk}{2\pi }\frac{1}{\sqrt{k^2+m^2}}$$
(14)
from Eq. (13) and taking $`\mathrm{\Lambda }\mathrm{}`$ I obtain the main formula, Eq. (19) of Ref. .

In the case in which $`𝒩=2+1`$, as explained above, I have to multiply the $`n=0`$ contribution for the free eigenfrequencies with a factor $`1/2`$ and the sum over the eigenfrequencies in the presence of the soliton must begin in $`n=+1`$. Then, Eq. (7) can be written as

$`\mathrm{\Delta }M_{bare}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _i+{\displaystyle \underset{n=+1}{\overset{N}{}}}\sqrt{k_n^2+m^2}{\displaystyle \underset{n=1}{\overset{N}{}}}\sqrt{q_n^2+m^2}{\displaystyle \frac{m}{2}}`$ (15)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _i{\displaystyle \underset{n=1}{\overset{}{}}}\sqrt{k_n^2+m^2}{\displaystyle \frac{m}{2}}+{\displaystyle \underset{n=1}{}}\left(\sqrt{k_n^2+m^2}\sqrt{q_n^2+m^2}\right).`$ (16)
Using Eq. (11), $`\delta _D(0^+)=(2+1)\pi `$ and Eq. (5) in Eq. (16) and taking $`L\mathrm{}`$ I obtain

$`\mathrm{\Delta }M_{bare}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _i\left(+{\displaystyle \frac{1}{2}}\right)m{\displaystyle \frac{\omega (k)}{2\pi }}\delta _D(k)|_0^\mathrm{\Lambda }+{\displaystyle _0^\mathrm{\Lambda }}{\displaystyle \frac{dk}{2\pi }}\omega (k){\displaystyle \frac{d}{dk}}\delta (k)`$ (17)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _i+{\displaystyle \frac{V}{4\pi }}+{\displaystyle \frac{1}{2}}{\displaystyle _\mathrm{\Lambda }^\mathrm{\Lambda }}{\displaystyle \frac{dk}{2\pi }}\omega (k){\displaystyle \frac{d}{dk}}\delta (k)+𝒪\left(\mathrm{\Lambda }^1\right).`$ (18)
Again, subtracting the tadpole graph contribution, given by expression (14), from Eq. (18) I recover the main formula of Ref. .

continuous phase shift In this case, since $`\delta _C(0)=0`$, from Eq. (4) I find $`n_0=0`$. On the other hand for the upper limit I have, from Eqs. (4) and (6)

$$\mathrm{\Lambda }=\frac{2N^{}\pi }{L}+\frac{𝒩\pi }{L}.$$
(19)
From Eqs. (8) and (19) I find that $`N^{}=N𝒩/2`$. In order to check that also in this case the main formula of Ref.  is obtained I consider for simplicity only the case in which $`𝒩=2`$. In this case Eq. (7) can be written as

$`\mathrm{\Delta }M_{bare}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _i+{\displaystyle \underset{n=0}{\overset{N}{}}}\sqrt{q_n^2+m^2}{\displaystyle \underset{n=0}{\overset{N}{}}}\sqrt{k_n^2+m^2}`$ (20)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _i{\displaystyle \underset{n=0}{\overset{1}{}}}\sqrt{q_n^2+m^2}+{\displaystyle \underset{n=0}{\overset{N}{}}}\left(\sqrt{q_n^2+m^2}\sqrt{k_n^2+m^2}\right).`$ (21)
From Eq. (4) I obtain

$$\sqrt{q_n^2+m^2}=\sqrt{k_n^2+m^2}\frac{k_n\left(\delta _C(k_n)+2\pi \right)}{L\sqrt{k_n^2+m^2}}+𝒪\left(L^2\right).$$
(22)
Using Eq. (22), $`\delta _C(0)=0`$ and Eq. (6) in Eq. (21), and taking $`L\mathrm{}`$ I obtain

$`\mathrm{\Delta }M_{bare}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _im{\displaystyle \frac{\omega (k)}{2\pi }}\left(\delta _C(k)+2\pi \right)|_0^\mathrm{\Lambda }+{\displaystyle \frac{1}{2}}{\displaystyle _\mathrm{\Lambda }^\mathrm{\Lambda }}{\displaystyle \frac{dk}{2\pi }}\omega (k){\displaystyle \frac{d}{dk}}\delta (k)`$ (23)
$`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{𝒩}{}}}\omega _i+{\displaystyle \frac{V}{4\pi }}+{\displaystyle \frac{1}{2}}{\displaystyle _\mathrm{\Lambda }^\mathrm{\Lambda }}{\displaystyle \frac{dk}{2\pi }}\omega (k){\displaystyle \frac{d}{dk}}\delta (k)+𝒪\left(\mathrm{\Lambda }^1\right).`$ (24)
Subtracting from the above expression the tadpole graph contribution as given by expression (14) and taking $`\mathrm{\Lambda }\mathrm{}`$, I obtain again Eq. (19) of Ref. . In the case in which $`𝒩`$ is odd, it is easy to show that also in this case it is recovered the main formula of Ref. .

It should be noted from Eqs. (10), (21) or (16) that the discretized sum over the eigenfrequencies in the presence of the soliton contains $`𝒩`$ terms less than the sum over the free eigenfrequencies, being $`𝒩`$, the number of discrete eigenfrequencies in the presence of the soliton. In the so called mode number regularization  this fact is used as the starting point. Here I have shown that it follows naturally from using adequately an energy momentum cut off regulator for the divergent sums in Eq. (1). In Ref. , in the case of the sine-Gordon model, in which there is only one discrete eigenvalue, it is not very clear why it is excluded the term with $`n=0`$ from the sum over the eigenfrequencies in the presence of the soliton. Here I have shown that such a term needs to be excluded since it does not satisfies the imposed periodic boundary conditions.

In Ref.  the authors also claim that the solitons in the field theoretical model, with density potential given by $`\varphi ^2\mathrm{cos}^2\mathrm{Ln}(\varphi ^2)`$, only exist semiclassically. The model $`\varphi ^2\mathrm{cos}^2\mathrm{Ln}(\varphi ^2)`$ was introduced in Ref. . The argument used by the authors of Ref.  is a one loop calculation performed in Ref. . In this reference the authors showed that the solitons of the triple degenerate $`\varphi ^6`$ model only exist at the classical level, at the quantum level an infrared divergence arises that can not be renormalized. This occurs because the potential of the Schrodinger equation for the eigenfrequencies in the presence of the soliton, takes different values at $`x\mathrm{}`$ and $`x\mathrm{}`$, that is, because the second derivative of the density potential at the neighbor minima are different. In the model, $`\varphi ^2\mathrm{cos}^2\mathrm{Ln}(\varphi ^2)`$ the second derivative in all minima are equal, but the third derivative in the neighbor minima are no longer equal. Since two loops quantum corrections are proportional to the third derivative of the density potential, the authors of Ref.  concluded trivially that at two loop order, something similar to what happens in the one loop calculation of Ref.  also will occur in the model $`\varphi ^2\mathrm{cos}^2\mathrm{Ln}(\varphi ^2)`$. In a future work I will show that this is not the case.

Acknowledgments
I would like to thank N. F. Svaiter for useful discussions and A. P. C. Malbouisson for reading the manuscript. Also a grant from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq-Brazil) is acknowledged.