text:
# Coupling of the lattice and superlattice deformations and hysteresis in thermal expansion for the quasi one-dimensional conductor TaS3

## Abstract

An original interferometer-based setup for measurements of length of needle-like samples is developed, and thermal expansion of o-TaS3 crystals is studied. Below the Peierls transition the temperature hysteresis of length $`L`$ is observed, the width of the hysteresis loop $`\delta L/L`$ being up to $`510^5`$. The behavior of the loop is anomalous: the length changes so that it is in front of its equilibrium value. The hysteresis loop couples with that of conductivity. The sign and the value of the length hysteresis are consistent with the strain dependence of the charge-density waves (CDW) wave vector. With lowering temperature down to 100 K the CDW elastic modulus grows achieving a value comparable with the lattice Young modulus. Our results could be helpful in consideration of different systems with intrinsic superstructures.

Internal degrees of freedom is a feature of a random system; in principle, they can give rise to metastable size states resulting, say, to hysteresis in thermal expansion. A special class form the compounds with intrinsic superstructures. Comprising two periodicities, generally incommensurate, the compounds occupy intermediate place between genuine aperiodic and truly periodic systems. In these systems, such as charge- and spin-density waves (CDW and SDW),  Wigner cristals, superconductors in magnetic fields, structurally incommensurate crystal phases,  the superstructure periodicities could be varied by external fields or temperature changes. The resulting metastable configurations can be reflected back onto the elastic properties and size of the underlying lattice, though this question is still poorly understood.

Quasi 1-dimensional conductors with CDW belong to a widely studied class of materials, in which intrinsic superstructure develops through the Peierls transition . When electrons condense into CDW they form a deformable medium, – an electronic crystal. Deformation of the CDW affects their main static and dynamic properties and gives rise to metastability and hysteresis.

The straightforward treatment of the CDW as a spring, whose strain is just applied to the crystal at the ends or via the impurities is not valid. Moreover, in the simple one-dimensional model the strains of the CDW and the crystal do not couple at all: if initially the CDW are relaxed, any change of the crystal length would not draw the CDW away from the equilibrium, i.e. give rise to a CDW deformation, as it was noticed in Refs. . Similarly, once the CDW is deformed, any change of the lattice constant, $`c`$, would neither decrease nor increase the deviation of the CDW wavelength $`\lambda `$ from the equilibrium value, $`\lambda _{eq}`$. So, within this model a CDW deformation would not give rise to a length change.

At the same time, the interaction of the CDW and the lattice is clearly seen from the elastic anomalies, including a drop of the Young modulus of the lattice,  $`Y_l`$, up to 4%, when the CDW become depinned. Mozurkewich  has concluded that the lattice deformation does give rise to a deformation of the CDW: when the CDW are at rest, they cannot relax, and so they contribute to the total elastic energy; the sliding CDW relax rapidly and their contribution drops out. For his model Mozurkewich  introduced an empirical parameter $`g`$ reflecting the deviation of $`\lambda _{eq}`$ from the simple one-dimensional model:

$$\delta \lambda _{eq}(c)/\lambda _{eq}(c_{eq})=(g+1)\delta c/c_{eq},$$
(1)
$`g0`$. However, no way was proposed to estimate $`g`$.

Another remarkable observation was reported in Ref. : the length $`L`$ of TaS3 samples as a function of electric field demonstrated hysteresis partly scaling with that of resistance. The values of $`L`$ obtained at different directions of the voltage sweep differed by $`10^6`$. This result was also treated in terms of coupling of the CDW strains with the deformation of the pristine lattice. Note that the field-induced length change (as well as the change of resistance) results only from an inhomogeneity of the sample properties: obviously, the length could depend on the voltage polarity only if the inversion symmetry of the sample is broken . Thus, study of the electric-field induced hysteresis of length can not provide complete understanding of the CDW-lattice coupling.

To study the effect of the CDW deformation on the crystal length it could be more fruitful to observe thermal hysteresis of $`L`$. With changing temperature the $`q`$-vector ($`q2\pi /\lambda `$) falls behind its equilibrium value, which is temperature dependent; the deviation of $`q`$ from the equilibrium is limited by the critical deformation at which the phase slippage (PS) begins. So, thermal cycling creates a CDW deformation, which is relatively uniform along the sample and is of the maximum possible value.

For our study we chose o-TaS3 as a representative CDW conductor. Among the quasi one-dimensional conductors TaS3 is one of the most widely studied, including the elastic properties. TaS3 demonstrates the Peierls transition at $`T_P=220`$ K, below which the resistance follows an activation law revealing the half-gap $`\mathrm{\Delta }=700`$ K. The dependence $`R(T)`$ demonstrates a pronounced hysteresis, the width of the loop $`\delta R/R`$ being up to 50 % at $`T`$ around 100 K.

In the present Letter we report temperature hysteresis of length for o-TaS3; the hysteresis couples with that of resistance having an anomalous sign. A quantitative treatment of the effect is proposed.

We studied samples of TaS3 with typical length 1 mm, width 15–35 $`\mu `$m, and thickness 5-10 times less. The samples were arranged on a transparent glass substrate (Fig. 1a). The contacts were fixed on the substrate with indium, while the central part of the sample formed an arc. The inner surface of the substrate played the role of a semi-transparent mirror. The laser beam with wavelength $`\lambda _l=635`$ nm fell down through the substrate, and partially reflected from the inner surface of the substrate and from the sample surface forming an interference pattern (similar with the wedge interference). The pattern was fixed with a video camera combined with a microscope. An example of such a pattern is shown in Fig. 1b. The neighboring dark or bright fringes correspond with change of the sample distance $`y`$ from the substrate surface by $`\pm \lambda _l/2`$. Making sections of such images (an example is shown in Fig. 1c) we obtained the profiles $`y(x)`$ of the samples, like the curve presented in Fig. 1d. Finally, the length of such curves was calculated. We detect a length change about 5 Å which corresponds to the relative length change $`\delta L/L=510^7`$, though in principle the sensitivity could be improved. Our estimates have shown that the contributions to $`\delta L`$ due to the changes of stress and $`Y_l`$ are negligible. Thus, the observed changes of $`L`$ with temperature were associated with the change of the equilibrium sample length with respect to the substrate. To have absolute results we have measured thermal expansion of the glass substrate in addition.

Figs. 2a) and b) show the temperature dependences of length (with respect to glass) and of the resistance for the representative sample of TaS3. In the inset to Fig. 2a) the absolute dependence of length is shown. The dependences of length and resistance both clearly show hysteresis below the Peierls transition temperature $`T_P=214`$ K. The length hysteresis loop opens immediately below $`T_P`$. The development of hysteresis is the dominating effect near $`T_P`$, on whose background it is difficult to distinguish a feature coupled to the transition in itself. The length hysteresis grows below $`T_P`$ achieving $`\delta L/L=510^5`$ (Fig. 2a), which by 1.5 orders of magnitude exceeds the maximum value reported in Ref. . At the same temperature the states approached from lower temperature have higher length than those approached from above. Correspondingly, application of electric field exceeding the threshold reduces the length after heating the sample, and increases after cooling. Thus, the main effect of electric field is the relaxation of the thermally induced states. The length change induced by the field in itself was $`10^6`$, like in Ref. . Note, that the sign of the hysteresis is non typical: the length is in front of its equilibrium value (inset to Fig. 2a), but not behind.

The length hysteresis apparently couples with the resistance hysteresis, higher length corresponding to higher resistance at a given temperature. The coupling between $`R`$ and $`L`$ is well seen from Fig. 2c, where the temperarure dependence of $`\delta L/L`$ is presented together with the loop of conductivity, $`\delta \sigma 1/R_{cool}1/R_{heat}`$. Note that in the unipolar approximation $`\delta \sigma \sigma _{300}(\delta \lambda /\lambda )(\mu /\mu _{300})`$,  where the index “300” marks the room–temperature values. So, the $`\sigma `$ hysteresis reflects that of the CDW wavelength to the accuracy of the temperature dependence of the mobility $`\mu `$, which in turn can be found from Ref. . Fig. 2d shows the ratio $`(\delta L/L)/(\delta \lambda /\lambda )`$ vs. $`T`$. At the low temperatures $`\delta L/L`$ is $`6`$ % of the CDW deformation.

As we have noticed in the beginning, the interactions of the CDW and the sample cannot be presented as that of two springs connected in parallel (the lower inset to Fig. 3). Though we do not see a way to obtain quantitatively the strain dependence of the $`q`$-vector for TaS3, one can recall experimental results, from which it is possible to estimate the value of $`g`$ (Eq. 1). It has been found  that at a certain value of uniaxial strain $`S^{}`$ the properties of the CDW change abruptly: e.g., the non-linear conduction nearly disappears. This has been attributed to the transition of the CDW to the 4-fold commensurability. Comparing the normalized dependences of the $`q`$-vector  and of the strain $`S^{}`$  on temperature (Fig. 3) one can see that they become similar  if one multiplies the value of $`S^{}`$ by 6. Thus, one can assume that a deformation of the lattice induces the change of the equilibrium CDW wavelength $`\lambda _{eq}(c)`$ in accordance with Eq. (1) with $`g6`$. Such a large value of $`g`$ indicates that the simple 1-dimensional model  cannot even roughly describe the strain-induced change of $`\lambda _{eq}`$. Evidently, the change of $`\lambda _{eq}`$ is dominated by transverse effects: the longitudinal strain decreases the thickness of the sample resulting in an increase of the intrachain coupling. This modifies the form of the Fermi surfaces, and consequently, the $`q`$ -vector. To comprehend the sign and the value of the $`q`$ change one should study the evolution of the Fermi surface in detail. In principle, we cannot exclude that $`g`$ is temperature dependent. However, as we shall see below, the assumption that $`g=`$ const, namely $`g6`$, is consistent with our experiment.

For the next step, we present the elastic energy density $`W`$ as a sum of the lattice and the CDW energies:

$$W=\frac{1}{2}(Y_l(\frac{cc_{eq}}{c_{eq}})^2+Y_c(\frac{\lambda \lambda _{eq}(c)}{\lambda _{eq}})^2),$$
(2)
where $`Y_l`$ and $`Y_c`$ are the elastic moduli of the lattice and the CDW respectively. Taking into account the condition (1) and assuming that no PS occurs ($`\lambda /c=\mathrm{const}`$), we can minimize $`W`$ and obtain the resulting length change:

$$\frac{\delta L}{L}=\frac{\delta c}{c_{eq}}=g\frac{Y_c}{Y_l+g^2Y_c}\frac{\delta \lambda }{\lambda }g\frac{Y_c}{Y_l}\frac{\delta \lambda }{\lambda },$$
(3)
where $`\delta \lambda `$ is the initial CDW deformation (at fixed $`c`$). The approximation implies that $`g^2Y_cY_l`$. Eq. 3 is quite transparent: the crystal deforms as a spring connected in parallel to the CDW, but with a factor $`g`$. $`g>0`$ means that, say, compressed CDW ($`\lambda <\lambda _{eq}`$) would result in a decrease of the sample length, which agrees with our observation. E.g., cooling corresponds to growth of $`\lambda _{eq}`$ , so the CDW is in compressed state ($`\lambda \lambda _{eq}<0`$), which results in decreased $`L`$ (Fig 2a). The CDW – crystal interaction could be illustrated with a scheme consisting of two springs connected via a lever (Fig. 3, upper inset).

From Eq. 3 it follows that the ratio $`(\delta L/L)/(\delta \lambda /\lambda )`$ (Fig. 2d) equals $`gY_c/Y_l`$. Note, however, that with $`g=6`$, $`g^2Y_c`$ at low temperatures is comparable with $`Y_l`$, and one should use the exact version of Eq. 3. The resulting value of $`gY_c/Y_l`$ is also shown in Fig. 2d. One can see the growth of the CDW elastic modulus with lowering temperature, though slower than $`Y_c\mathrm{exp}(\mathrm{\Delta }/T)`$  (see the solid line in Fig. 2d). Note that the curve $`(\delta L/L)/(\delta \lambda /\lambda )`$ vs. $`T`$ resembles the temperature dependences of the lattice softening in electric field . This could be expected, since from Eqs. 1 and 2 it follows that the value $`g^2Y_c`$ is added to $`Y_l`$ when the CDW are at rest, while in the sliding state the CDW contribution drops out .

From the condition of neutrality  accompanying the CDW deformation we can obtain for the CDW Young modulus: $`Y_c=\frac{d\zeta }{dq}\frac{q^2}{\pi s}`$ , where $`\frac{d\zeta }{dq}`$ is the derivative of the chemical potential level by the CDW $`q`$-vector, and $`s`$ is the area per chain. Substituting this in Eq. 3, we obtain:

$$\frac{\delta L}{L}g\frac{2\delta \zeta }{\lambda sY_l}.$$
(4)
For the unipolar (p-type for TaS3) conduction  $`\delta \zeta =T\mathrm{ln}(R_{heat}/R_{cool})`$. Substituting into Eq. 4 $`g=6`$, $`\delta \zeta =50`$ K (for $`T=120`$ K, Fig. 2b), $`\lambda =10`$ Å, $`s=80`$ Å2 (Ref. ), $`Y_l=380`$ GPa , we obtain $`\delta L/L=310^5`$, in agreement with the experiment (Fig. 2b). This agreement supports the approach we used to estimate $`g`$, and proves that under a certain uniaxial strain ($`S^{}`$) the CDW in fact transit to the 4-fold commensurability.

In conclusion, we have observed thermal hysteresis of length for the quasi one-dimensional conductor TaS3. In the framework of the model proposed the observed anomalous sign and value of the hysteresis loop are consistent with the dependence of the $`q`$-vector on the lattice strain. The results could be helpful in understanding the behavior of other systems with intrinsic superstructures, whose parameters could be rearranged due to the temperature effects or external fields.

We are thankful to G.P. Vorob’ev, S.G. Zybtsev, and S.V. Zaitsev-Zotov for help in the experiment and discussions, to S.N. Artemenko, A.A. Sinchenko, V.B. Preobrazhenskii, G.V. Stepanov, J.W. Brill, and M.E. Itkis for helpful discussions. We acknowledge RFBR (02-02-17301, 00-02-22000 CNRS), the State programs ”Physics of Solid-State nanostructures” and ”Low-Dimensional quantum structures”, NWO and INTAS-01-0474.


format: text

text:
# Granularity-induced gapless superconductivity in NbN films: evidence of thermal phase fluctuations

## Abstract

Using a single coil mutual inductance technique, we measure the low temperature dependence of the magnetic penetration depth in superconducting NbN films prepared with similar critical temperatures $``$16 K but with different microstructures. Only (100) epitaxial and weakly granular (100) textured films display the characteristic exponential dependence of conventional BCS s-wave superconductors. More granular (111) textured films exhibit a linear dependence, indicating a gapless state in spite of the s-wave gap. This result is quantitatively explained by a model of thermal phase fluctuations favored by the granular structure.

Distinctive property of simple superconducting materials like metals and alloys is the absence of thermally activated electronic excitations at low temperature $`TT_c`$, where $`T_c`$ is the critical temperature. This property is consistent with a model of superconducting gap $`\mathrm{\Delta }_\text{k}`$ with s-wave symmetry, as predicted by the BCS theory. In contrast, in complex materials like heavy fermions gro , cuprates har  or rhutenates bon , thermal excitations of arbitrarily small energies are experimentally observed. This gapless property has been commonly explained by an excitation spectrum with nodal points or lines, characteristic of $`\mathrm{\Delta }_\text{k}`$ with p\- or d-symmetry respectively. According to theoretical considerations, non-s symmetries are favored by non-phonon mediated pairing mechanisms like antiferromagnetic spin fluctuations mon . However, experimental and theoretical studies show that a gapless state can be induced in s-wave superconductors by residual interactions coexisting with the pairing interaction, such as magnetic impurities abr , spin density waves ove , proximity effects cla ; pam  or phase fluctuations of the order parameter cof ; rod ; eme . In conclusion, a non-s symmetry implies a gapless state, but the opposite does not hold.

In this work we investigate whether the above predictions of gapless state apply to the s-wave superconductor NbN. By measuring the variation of the magnetic penetration depth $`\lambda (T)\lambda (0)\delta \lambda (T)`$ at $`T`$ 1.5 K, we study the low energy excitations on films with different microstructures obtained by varying growth conditions. NbN is a model system for our purpose taking into account the following: i) the simple cubic structure with fcc lattice bec ; ii) the conventional BCS s-wave superconducting properties pam2 ; iii) the relatively high critical temperature $`T_c16`$ K enables to study the low energy excitations in a wider temperature region; iv) various materials preparation routes have been developed and optimized, especially in view of electronics applications bac ; vil1 ; vil2 .

For this work we have selected six NbN films prepared by DC magnetron sputtering with similar $`T_c`$ close to the bulk value $``$16 K but with different degrees of granularity obtained by varying growth conditions. These conditions are reported in Table I (see vil1 ; vil2  for more details): growth temperatures are either room temperature or 620 C; substrates are (1$`\overline{1}`$02) Al2O3 or (100) MgO single crystals or (1$`\overline{1}`$02) Al2O3 single crystals with a 12 nm thick (100) MgO buffer layer; NbN thicknesses range from 10 to 1400 nm. Three different types of microstructure are reproducibly obtained: i) (100) epitaxial; ii) highly (100) textured; iii) (111) textured. Using (100) MgO buffer layers, a high degree of (100) texture is obtained only in sufficiently thin ($`\text{d}`$ 150 nm) NbN films vil2 .

The film microstructure is characterized by X-ray diffraction in the Bragg-Brentano geometry, X-ray pole figure measurements and transmission electron and atomic force microscopies (TEM, AFM). X-ray rocking curve analysis indicates that samples A1043 and A1013 are epitaxial with FWHM of the (420) Bragg refelection less than 1. The former sample is thin (10 nm), whereas the latter is thick (950 nm). The film on the MgO buffer layer (A1044) is thin (10 nm) and highly (100) textured. The remaining three films on sapphire (A1057, A1060 and A1063) are all thick (280-1400 nm) and display predominant (111) texture. The different degrees of granularity in these films are confirmed by AFM and TEM. For example, in sample A1043, we find an atomically flat surface lam  and a crystal structure with apparently no grain boundaries (see fig. 1a), which indicates that this film is monocrystalline. A granular structure is evident in the (111) textured sample A1060 from the TEM cross-sectional micrograph of fig. 1b. We note a columnar structure with grains of size 15-35 nm and well-defined grain boundaries, as previously reported for films grown under the same conditions vil1 .

The microstructural properties of the six films are correlated with DC resistivity $`\varrho `$ and DC transport critical current density $`j_c`$ measurements. The results of these measurements are summarized in Table I. In all films, $`T_c`$ is close to the bulk value $``$16 K, which indicates that the microstructure does not appreciably perturb the superconducting phase. Smaller residual resistivities $`\varrho `$(20 K) $``$ 50$`\mu \mathrm{\Omega }`$cm and larger resistivity ratios $`\varrho (300\text{ K})/\varrho (20\text{ K})`$ 1.2 are found in the epitaxial films. The highest values of $`j_c`$ (4.2 K) $``$ 3 MA/cm2 are found in the epitaxial or thin (100) textured samples, while the three thick (111) textured ones exhibit values $``$ 1 MA/cm2 that decrease with increasing NbN thickness. In agreement with earlier studies vil1  and with our TEM results, we conclude the following: i) the epitaxial samples have higher crystalline quality; ii) in the (111) textured samples, the degree of granularity - defined as the degree of grain decoupling - increases with increasing thickness.

To correlate the above structural and transport properties with the low energy excitations, we have carried out $`\delta \lambda (T)`$ measurements in the range 1.5 K-$`T_c`$ using a single-coil mutual inductance technique described elsewhere gau . Typical frequency and magnitude of the e.m. field perpendicular to the film surface are 2-4 MHz and $`<`$ 0.1 mT respectively. Care has been taken to always operate in the linear response regime, i.e. below the Josephson critical field $`B_{c1}`$, monitored by varying the applied magnetic field. Note that we measure an effective magnetic penetration depth that, in the dirty limit, differs from the London penetration depth by a correction factor of the order of unity tink . This limit is appropriate for our films, for which we estimate a mean free path $``$1 nm and a coherence length $`\xi 35`$ nm vil1 ; lam . As a good approximation, the above correction factor is temperature independent in the limited temperature range of our measurements. Therefore the temperature dependence of the variations $`\delta \lambda (T)`$ is not affected by dirty limit effects.

In figs. 2a-b, we report the experimental $`\delta \lambda (T)`$ curves obtained. Only the epitaxial and thin (100) textured films display the exponential behavior $`\delta \lambda \text{exp}(\mathrm{\Delta }/k_BT)`$, characteristic of BCS s-wave superconductors bar  (see fig. 2a). All of the three thick (111) textured samples exhibit a linear behavior which can not be approximated by any power law (see fig. 2b); the range of linearity is limited below $``$ 2.7 K in the thinner film A1057 and extends up to $``$4 K in the thicker films A1060 and A1063. A BCS s-wave fit perfectly accounts for the experimental data of the three samples of fig. 2a. $`\lambda (0)`$ is the only free parameter of this fit, whereas the ratio 2$`\mathrm{\Delta }(0)/k_BT_c`$ is kept fixed at 4.6, taken from tunneling data at 1.5 K on sample A1043. The fitted $`\lambda `$(0) values are in the range 200-400 nm, in agreement with previous reports pam2 ; vil2 . A BCS fit is satisfactory also for the thinner (111) textured film A1057 of fig. 2b; indeed, the behavior of this sample is intermediate between exponential and linear. For the remaining two samples A1060 and A1063, no BCS fit is possible for any plausible set of parameters.

These results give evidence of a gapless superconducting state in the (111) textured films. We recall the mechanisms that may induce a gapless state in a conventional BCS s-wave superconductor: i) magnetic impurities; ii) proximity effects; iii) phase fluctuations of the order parameter. In our case, the first two possibilities are excluded: i) the purity of the Nb target used is 99.995 % and traces of magnetic impurities are below 10 ppm. We conclude that the amount of these impurities in the films is negligible. In any case, the gapless samples contain no more impurities than the others, since all of them are prepared with the same Nb target; ii) a previous TEM and tunneling study vil1  indicates that the intergranular phases are insulating, thus no proximity effects are expected. The third picture remains plausible: as suggested earlier deu ; sim1 ; cof ; rod ; eme ; cha , both quantum and thermal phase fluctuations can be relevant in granular superconductors with short coherence length $`\xi `$ like our NbN films.

To verify the validity of the above picture in our case, we consider an array of superconducting grains with intergranular Josephson coupling energy $`E_J`$ and charging energy U, as first proposed by Abeles abe , where U is associated with the excess of Cooper pairs n in each grain. U and $`E_J`$ are related to the intergranular capacitance C and the critical current $`I_c`$ respectively, according to the relations $`U=2e^2/C`$ and $`E_J=\mathrm{}I_c/2e`$. Following previous work sim1 ; mcl , we assume the array to form a square lattice of constant a equal to the average grain size. The resulting Hamiltonian is:

$`H={\displaystyle \frac{U}{2}}{\displaystyle \underset{i}{}}n_{i}^{}{}_{}{}^{2}+E_J{\displaystyle \underset{ij}{}}[1\text{cos}(\vartheta _i\vartheta _j)]`$ (1)
where $`\vartheta _i`$ is the phase of the order parameter in the i-th grain and $`n_i=2i/\vartheta _i`$ is its conjugate operator.

Eq. (1) describes both quantum and classical phase fluctuations. Two arguments suggest a classical description of our results in the range $`T`$ 1.5 K of our experiment. i) According to theoretical studies cof ; cha  supported by numerical calculations rod , quantum fluctuations lead to a reduction of the renormalized superfluid density, implying a progressive decrease of $`\lambda /T`$ as $`T`$ 0 K. This effect is not observed in the curves of fig. 2b which remain linear above 1.5 K. ii) In our (111) textured films, we estimate $`UE_J`$, which corresponds to a purely classical regime sim1 ; cof ; cha . Indeed, typical figures are: $`a`$ 20 nm, $`\lambda (0)`$ 300 nm, $`j_c`$ 0.5-1 MA/cm2, thickness and relative dielectric constant of the intergranular regions respectively $`t`$ 2 nm and $`ϵ_r`$ 2. Thus, we obtain $`E_J`$ 50-100 meV and $`U`$ 10 $`\mu `$eV.

In the classical regime, the fluctuation strength is determined by the ratio $`E_J/k_BT_c`$, where $`T_c`$ is the mean field critical temperature. Using the above estimate, we find $`E_Jk_BT_c`$ in our granular films. Thus, according to previous studies sim1 ; mcl , we expect small fluctuations which do not inhibit long range phase order at any temperature. Indeed, we do not observe any appreciable decrease of $`T_c`$ induced by granularity ($`T_c`$ is even slightly larger in the (111) textured films). These fluctuations coexist with the usual quasiparticle excitations predicted for BCS s-wave superconductors. Both excitations lead to a progressive increase of $`\lambda `$ as temperature increases. In the Gaussian approximation, valid at low temperatures $`TE_J/k_B`$ and in the absence of dissipation, only long wavelength longitudinal fluctuations of arbitrarily small energy are created (Goldstone modes). This produces a linear increase $`\delta \lambda T`$ cof ; rod ; eme  experimentally observed in fig. 2b. On the other hand, the creation of quasiparticles requires a finite thermal energy $`\mathrm{\Delta }`$, which leads to the exponential dependence of fig. 2a.

We now analyze in quantitative form the experimental curves of fig. 2b within the framework of the above picture. Using Table I, we estimate first the fluctuation contribution and then the temperature at which the quasiparticle contribution becomes significant. Following previous calculations rod ; cof , in the Gaussian regime, our lattice model with U=0 leads to:

$`{\displaystyle \frac{\lambda }{T}}={\displaystyle \frac{ek_B}{\mathrm{}}}{\displaystyle \frac{1}{j_ca}}`$ (2)
where $`I_cj_ca\lambda `$ at 0 K and we neglect the weak temperature dependence of this product at low temperature. In Table I we report the values of $`\lambda /T`$ estimated using eq. (2) and compare them with the experimental values. Taking into account the uncertainty in the estimate of a, we notice a good agreement. Interestingly, the trend of $`j_c`$ in the (111) textured films follows the trend of $`\lambda /T`$, as predicted by eq. (2): the thicker the film, the smaller $`j_c`$ and the larger the slope $`\lambda /T`$. In the films of fig. 2a, exhibiting a BCS s-wave behavior, $`j_c`$ is larger. Thus, the fluctuation effect vanishes and only the quasiparticle effect is observed. To estimate the temperature at which the two effects become comparable, we separately evaluate their contribution to $`\lambda /T`$ using Table I and the BCS fit of fig. 2b. We find that the linear increase of $`\lambda `$ produced by thermal phase fluctuations is expected to dominate at $`T`$ 4 K. This is in agreement with the point of departure from linearity in samples A1060 and A1063. As discussed before, in sample A1057 a gapless behavior is visible but does not dominate yet above 1.5 K, since a BCS fit is still satisfactory.

In conclusion, we have reported for the first time experimental evidence of gapless state in a conventional BCS s-wave superconductor, NbN. Evidence is found in granular films exhibiting a clear linear dependence of $`\lambda `$ in the range 1.5-4 K. Epitaxial or less granular films exhibit the conventional dependence of s-wave superconductors and no trace of gapless state is found. A simple model of granularity-induced thermal phase fluctuations of the order parameter (Goldstone modes) quantitatively accounts for the above linear dependence. No signature of quantum fluctuations is observed above 1.5 K, as expected taking into account the relatively large grains of our films.

Our evidence of granularity-induced gapless superconductivity in NbN raises the question whether this phenomenon is observable in other granular superconductors. In cuprates, which display gapless properties in agreement with a d-wave symmetry model, it should be investigated if these properties could be - at least partially - explained by intrinsic granular properties associated with chemical or electronic inhomogeneities in the coherence length scale $``$ 3 nm. In any case, our results suggest that the only reliable methods for testing non-s symmetries are those based on phase-sensitive measurements.

###### Acknowledgements.

We acknowledge A. Andreone, A. Barone, B. K. Chakraverty and G. Deutscher for stimulating discussions and A. Dubon for assistance in the TEM work.

Corresponding author. E-mail address: gauzzi@maspec.bo.cnr.it

Present address is: INFM and Dipartimento di Scienze Fisiche, Università di Napoli Federico II, 80125 Napoli, Italy.

§Present address is: Laboratoire Pierre Süe, DRECAM, CEA-Saclay, 91191 Gif sur Yvette, France.


text:
# Progress in lattice algorithms

## 1 INTRODUCTION

This conference illustrated that most large collaborations actively investigating QCD physics are performing simulations of the theory with dynamical quark fields. Dynamical staggered fermion simulations have been pursued for many years now, however it is clear that the community is now also firmly in the era of full Wilson fermion studies. These calculations will address one of the biggest systematic error in many non-perturbative QCD predictions from the lattice; use of the “quenched approximation”. At present, most simulations are not being performed at realistic values of the quark masses, or with the correct number of flavours of light fermions. A technically challenging extrapolation of data to the up and down quark masses must subsequently be performed.

At the conference, statements were given  by the biggest collaborations. Estimated costs for reproducing the state-of-the-art in quenched spectrum data in the full theory were in the range 10-100 Tflop-years. In this review, efforts to develop better algorithms to simulate closer to the physical parameters of QCD are discussed.

## 2 ALGORITHMS FOR WILSON FERMION SIMULATIONS

In this section, the main features of the techniques for simulating dynamical Wilson fermions are summarised. The simulations of dynamical staggered fermions are more mature; the principle difficulty at present is one of controlling the choice of discretisation . The new challenge is to find efficient methods for including recently proposed improved discretisations. Developing algorithms incorporating “fat-link” gauge fields, designed to enhance the flavour symmetry of the action are presented later. It is worth noting that the exact algorithms developed for odd-flavour simulations with Wilson fermions are being applied to staggered fermion simulations.

The Hybrid Monte Carlo (HMC) algorithm  is still the work-horse powering most large Wilson fermion simulations. Over the past few years there have been further advances in our understanding of the issue of finite precision arithmetic and the reversibility and stability of the molecular dynamics (MD) integrator component. HMC is an exact algorithm if the MD integrator is reversible and area preserving. It has long been understood that the Hamiltonian evolution of QCD is chaotic, and the Liapunov exponents were studied in detail . More recently, the stability of the MD leap-frog integrator  has been investigated in detail by the UKQCD collaboration . Their findings were presented to the conference by Joó . Analysis of a leap-frog integrator of the simple harmonic oscillator reveals the solution diverges once the step-size exceeds a critical value. The problem is exacerbated for higher-order integrators. This same behaviour is observed in QCD, and the critical step-size is seen to fall as the quark mass, $`m_q0`$. The suggestion is that the problem of simulating light Wilson fermions with HMC will amount to managing this stability problem.

There has been progress in implementations of the local boson algorithm (LBA) . This algorithm relies on first constructing a polynomial approximation to the inverse of the fermion matrix, then expressing the polynomial as a product of roots. The determinant is expressed as a local bosonic partition function involving a large number (the degree of the polynomial) of auxiliary fields. The beauty of the method is that a local change in the gauge fields leads to a local computation for the change in the action. Unfortunately, the simplest implementations of this method suffer from long autocorrelations ; a large number of auxiliary fields are coupled to the gauge fields, and so only small changes in these degrees of freedom can be made in each sweep. The problem can be circumvented by using lower-order polynomials, and correcting for any discrepancy between the approximation and the true inverse using a noisy accept-reject step . This idea lead to the development of the two-step multi-boson (TSMB) algorithm , where two polynomials are constructed, $`𝒫_1(x)`$ (low-order) and $`𝒫_2(x)`$ (high-order) such that

$$𝒫_1(x)𝒫_2(x)\frac{1}{x^{N_f/2}}.$$
(1)
$`𝒫_1(x)`$ is replaced by a partition function of a small number of boson fields, and $`𝒫_2(x)`$ is used in the acceptance test.

Polynomial approximations have been incorporated into other techniques. The Polynomial Hybrid Monte Carlo (PHMC) algorithm  proposed by Frezzotti and Jansen is one such development. A single auxiliary field, $`\varphi `$ is introduced, rather than the large number of LBA, and the action is

$$S_\varphi =S_G[U]+\varphi ^{}𝒫(M^{}M)\varphi $$
(2)
This avoids the stiff update dynamics, however the action is now non-local, so HMC is used to update the gauge fields. The algorithm can be made exact either via an accept-reject step or by a re-weighting of expectation values by a correction factor:

$$𝒪_{\mathrm{true}}=\frac{𝒪detQ^2𝒫(Q^2)_\varphi }{detQ^2𝒫(Q^2)_\varphi }.$$
(3)
Frezzotti and Jansen advocate this re-weighting and the break-down of the approximation close to zero is now regarded as advantageous. The polynomial weight generates more configurations with low eigenvalues (which are then assigned a lower weight in the ensemble average), and this over-sampling should lead to a more reliable determination of quantities that are sensitive to the details of the lowest eigenvalues.

The ALPHA collaboration have made a comparison  between PHMC and pseudofermionic HMC. In the range of parameters of interest, they find there is a slight advantage to PHMC and favour this method for future simulations, on account of its over-sampling property. Schroers presented results to the conference  that suggest the performance of the tuned LBA is faster than the SESAM production HMC code. Performance comparisons between TSMB and HMC in compact QED were presented to the conference by Zverev . The study finds the two algorithms are equally competitive, although the authors emphasise their TSMB implementation can still be optimised by changing the local gauge update scheme and tuning the polynomials $`𝒫_1`$ and $`𝒫_2`$ more carefully.

One remarkable feature of these comparisons is that the results are, in fact, so close. de Forcrand emphasises that the LBA algorithm turns into a standard local heatbath/over-relaxation algorithm as the fermions are made heavier, and the expectation there is that the over-relaxation algorithms should out-perform (global) HMC significantly.

## 3 ODD-FLAVOUR SIMULATIONS

Since QCD has three flavours of light fermions, there is a natural incentive to study algorithms for simulating an odd number of fermion flavours. The standard pseudofermion formulation is unsuitable, since the Wilson matrix can have eigenvalues with negative real parts, meaning there are regions of configuration space where the gaussian integral is not defined. One extremely useful development arising from the local boson method has been the possibility of performing these odd-flavour simulations, even maintaining an exact update algorithm. Since these simulations still rely on gaussian integrals, they in fact generate configurations with probability measure defined from the modified partition function,

$$Z_+=𝒟U|detM|e^{S_G}.$$
(4)
The sign of the determinant can be included in a post-hoc reweighting of observables:

$$𝒪_{\mathrm{true}}=\frac{𝒪\text{ sgn}(detM)_+}{\text{sgn}(detM)_+}.$$
(5)
This may lead to a “sign problem”, where statistical estimates of reweighted observables fluctuate wildly if a mixture of configurations with both positive and negative sign determinants occur in the ensemble. Empirical studies of this problem for new simulations of QCD with three quarks (with a target of $`m_qm_s/4`$) were presented to the conference by Gebert . They expect the sign problem to be very mild in this region.

The local boson algorithm can approximate the partition function of Eqn. 4; starting from a polynomial approximation to $`1/\sqrt{x}`$ in the range $`[ϵ,1]`$ leads directly to the one-flavour partition function. Polynomials of the non-hermitian fermion matrix $`M`$ also present a solution, as suggested in Ref.  and developed in Ref. . A polynomial approximation to $`M^1`$ of order $`2n`$ is constructed,

$$M^1𝒫(M)=\underset{k=1}{\overset{2n}{}}(Mz_k).$$
(6)
For a suitably chosen polynomial, the roots, $`z_k`$ come in complex-conjugate pairs, so the polynomial can be written

$$𝒫(M)=\underset{k=1}{\overset{n}{}}(Mz_k)(Mz_k^{}).$$
(7)
The $`\gamma _5`$ hermiticity of $`M`$ means

$$det(Mz_k^{})=det(M^{}z_k^{})$$
(8)
so $`|detM|`$ can be represented approximately by a gaussian integral,

$$|detM|𝒟\varphi 𝒟\varphi ^{}\mathrm{exp}\left\{\varphi ^{}T_n^{}(M)T_n(M)\varphi \right\}$$
(9)
with

$$T_n(M)=\underset{k=1}{\overset{n}{}}(Mz_k).$$
(10)
The scheme can be made exact with a Metropolis test; two alternatives were proposed in Refs , both of which rely on a noisy Kennedy-Kuti  acceptance test. The three-flavour simulations of QCD are being carried out by the JLQCD collaboration and results were presented to the conference .

## 4 ACCELERATING HMC

HMC remains the most popular technique for generating ensembles of gauge fields with two flavours of dynamical Wilson (and Sheikholeslami-Wohlert improved) fermions. Most of the simulations use this well-established method, so it is useful to find simple schemes to enhance the method.

### 4.1 ILU preconditioning

Replacing the fermion determinant by a gaussian integral couples extra stochastic degrees of freedom to the gauge fields. This in turn leads to a lower Metropolis acceptance rate in HMC. It has been known for a long time that coupling the gauge fields and pseudofermions with the even-odd preconditioned fermion matrix (whose determinant is equal to that of the original fermion matrix) leads to a higher acceptance rate for a given molecular dynamics step-size . The SESAM collaboration  emphasised that the even-odd preconditioning scheme is an example of an incomplete LU factorisation technique for matrix preconditioning. ILU preconditions a matrix by left and right multiplication by two readily invertible matrices;

$$\stackrel{~}{M}=(IL)^1M(IU)^1$$
(11)
where L and U are the lower and upper sections of the matrix. Defining L and U first requires sites on the lattice, $`𝐱`$ are assigned an index, $`s(𝐱)`$ then site ordering is defined as

$$𝐱>𝐲\text{ if }s(𝐱)>s(𝐲).$$
(12)
SESAM developed a site ordering that is well-suited to accelerate matrix inversion on parallel computers. Since this preconditioning leaves the determinant of the matrix invariant, $`\stackrel{~}{M}`$ can be used to couple the pseudofermions to the gauge fields. This holds for any choice of ordering function, $`s(𝐱)`$. In Ref. , different site orderings were tested in HMC simulations of the two-flavour Schwinger model. An unusual pattern emerged; the optimal ordering for speeding up matrix inversion was the global lexicographic scheme, while the best ordering for improving Metropolis acceptance was the usual even-odd preconditioner. In Ref. , it was noted that the even-odd matrix itself can be ILU factorised again, and this both accelerates matrix inversion and improves the Metropolis acceptance rate. Ref.  tested a range of two step (eo-ILU) preconditioned matrices (defined now by different choices of an indexing function, $`s(𝐱_e)`$ on the even sub-lattice only). These simulations found the same pattern emerged; the optimal ordering for inversion is a globally lexicographic one, while using a “locally ordered” scheme is best for improving Metropolis acceptance. Using an eo-ILU ordered interaction allows a three times larger step-size for a given Metropolis acceptance rate. The naive expectation that the autocorrelations of the update scheme are independent of the integrator step-size at a fixed acceptance is supported by simulation data.

### 4.2 Splitting the pseudofermions

Hasenbusch  recently proposed modifying the pseudofermion sector of HMC to improve the acceptance rate for a fixed molecular dynamics step-size. The key observation is (as above) that coupling the gauge fields to pseudofermions via a well-conditioned matrix allows a larger molecular dynamics step-size to be taken for a fixed Metropolis acceptance. The fermion determinant is first split into two pieces, namely

$$detM=detM\overline{M}^1det\overline{M}.$$
(13)
The two-flavour determinant is then represented by a gaussian integral,

$$detM^2=𝒟[\varphi ,\varphi ^{}\psi ,\psi ^{}]\mathrm{exp}\left\{S_\varphi S_\psi \right\}$$
(14)
with

$$S_\varphi =|\overline{M}M^1\varphi |^2\text{ and }S_\psi =|\overline{M}^1\psi |^2.$$
(15)
Hasenbusch then chooses $`\overline{M}=I\stackrel{~}{\kappa }\mathrm{\Delta }`$ and recognises that for $`0<\stackrel{~}{\kappa }<\kappa `$, the two matrices $`M\overline{M}^1`$ and $`\overline{M}`$ are better conditioned than $`M`$. This leads to improved acceptance.

In Ref , Hasenbusch tested the idea in the two-flavour Schwinger model and at the conference , results for preliminary studies in QCD were presented. The cost of HMC simulations (defined in units of matrix $`\times `$ vector operations) as a function of the splitting parameter, $`\stackrel{~}{\kappa }`$ is summarised in Fig. 2 In the Schwinger model study, simulations are accelerated by a factor of 1.7 when the pseudoscalar meson mass, $`m_P=0.210(3)`$, rising to a speed-up of 2 when $`m_P=0.124(5)`$. This increase at lighter fermion mass is an encouraging result. The studies in QCD  are more exploratory, but effects of similar magnitude are seen.

### 4.3 Multiple time-scale integration

Based on the idea of separating the infra-red and ultra-violet sectors of the fermion matrix, Jim Sexton and I are working on a practical implementation of a multiple time-scale molecular dynamics integrator. A low-order polynomial approximation to the fermion matrix inverse, which captures the short distance scale effects of quark loops is used, along with a pseudofermion action to reproduce long-range physics. The force on the gauge field from the polynomial term in the action is integrated (along with the easy-to-compute force from the Yang-Mills discretisation) using a finer time-step, and the long-range modes, contained in the pseudofermion action are integrated using a larger time step .

### 4.4 Are these schemes compatible?

If these ideas could be combined, the possibility of an order-of-magnitude acceleration of Wilson fermion simulations would open up. ILU preconditioned pseudofermions (Sec. 4.1) can be used straightforwardly in the multiple time-scale scheme of Sec. 4.3. This idea is under investigation, and it seems likely the benefits from each component will be combined. Hasenbusch demonstrated directly that his scheme works with even-odd preconditioned pseudofermions.

## 5 NEW DEVELOPMENTS

The Kentucky group  are developing a new dynamical fermion updater based on a Kennedy-Kuti noisy acceptance test on changes in the determinant. The ratio of the determinants evaluated on the new and old configurations is computed using a $`Z_2`$ stochastic estimate of the trace of a Padé approximation to the logarithm of the fermion matrix. The scheme is still being developed. One interesting feature  is the possibility of using the algorithm to simulate finite density QCD, by including a projection onto states with definite particle number.

Bakeyev  has proposed a new scheme for performing exact, finite updates of the gauge fields interacting with dynamical fermions. The update involves finding a new configuration, $`\stackrel{~}{U}`$ that solves

$$M(\stackrel{~}{U})\eta =\omega M(U)\eta $$
(16)
for a stochastic parameter, $`\omega `$ and gaussian noise source, $`\eta `$. The idea is in its infancy as yet, and may prove to be most useful as an additional step in an HMC simulation.

The flavour-symmetry breaking of staggered fermions can be reduced by using a “fat link” to induce quark-gluon interactions. Simulating these staggered actions poses a new problem in algorithm design. The fat link is constructed following a similar procedure to APE smearing, which involves summing the six staples around a link, then projecting the resulting matrix back into the gauge group. The fattening step is repeated iteratively. The challenge for simulation is that the force term on the underlying gauge fields is difficult to compute. New algorithms are being developed and tested to circumvent this problem. Knechtli  presented these ideas to the conference. HMC fat-link simulations with four flavours when the lightest pion has mass $`r_0m_\pi =2.0`$ are seven times more costly than runs where the quarks interact directly with the link variables. This factor is reduced to about 3 once $`r_0m_\pi `$ reaches 1.6 This extra overhead is more than offset by the much better flavour symmetry of the fat-link actions, and the quark mass accessible to these studies should be much lighter.

## 6 HIERARCHICAL NOISE REDUCTION IN YANG-MILLS THEORIES

Recently, Lüscher and Weisz  proposed a novel algorithm for reducing the stochastic noise in large Wilson loop expectation values by a factor that grows exponentially with the loop size. The scheme involves building a hierarchy of sub-domains within the lattice, and performing a reduced Monte-Carlo integration over the degrees of freedom inside each domain. The results at each level are passed up the hierarchy to construct variance reduced averages of larger and larger objects. The scheme is a natural extension of the temporal link “multi-hit” method . They illustrate their proposal with a calculation of the correlator of two spatially-separated Polyakov lines (see Fig. 3), and show a reduction of many orders of magnitude in the statistical error over the standard technique.

The scheme can be extended to include matter fields. If the fermion determinant is included using the LBA, the action is localised within a small neighbourhood. With Wilson fermions (and the Wilson parameter set to unity) the hierarchical measurement scheme can be incorporated without modification. The difficulty is including the stochastic accept-reject schemes usually employed with the LBA to make the method exact. These rely on lattice-wide observations, involving all the gauge degrees of freedom and the domain decompositions of Ref  would no longer apply.

It remains an unsolved challenge to include the highly developed tools for optimising ground-state overlap into the hierarchy scheme. Methods such as APE smearing and Teper blocking have proved crucial in the accurate calculations of the inter-quark potential and its excitations in Yang-Mills theory. Variational techniques, which build a large basis sets of creation operators for the states of interest, then find an optimal ground-state operator by diagonalisation have also dramatically improved these calculations (for an example, see Ref. ). While the theoretical aspect of including these ideas is readily resolved, it is not obvious how well this toolkit can be used within the hierarchy method in practice. An empirical test seems the only way to resolve this issue. If the noise reduction of the method of Ref.  can be married to the excellent ground-state construction of the smoothing and operator basis methods, then extremely precise calculations of the detailed nature of the confining string spectrum at large separations could be performed.

## 7 ALL-TO-ALL FERMION PROPAGATORS

Having expended such a considerable number of cycles on some of the world’s largest supercomputers in generating an ensemble of gauge field configurations, it is clear we are duty bound to make the most of the information they contain.

For many observables, believed to be sensitive to vacuum fluctuations in the quark fields, traditional point-source propagator methods are inefficient. To extract a useful signal from the ensemble, the propagator from all (or many) points to all points on the lattice is needed. There has been a good deal of interest in improving these techniques in recent years.

### 7.1 Gaussian and $`Z_2`$ stochastic estimators

A gaussian representation of a fermion matrix entry can be written :

$$Q_{ij}=𝒟\varphi 𝒟\varphi ^{}\varphi _i(\varphi ^{}Q)_j\mathrm{exp}\left\{\varphi ^{}Q^2\varphi \right\}.$$
(17)
Thus a propagator from any point on the lattice to any other can be estimated by performing a sub-Monte Carlo simulation on each configuration.

At the conference, Duncan  presented simulation data using a gaussian all-to-all scheme. He emphasised that for many interesting physics applications, this simplest of methods is adequate to ensure the dominant source of statistical error is from fluctuations within the gauge ensemble, rather than the stochastic estimator. The stochastic degrees of freedom can be updated in a global heatbath step, but this requires a matrix inversion. Ref.  gives a scheme to accelerate the process; the iterative matrix solver is stopped after a small number of steps, and a Metropolis test is used to correct for the error in this step.

The Kentucky group  noted that other stochastic variables can be used to compute matrix traces, and exploited random elements of $`Z_2`$. These estimators have been employed in a range of calculations involving disconnected diagrams.

### 7.2 Improving stochastic estimators

Ref.  suggests reducing the noise in a gaussian estimator by dividing the lattice into two segments, then computing the independent sets of integrals in each segment exactly in a fixed background of the degrees of freedom on the interface between the segments. Conceptually, the new Lüscher-Weisz hierarchy method, described in Sec. 6 is reminiscent of this scheme. The method has one restriction; the source and sink points of the propagator can not be in the same domain, so closed fermion loops can not be computed. In Ref.  a further enhancement of the method is presented.

At the conference, a method was presented by Wilcox  to mix gaussian and $`Z_2`$ noise, allowing an approximate heatbath with Metropolis correction to be used.

The SESAM collaboration  improved $`Z_2`$ estimators by breaking the vector space into sub-spaces (in their study, the 4 sub-spaces corresponding to the spin indices) and performing separate stochastic estimators in each sub-space. The idea was tested further in Ref . While this means more matrix inversion must be performed for a fixed stochastic ensemble size, the resulting estimators can have a significantly smaller variance. This idea of thinning can be extended beyond just spin and colour indices to spatial sites as well. If the scheme were to be carried to the extreme, where the $`N`$ dimensional vector space is decomposed into $`N`$ 1-dimensional spaces, a $`Z_2`$ estimator would reproduce the exact result for any ensemble size. This is not such a surprising result, since this computation would require $`N`$ inversions to be performed and it is a trivial result that the trace of $`M^1`$ can be computed in $`N`$ inversions. Without thinning however the error falls off, as any statistical estimate should, like $`1/\sqrt{N}`$ so eventually thinning must become more efficient.

The thinning procedure has not been explored fully for computations of all-to-all propagators in QCD, and it would be worthwhile to investigate this more completely.

### 7.3 Eigenvalue decomposition

If all the eigenvectors (and eigenvalues) of the fermion matrix are known, then any propagator entry can be computed straightforwardly. It is natural to expect that the long-range physics of QCD are contained in the lowest-lying eigenvectors, and these will then dominate in the spectral representation of a mesonic correlator (for example). The importance of understanding the physics of the lowest modes of the fermion matrix was emphasised in the plenary talk of Edwards .

Neff presented the conference  with an investigation of methods for computing the low-lying eigenvectors efficiently. The acceleration comes from using a polynomial to isolate the region of eigenvalues of interest, and this leads to more rapid converges of the Arnoldi method. At the conference, Schilling  discussed the dominance of the $`\pi `$ and $`\eta `$ correlators by the low-lying eigenmodes. Truncating the spectral expansion introduces a systematic error, but this can be avoided by computing all-to-all propagators in a hybrid scheme (as noted by Edwards ).

### 7.4 Hybrid schemes

The two schemes can be combined straightforwardly to construct a hybrid all-to-all estimator method . A number of low-lying eigenvalues of the fermion matrix are computed and incorporated exactly and the effect of all the remaining eigenvalues is estimated stochastically. This involves building a projection operator onto the reduced sub-space spanned by the explicitly computed eigenvectors.

For illustration, consider the observable

$$\text{Tr }Q^1\mathrm{\Gamma }=\underset{i}{\overset{N}{}}\frac{1}{\lambda _i}v_i^{}\mathrm{\Gamma }v_i.$$
(18)
If the lowest $`m`$ ($`mN`$) eigenvectors of the hermitian matrix $`Q=\gamma _5M`$ are computed, then it is natural break the the spectral representation of $`Q`$ into two parts,

$$Q=\underset{i=1}{\overset{m}{}}\lambda _iv_iv_i^{}+\underset{i=m+1}{\overset{N}{}}\lambda _iv_iv_i^{}.$$
(19)
Defining the (sub-space) inverses and projectors,

$$\overline{Q}_{(0)}=\underset{i=1}{\overset{m}{}}\frac{1}{\lambda _i}v_iv_i^{}\text{}𝒫_{(0)}=\underset{i=1}{\overset{m}{}}v_iv_i^{}$$
(20)
the trace becomes

$$\text{Tr }Q^1\mathrm{\Gamma }=\text{Tr }\overline{Q}_{(0)}\mathrm{\Gamma }+\text{Tr }\overline{Q}_{(1)}\mathrm{\Gamma }.$$

The first term is the truncated estimate, and the remaining term can be estimated stochastically;

$`\text{Tr }\overline{Q}_{(1)}\mathrm{\Gamma }`$ $`=`$ $`\eta ^{}\mathrm{\Gamma }\overline{Q}_{(1)}\eta `$ (21)
$`=`$ $`\eta ^{}\mathrm{\Gamma }(Q^1\overline{Q}_{(0)})\eta `$
$`=`$ $`\eta ^{}\mathrm{\Gamma }Q^1(I𝒫_{(0)})\eta .`$
The operation of $`I𝒫_{(0)}`$ is an orthogonalisation of the noise vector, $`\eta `$ with respect to all $`m`$ known eigenvectors. The matrix inversion in this last step is accelerated, since the eigenvectors corresponding to the low-lying eigenvalues have been removed from $`\eta `$. It will be interesting to see hybrid schemes tested in practical applications. The developments of GMRES algorithms presented to the conference  look particularly interesting in this context.

## 8 CONCLUSIONS

For the next few years, most large-scale simulations of full QCD will use the Wilson or staggered fermion formulation, while the chiral lattice actions  are beyond the reach of current computing resources; Ginsparg-Wilson fermions are a factor of 100 times more expensive. As pointed out by Karl Jansen, the Hybrid Monte-Carlo algorithm is not far off its twentieth birthday and remains (largely unmodified) the most popular tool for dynamical Wilson fermion simulations.

Current predictions  are that resources available in the next few years are an order of magnitude below the requirements for reliable QCD simulations. Including the newest ideas presented, it might be hoped to get close to the required factor of ten improvement in performance.

Realistic comparative tests of HMC vs LBA have been carried out in recent years, leading to the observation that the two methods have identical costs (for all practical purposes). These comparisons are always very difficult to make. It would be helpful to see a standard benchmark emerge on which direct comparisons of algorithms can be made. Defining a standard is difficult but with this in place, fair races of new methods against “thoroughbred” algorithms rather than the simplest implementations could be made.

Excitement at the conference generated by the recent idea of Lüscher and Weisz suggests reconsidering the traditional Monte-Carlo analysis pathway; generating configurations then subsequently making measurements. Their hierarchy scheme, which knits these two processes together can lead to orders of magnitude increases in efficiency. The challenge is to apply this philosophy to a wider range of applications.

I am grateful to Ph. de Forcrand, C. Gebert, M. Hasenbusch, B. Joó, A. Kennedy, J. Kuti, M. Lüscher, I. Montvay, H. Neff, W. Schroers, J. Sexton and N. Zverev for many stimulating discussions, observations and helpful correspondence.


text:
# Quantum Zeno subspaces

## Abstract

The quantum Zeno effect is recast in terms of an adiabatic theorem when the measurement is described as the dynamical coupling to another quantum system that plays the role of apparatus. A few significant examples are proposed and their practical relevance discussed. We also focus on decoherence-free subspaces.

If very frequent measurements are performed on a quantum system, in order to ascertain whether it is still in its initial state, transitions to other states are hindered and the quantum Zeno effect takes place Beskow ; Misra . This phenomenon stems from general features of the Schrödinger equation that yield quadratic behavior of the survival probability at short times strev ; PIO . The first realistic test of the quantum Zeno effect (QZE) for oscillating (two-level) systems was proposed about 15 years ago Cook . This led to experiments, discussions and new proposals Itano . A few years ago, the presence of a short-time quadratic region was experimentally confirmed also for a bona fide unstable system Wilkinson . The same experimental setup has been used very recently raizenlatest  in order to prove the existence of the Zeno effect (as well as its inverse antiZeno ; heraclitus ) for an unstable quantum mechanical system, leading to new ideas Agarwal01 ; Frishman01 .

It is important to stress that the quantum Zeno effect does not necessarily freeze everything. On the contrary, for frequent projections onto a multi-dimensional subspace, the system can evolve away from its initial state, although it remains in the subspace defined by the measurement. This continuing time evolution within the projected subspace (“quantum Zeno dynamics”) has been recently investigated compactregularize . It has peculiar physical and mathematical features and sheds light on some subtle mathematical issues Friedman72 ; Gustafson ; Misra .

All the above-mentioned investigations deal with what can be called “pulsed” measurements, according to von Neumann’s projection postulate von . However, from a physical point of view, a “measurement” is nothing but an interaction with an external system (another quantum object, or a field, or simply another degree of freedom of the very system investigated), playing the role of apparatus. In this respect, if one is not too demanding in philosophical terms, von Neumann’s postulate can be regarded as a useful shorthand notation, summarizing the final effect of the quantum measurement. This simple observation enables one to reformulate the QZE in terms of a (strong) coupling to an external agent. We emphasize that in such a case the QZE is a consequence of the dynamical features (i.e. the form factors) of the coupling between the system investigated and the external system, and no use is made of projection operators (and non-unitary dynamics). The idea of “continuous” measurement in a QZE context has been proposed several times during the last two decades Peres80 ; PeresKraus , although the first quantitative comparison with the “pulsed” situation is rather recent Schulman98 .

The purpose of the present article is to cast the quantum Zeno evolution in terms of an adiabatic theorem and study possible applications. We will see that the evolution of a quantum system is profoundly modified (and can be tailored in an interesting way) by a continuous measurement process: the system is forced to evolve in a set of orthogonal subspaces of the total Hilbert space and a dynamical superselection rule arises in the strong coupling limit. These general ideas will be corroborated by some simple examples.

We start by considering the case of “pulsed” observation. We first extend Misra and Sudarshan’s theorem Misra  in order to accomodate multiple projectors. Let Q be a quantum system, whose states belong to the Hilbert space $``$ and whose evolution is described by the unitary operator $`U(t)=\mathrm{exp}(iHt)`$, where $`H`$ is a time-independent lower-bounded Hamiltonian. Let $`\{P_n\}_n(P_nP_m=\delta _{mn}P_n,_nP_n=1)`$ be a (countable) collection of projection operators and $`\text{Ran}P_n=_{P_n}`$ the relative subspaces. This induces a partition on the total Hilbert space $`=_n_{P_n}`$. Let $`\rho _0`$ be the initial density matrix of the system. We “prepare” the system by performing an initial measurement, described by the superoperator

$$\widehat{P}\rho =\underset{n}{}P_n\rho P_n=\rho _0.$$
(1)
The free evolution reads

$$\widehat{U}_t\rho _0=U(t)\rho _0U^{}(t),U(t)=\mathrm{exp}(iHt)$$
(2)
and the Zeno evolution after $`N`$ measurements in a time $`t`$ is governed by the superoperator

$$\widehat{V}_t^{(N)}=\widehat{P}\left(\widehat{U}\left(t/N\right)\widehat{P}\right)^{N1},$$
(3)
which yields

$$\rho (t)=\widehat{V}_t^{(N)}\rho _0=\underset{n_1,\mathrm{},n_N}{}V_{n_1\mathrm{}n_N}^{(N)}(t)\rho _0V_{n_1\mathrm{}n_N}^{(N)}(t),$$
(4)
where

$`V_{n_1\mathrm{}n_N}^{(N)}(t)=P_{n_N}U\left(t/N\right)P_{n_{N1}}\mathrm{}P_{n_2}U\left(t/N\right)P_{n_1}.`$ (5)
We follow Misra  and assume for each $`n`$ the existence of the strong limits ($`t>0`$)

$$\underset{N\mathrm{}}{lim}V_{n\mathrm{}n}^{(N)}(t)𝒱_n(t),\underset{t0^+}{lim}𝒱_n(t)=P_n.$$
(6)
Then $`𝒱_n(t)`$ exist for all real $`t`$ and form a semigroup Misra , and $`𝒱_n^{}(t)𝒱_n(t)=P_n`$. Moreover it is easy to show that

$$\underset{N\mathrm{}}{lim}V_{n\mathrm{}n^{}\mathrm{}}^{(N)}(t)=0,\text{for}n^{}n.$$
(7)
Therefore the final state is

$`\rho (t)=\widehat{𝒱}_t\rho _0={\displaystyle \underset{n}{}}𝒱_n(t)\rho _0𝒱_n^{}(t),`$ (8)
$`\text{with}{\displaystyle \underset{n}{}}𝒱_n^{}(t)𝒱_n(t)={\displaystyle \underset{n}{}}P_n=1.`$
The components $`𝒱_n(t)\rho _0𝒱_n^{}(t)`$ make up a block diagonal matrix: the initial density matrix is reduced to a mixture and any interference between different subspaces $`_{P_n}`$ is destroyed (complete decoherence). In conclusion,

$$p_n(t)=\text{Tr}\left(\rho (t)P_n\right)=\text{Tr}\left(\rho _0P_n\right)=p_n(0),n.$$
(9)
In words, probability is conserved in each subspace and no probability “leakage” between any two subspaces is possible. The total Hilbert space splits into invariant subspaces and the different components of the wave function (or density matrix) evolve independently within each sector. One can think of the total Hilbert space as the shell of a tortoise, each invariant subspace being one of the scales. Motion among different scales is impossible. (See Fig. 1 in the following.) The study of the Zeno dynamics within a given infinite-dimensional subspace is an interesting problem compactregularize  that will not be discussed here. The original formulation of the Zeno effect is reobtained when $`p_n=1`$ for some $`n`$, in (9): the initial state is then in one of the invariant subspaces and the survival probability in that subspace remains unity.

The previous theorem hinges upon von Neumann’s projections von . However, as we explained in the introduction, a QZE can also be obtained by performing a continuous measurement on a system. For example, consider the Hamiltonian

$$H_{\text{3lev}}=\mathrm{\Omega }\sigma _1+K\tau _1=\left(\begin{array}{ccc}0& \mathrm{\Omega }& 0\\ \mathrm{\Omega }& 0& K\\ 0& K& 0\end{array}\right),$$
(10)
describing two levels (system), with Hamiltonian $`H=\mathrm{\Omega }\sigma _1=\mathrm{\Omega }(|12|+|21|)`$, coupled to a third one, that plays the role of measuring apparatus: $`KH_{\text{ meas}}=K\tau _1=K(|23|+|32|)`$. This model, first considered in Peres80 , is probably the simplest way to include an “external” apparatus in our description: as soon as the system is in $`|2`$ it undergoes Rabi oscillations to $`|3`$. We expect level $`|3`$ to perform better as a measuring apparatus when the strength $`K`$ of the coupling becomes larger. Indeed, if initially the system is in state $`|1`$, the survival probability reads

$$p(t)=\left[K^2+\mathrm{\Omega }^2\mathrm{cos}(K_1t)\right]^2/K_1^4\stackrel{K\mathrm{}}{}1,$$
(11)
where $`K_1=\sqrt{K^2+\mathrm{\Omega }^2}`$. This simple model captures many interesting features of a Zeno dynamics (and will help clarify our general approach). Many similar examples can be considered: in general PIO ; Napoli , one can include the detector in the quantum description, by considering the Hamiltonian

$$H_K=H+KH_{\text{meas}},$$
(12)
where $`H`$ is the Hamiltonian of the system under observation (and can include the free Hamiltonian of the apparatus) and $`H_{\text{meas}}`$ is the interaction Hamiltonian between the system and the apparatus.

We now prove a theorem, which is the exact analog of Misra and Sudarshan’s theorem for a dynamical evolution of the type (12). Consider the time evolution operator

$`U_K(t)=\mathrm{exp}(iH_Kt).`$ (13)
We will prove that in the “infinitely strong measurement” limit $`K\mathrm{}`$ the evolution operator

$$𝒰(t)=\underset{K\mathrm{}}{lim}U_K(t),$$
(14)
becomes diagonal with respect to $`H_{\text{meas}}`$:

$$[𝒰(t),P_n]=0,\text{where}H_{\text{meas}}P_n=\eta _nP_n,$$
(15)
$`P_n`$ being the orthogonal projection onto $`_{P_n}`$, the eigenspace of $`H_{\text{meas}}`$ belonging to the eigenvalue $`\eta _n`$. Note that in Eq. (15) one has to consider distinct eigenvalues, i.e., $`\eta _n\eta _m`$ for $`nm`$, whence the $`_{P_n}`$’s are in general multidimensional.

The theorem is easily proven by recasting it in the form of an adiabatic theorem. In the $`H`$ interaction picture,

$$H_{\text{meas}}^\text{I}(t)=e^{iHt}H_{\text{meas}}e^{iHt},$$
(16)
the Schrödinger equation reads

$$i_tU_K^\text{I}(t)=KH_{\text{meas}}^\text{I}(t)U_K^\text{I}(t).$$
(17)
This has exactly the same form of an adiabatic evolution $`i_sU_T(s)=TH(s)U_T(s)`$ Messiah61 : the large coupling $`K`$ limit corresponds to the large time $`T`$ limit and the physical time $`t`$ to the scaled time $`s=t/T`$. In the $`K\mathrm{}`$ limit, by considering a spectral projection $`P_n^\text{I}(t)=e^{iHt}P_ne^{iHt}`$ of $`H_{\text{meas}}^\text{I}(t)`$, the limiting operator $`𝒰^\text{I}(t)=lim_K\mathrm{}U_K^\text{I}(t)`$ satisfies the intertwining property $`𝒰^\text{I}(t)P_n^\text{I}(0)=P_n^\text{I}(t)𝒰^\text{I}(t)`$, i.e. maps $`_{P_n^\text{I}(0)}`$ onto $`_{P_n^\text{I}(t)}`$:

$$\psi _0^\text{I}_{P_n^\text{I}(0)}\psi ^\text{I}(t)_{P_n^\text{I}(t)}.$$
(18)
In the Schrödinger picture

$$\psi _0_{P_n}\psi (t)_{P_n},$$
(19)
whence

$`\rho (t)=e^{iHt}𝒰^\text{I}(t)\rho _0𝒰^\text{I}(t)e^{iHt}=𝒰(t)\rho _0𝒰^{}(t),`$ (20)
where $`𝒰(t)`$ has the property (15) and the probability to find the system in $`_{P_n}`$ satisfies Eq. (9) and is therefore constant: if the initial state of the system belongs to a given sector, it will be forced to remain there forever (QZE).

Even more, by exploiting the features of the adiabatic theorem in greater details, it is possible to show that, for time independent Hamiltonians, the limiting evolution operator has the explicit form tokyo

$$𝒰(t)=\mathrm{exp}[i(H_{\text{diag}}+KH_{\text{meas}})t],$$
(21)
where

$$H_{\text{diag}}=\underset{n}{}P_nHP_n$$
(22)
is the diagonal part of the system Hamiltonian $`H`$ with respect to the interaction Hamiltonian $`H_{\text{meas}}`$.

Let us briefly comment on the physical meaning. In the $`K\mathrm{}`$ limit, due to (15), the time evolution operator becomes diagonal with respect to $`H_{\text{meas}}`$, $`[𝒰(t),H_{\text{meas}}]=0`$, an effective superselection rule arises and the total Hilbert space is split into subspaces $`_{P_n}`$ that are invariant under the evolution. These subspaces are defined by the $`P_n`$’s, i.e., they are eigenspaces belonging to distinct eigenvalues $`\eta _n`$: in other words, subspaces that the apparatus is able to distinguish. On the other hand, due to (22), the dynamics within each Zeno subspace $`_{P_n}`$ is governed by the diagonal part $`P_nHP_n`$ of the system Hamiltonian $`H`$. This bridges the gap with the description (1)-(9) and clarifies the role of the detection apparatus. In Fig. 1 we endeavored to give a pictorial representation of the decomposition of the Hilbert space as $`K`$ is increased. It is worth noticing that the superselection rules discussed here are de facto equivalent to the celebrated “W3” ones WWW , but turn out to be a mere consequence of the Zeno dynamics.

Four examples will prove useful. First example: reconsider $`H_{\text{3lev}}`$ in Eq. (10). As $`K`$ is increased, the Hilbert space is split into three invariant subspaces (the three eigenspaces of $`H_{\text{meas}}=\tau _1`$): (level $`|1`$) $``$ (level $`|2+|3`$) $``$ (level $`|2|3`$).

Second example: consider

$$H_{\text{4lev}}=\left(\begin{array}{cccc}0& \mathrm{\Omega }& 0& 0\\ \mathrm{\Omega }& 0& K& 0\\ 0& K& 0& K^{}\\ 0& 0& K^{}& 0\end{array}\right),$$
(23)
where level $`|4`$ “measures” whether level $`|3`$ is populated. If $`K^{}K\mathrm{\Omega }`$, the total Hilbert space is divided into three subspaces: (levels $`|1`$ and $`|2`$) $``$ (level $`|3+|4`$) $``$ (level $`|3|4`$). Notice that the $`\mathrm{\Omega }`$ oscillations are restored as $`K^{}K`$ (in spite of $`K\mathrm{\Omega }`$). A watched cook can freely watch a boiling pot.

Third example (decoherent-free subspaces Viola99  in quantum computation). The Hamiltonian Beige00

$$H_{\text{meas}}=ig\underset{i=1}{\overset{2}{}}[b|2_{ii}1|\text{H.c.}]i\kappa b^{}b$$
(24)
describes a system of two ($`i=1,2`$) three-level atoms in a cavity. The atoms are in a $`\mathrm{\Lambda }`$ configuration with split ground states $`|0_i`$ and $`|1_i`$ and excited state $`|2_i`$, while the cavity has a single resonator mode $`b`$ in resonance with the atomic transition 1-2. Spontaneous emission inside the cavity is neglected, but a photon leaks out through the nonideal mirrors with a rate $`\kappa `$. The (5-dimensional) eigenspace $`_{P_0}`$ of $`H_{\text{meas}}`$ belonging to the eigenvalue $`\eta =0`$ is spanned by

$$\{|000,|001,|010,|011,(|021|012)/\sqrt{2}\},$$
(25)
where $`|0j_1j_2`$ denotes a state with no photons in the cavity and the atoms in state $`|j_1_1|j_2_2`$. If the coupling $`g`$ and the cavity loss $`\kappa `$ are sufficiently strong, any other weak Hamiltonian $`H`$ added to (24) reduces to $`P_0HP_0`$ and changes the state of the system only within the decoherence-free subspace (25).

Fourth example. Let

$$H=\left(\begin{array}{ccc}0& \tau _\text{Z}^1& 0\\ \tau _\text{Z}^1& i2/\tau _\text{Z}^2\gamma & K\\ 0& K& 0\end{array}\right).$$
(26)
This describes the spontaneous emission $`|1|2`$ of a system into a (structured) continuum, while level $`|2`$ is resonantly coupled to a third level $`|3`$ PIO . This case is also relevant for quantum computation, if one is interested in protecting a given subspace (level $`|1`$) from decoherence Viola99 ; Beige00  by inhibiting spontaneous emission Agarwal01 . Here $`\gamma `$ represents the decay rate to the continuum and $`\tau _\text{Z}`$ is the Zeno time (convexity of the initial quadratic region). As the Rabi frequency $`K`$ is increased one is able to hinder spontaneous emission from level $`|1`$ (to be protected) to level $`|2`$. However, in order to get an effective “protection” of level $`|1`$, one needs $`K>1/\tau _\text{Z}`$. More to this, when the presence of the inverse Zeno effect is taken into account, this requirement becomes even more stringent heraclitus  and yields $`K>1/\tau _\text{Z}^2\gamma `$. Both these conditions can be very demanding for a real system subject to dissipation PIO ; Napoli ; heraclitus . For instance, typical values for spontaneous decay in vacuum are $`\gamma 10^9`$s-1, $`\tau _\text{Z}^210^{29}`$s2 and $`1/\tau _\text{Z}^2\gamma 10^{20}`$s-1 hydrogen .

The formulation of a Zeno dynamics in terms of an adiabatic theorem is powerful. Indeed one can use all the machinery of adiabatic theorems in order to get results in this context. An interesting extension would be to consider time-dependent measurements

$$H_{\text{meas}}=H_{\text{meas}}(t),$$
(27)
whose spectral projections $`P_n=P_n(t)`$ have a nontrivial time evolution. In this case, instead of confining the quantum state to a fixed sector, one can transport it along a given path (subspace) $`_{P_n(t)}`$. One then obtains a dynamical generalization of the process pioneered by Von Neumann in terms of projection operators von ; AV . The influence of non-adiabatic corrections, for $`K`$ large but finite, as well as practical estimates and protection from decoherence effects will be considered in a future article.

###### Acknowledgements.

We thank L. Neglia and N. Cillo for the drawing.


text:
# Theory of Systematic Computational Error in Free Energy Differences

## Abstract

Systematic inaccuracy is inherent in any computational estimate of a non-linear average, due to the availability of only a finite number of data values, $`N`$. Free energy differences $`\mathrm{\Delta }F`$ between two states or systems are critically important examples of such averages in physical, chemical and biological settings. Previous work has demonstrated, empirically, that the “finite-sampling error” can be very large — many times $`k_BT`$ — in $`\mathrm{\Delta }F`$ estimates for simple molecular systems. Here, we present a theoretical description of the inaccuracy, including the exact solution of a sample problem, the precise asymptotic behavior in terms of $`1/N`$ for large $`N`$, the identification of universal law, and numerical illustrations. The theory relies on corrections to the central and other limit theorems, and thus a role is played by stable (Lévy) probability distributions.

Introduction. Free energy difference calculations have a tremendous range of applications in physical, chemical, and biological systems; examples include computations relating magnetic phases, estimates of chemical potentials, and of binding affinities of ligands to proteins (e.g., ). Since the work of Kirkwood , it has been appreciated that the free energy difference, $`\mathrm{\Delta }F\mathrm{\Delta }F_{01}`$, of switching from a Hamiltonian $`_0`$ to $`_1`$ is given by a *non-linear average*,

$$\mathrm{\Delta }F=k_BT\mathrm{log}\left[\mathrm{exp}(W_{01}/k_BT)_0\right],$$
(1)
where $`k_BT`$ is the thermal unit of energy at temperature $`T`$ and $`W_{01}`$ is the work required to switch the system from $`_0`$ to $`_1`$. The angled brackets indicate an average over switches starting from configurations drawn from the equilibrium distribution governed by $`_0`$. In instantaneous switching the work is defined by $`W_{01}=_1(\text{x})_0(\text{x})`$ for a start (and end) configuration x; however, gradual switches requiring a “trajectory”-based work definition may also be used as was pointed out by Jarzynski .

Whenever a convex, nonlinear average such as (1) is estimated computationally, that result will *always* be systematically biased  because one has only a finite amount of data — say, $`N`$ work values. The bias results from incomplete sampling of the smallest (or most negative) $`W_{01}`$ values: these values dominate the average (1) and cannot be sampled perfectly for finite $`N`$, regardless of the $`W_{01}`$ distribution. Thus, a running estimate of $`\mathrm{\Delta }F`$ will typically decline as data is gathered. Such considerations led Wood *et al.*  to consider the block-averaged $`n`$-data-point estimate of the free energy based on $`N=mn`$ total work values $`\{W^{(k)}\}`$, namely,

$$\mathrm{\Delta }F_n=\frac{1}{m}\underset{j=1}{\overset{m}{}}k_BT\mathrm{log}\left[\frac{1}{n}\underset{k=(j1)n+1}{\overset{jn}{}}\mathrm{exp}(W^{(k)}/k_BT)\right].$$
(2)
In the limit $`m,N\mathrm{}`$, $`\mathrm{\Delta }F_n`$ is mathematically well-defined and amenable to analysis; it represents the expected value of a free energy estimate from $`n`$ data points — that is, of

$$_n=k_BT\mathrm{log}\left[\left(e^{W_1/k_BT}+\mathrm{}+e^{W_n/k_BT}\right)/n\right].$$
(3)
See Fig. 1. Wood *et al.* estimated the lowest order correction to $`\mathrm{\Delta }F\mathrm{\Delta }F_{\mathrm{}}`$ as $`\sigma _w^2/2nk_BT`$, where $`\sigma _w^2`$ is the variance in the distribution of work values, $`W`$ . Ferrenberg, Landau and Binder discussed analogous issues for the magnetic susceptibility .

More recently, Zuckerman and Woolf  suggested a means by which a range of $`\mathrm{\Delta }F_n`$ values for $`n<N`$ could be used to extrapolate to the true, infinite-data answer, $`\mathrm{\Delta }F`$. The authors also observed that, for large $`m=N/n`$, the free energy is bounded according to

$$\mathrm{\Delta }F\mathrm{\Delta }F_n,\text{any }n.$$
(4)
This inequality results from the convexity of the exponential function, as will be demonstrated explicitly in a fuller account of the theory. Finally, Zuckerman and Woolf noted that the leading behavior of $`\mathrm{\Delta }F_n`$ appeared to be *not always linear* in $`1/n`$ but, rather, seemed to behave as $`(1/n)^{\tau _1}`$ for $`\tau _11`$.

This Letter presents the theory — apparently for the first time — describing the finite-sampling inaccuracy for $`\mathrm{\Delta }F`$ estimates. Previous work discussing $`\mathrm{\Delta }F_n`$ has been, primarily, empirical . Our report includes (i) the formal analytic expression for the expected value of the error from $`N`$ work values, $`\mathrm{\Delta }F_N\mathrm{\Delta }F`$; (ii) an exact solution, for all $`N`$, of this expected value when the Boltzmann factor of the work value $`ze^{W/k_BT}`$ follows a gamma distribution; (iii) exact asymptotic expressions for $`\mathrm{\Delta }F_n`$ and the variance of $``$ as $`n\mathrm{}`$ for arbitrary $`W`$ distribtions, including non-analytic behavior in the case when the variance and higher moments of $`z`$ diverge; and (iv) discussion and numerical illustrations based on Gaussian distributions of $`W`$, plus corrections expected from skewed Gaussian distributions. The present discussion makes use of mathematical results regarding the convergence — to “stable” limiting distributions , also known as Lévy processes (e.g., ) — of the distributions of sums of variables. The results are expected to have practical application in the extrapolation process outlined in .

Formal Development of $`\mathrm{\Delta }F_n`$. The derivation proceeds via continuum expressions simplified by the definitions $`wW/k_BT`$, $`f\mathrm{\Delta }F/k_BT`$, and $`f_n\mathrm{\Delta }F_n/k_BT`$. First, in terms of the probability density $`\rho _w`$ of work values, which is normalized by $`dw\rho _w(w)=1`$, the free energy is given by the continuum analog of (1),

$$f=\mathrm{\Delta }F/k_BT=\mathrm{log}\left[dw\rho _w(w)e^w\right].$$
(5)
The finite-data average free energy, following (2) must apply the logarithm “before” the average of the $`n`$ Boltzmann factors, and one has

$`f_n`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}\left[\mathrm{d}w_i\rho _w(w_i)\right]\mathrm{log}\left[\frac{1}{n}\underset{i=1}{\overset{n}{}}e^{w_i}\right]},`$ (6)
Now, motivated by the central and related limit theorems  for the sum of the $`e^w`$ variables, we introduce a change of variables which will permit the development of a $`1/n`$ expansion for $`f_n`$. In particular, we define

$$y=(e^{w_1}+\mathrm{}+e^{w_n}ne^f)/b_1n^{1/\alpha },$$
(7)
where $`b_1`$ is a constant and $`\alpha 2`$ is an exponent characterizing the distribution of the variable $`e^w`$. In fact, the requirement that $`\mathrm{\Delta }F`$ be finite in (5) further implies $`\alpha >1`$. The finite-data free energy difference can now be written

$$f_n=_{cn^a}^{\mathrm{}}dy\rho _n(y)\mathrm{log}\left(e^f+\frac{b_1}{n^a}y\right)$$
(8)
where $`c=\mathrm{exp}(f)/b_1`$, $`a(\alpha 1)/\alpha <1/2`$, and $`\rho _n`$ is the probability density of the variable $`y`$. Note that $`a`$ is always positive because $`\alpha >1`$.

To continue, we must call upon some mathematical results regarding the approach, with increasing $`n`$, to general stable limit distributions (of which the Gaussian, for $`\alpha =2`$, is the best known ). More precisely, the sum of *any* set of random variables, suitably normalized as in (7), has a distribution with zero mean which may be expressed as a stable distribution function multiplied by a large-$`n`$ asymptotic expansion .

Finite-Moments Case and An Exact Solution. To illustrate the case of a Gaussian limit ($`\alpha =2`$), assume the variable $`e^w`$ possesses finite “Boltzmann moments” — a mean $`\widehat{\mu }=e^f`$, variance $`\widehat{\sigma }^2`$, and third moment $`\widehat{\mu }_3`$ — not to be confused with the moments of the distribution of $`w`$. The finite-$`n`$ corrections to the central limit theorem indicate that the variable $`y=(^ne^{w_i}n\widehat{\mu })/\sqrt{n}\widehat{\sigma }`$ \[cf. (7)\] is distributed according to

$$\rho _n(y)=\rho _G(u;1)\left[1+\nu _1(y)/\sqrt{n}+\nu _2(y)/n+\mathrm{}\right],$$
(9)
for large $`n`$, where the remaining terms are higher integer powers of $`1/\sqrt{n}`$ and the Gaussian density is

$$\rho _G(y;\sigma )=\mathrm{exp}(y^2/2\sigma ^2)/\sqrt{2\pi }\sigma ,$$
(10)
The $`\nu _i`$ depend on the original distribution of $`e^w`$; for instance, $`\nu _1(y)=(\widehat{\mu }_3/6\widehat{\sigma }^3)(y^33y)`$ . Moreover, the $`\nu `$ functions are odd or even according to whether $`i`$ is odd or even, in this $`\alpha =2`$ case.

One arrives at the explicit form of the finite-data-corrected free energy for the case of finite $`\widehat{\sigma }^2`$ and $`\widehat{\mu }_3`$ by substituting (9) into (8), along with an expansion of the logarithm about $`y=0`$. (More careful consideration of series convergence for large $`y`$ yields the same final result for $`f_n`$, as will be elucidated in future work.) Because of the odd- and even-ness of the factors to be integrated, one finds an expansion consisting *solely of integer powers of $`1/n`$*, namely,

$$f_n=f+\phi _1/n+\phi _2/n^2+\mathrm{},$$
(11)
with $`\phi _1=\widehat{\sigma }^2/2\widehat{\mu }^2`$ and $`\phi _2=(4\widehat{\mu }\widehat{\mu }_39\widehat{\sigma }^4)/12\widehat{\mu }^4`$ . To compare this with the finding of Wood *et al.* for $`f_nf`$, one can consider a Gaussian distribution of $`W=k_BTw`$ with variance $`\sigma _w^2`$: expanding the resulting Boltzmann moments of previous result for small $`\sigma _w`$ yields $`\phi _1=k_BT[\mathrm{exp}[(\sigma _w/k_BT)^2]1]/2\sigma _w^2/2k_BT`$, which yields precisely the first-order precdiction of Wood *et al.* .

Figure 1 illustrates the behavior of the finite-data free-energy for a Gaussian distribution of work values, based on numerical block averages (2) and the asymptotic behavior given in (11). Although the leading term in $`f_nf`$ is linear in $`1/n`$, the leading coefficient is exponential in the *square* of the distribution’s width, while the next coefficient depends on the *cube* of the width. The asymptotic expressions (11) thus represent viable approximations only for a very small window about $`1/n=0`$ for large widths. Fig. 1 shows that such behavior is easily mistaken for non-analytic (e.g., power-law) behavior.

An exactly solvable case occurs when the Boltzmann factor $`e^wz`$ is distributed according to a gamma distribution, namely,

$$\rho _\mathrm{\Gamma }(z;b,q)=(z/b)^{q1}\mathrm{exp}(z/b)/b\mathrm{\Gamma }(q).$$
(12)
Because this density is “infinitely divisible” (see, e.g., ) the required sums in (3) also follow gamma distributions, and after performing the integration described in (8), one finds

$$f_n(n;b,q)=\mathrm{log}(n/b)\psi (nq)$$
(13)
where the digamma function is defined by $`\psi (x)=(\mathrm{d}/\mathrm{d}x)\mathrm{\Gamma }(x)`$. The exact solution is illustrated in Fig. 2 for $`b=10`$, $`q=2`$.

When asymmetry is added to a Gaussian distribution via the first Edgeworth correction (see, e.g., ), one finds that the exponential dependence of the $`\phi _i`$ on $`\sigma _w`$ is only corrected linearly by the now non-zero third moment of the $`W`$ distribution.

Divergent Moments Case. When the variable $`e^wz`$ in (7) possesses a long-tailed distribution $`\rho _z`$, the limiting distribution is not a Gaussian and the results (9) and (11) no longer hold. In particular, if one of the tails of $`\rho _z(z)`$ decays as $`z^{(1+\alpha )}`$ with $`\alpha <2`$ (implying an infinite Boltzmann variance, $`\widehat{\sigma }^2`$), then the distribution of the variable $`y`$ in (7) approaches a non-Gaussian stable law for large $`n`$ . Note that such power-law behavior in $`z`$ corresponds to *simple exponential decay in the work distribution.* Further, because the mean of $`e^w`$ must be finite for $`\mathrm{\Delta }F`$ to exist \[recall (5)\], we also have $`\alpha >1`$. Unfortunately, no explicit forms for stable distributions are known in the range $`1<\alpha <2`$ .

A long-tailed $`z`$ distribution $`\rho _z\rho _1`$ also alters the *form* of the asymptotic expansion of the sum-variable $`y`$ distribution and, hence, the expansion of $`f_n`$. Instead of (9), the distribution of the variable $`y`$ (7) now takes the more complicated form

$$\rho _n(y)=\rho _\alpha (y)\left[1+\stackrel{}{}\nu _{uv}(y)/n^{\theta (u,v)}\right],$$
(14)
where $`\rho _\alpha `$ is the appropriate stable probability density with exponent $`\alpha `$. The functions $`\nu _{uv}`$, which are not available analytically, depend on the original distribution of $`e^w`$ and partial derivatives of the stable distribution. The exponents are given by $`\theta (u,v)=(u+\alpha v)/\alpha `$, and the summation $`^{}`$ includes $`u0`$ and $`vu/2`$, where $`x`$ denotes the integer part of $`x`$. Note that we have omitted an asymmetry parameter, $`\beta =1`$, of the stable laws  which will be discussed in future work; it does not, however, affect the form of the expansions.

Development of the expansion of $`f_n`$ for large $`n`$ in the case of diverging Boltzmann moments is more complicated, and will only be sketched here. The basic strategy is to ensure that the coefficients of the powers of $`1/n`$ are all rendered in terms of convergent integrals, which requires both an expansion of the logarithm in (8), as well as series *and* asymptotic expansions of $`\rho _\alpha `$ in (14) available from . The asymptotic result for $`n\mathrm{}`$ takes a reasonably simple form, namely,

$$f_nf\phi _{\alpha 1}(1/n)^{(\alpha 1)},$$
(15)
where $`\phi _{\alpha 1}(\alpha )`$ depends on $`\alpha `$ and on the distribution $`\rho _1`$ in a complicated way; details will be presented in a future publication.

Fluctuations and a Universal Law. The fluctuations in the finite-data free energy, $`f_n=\mathrm{\Delta }F_n/k_BT`$, as measured by the variance $`\sigma _n`$ of $`_n`$ of (3), are of considerable interest because of their potential to provide parameter-free extrapolative estimates of $`f_{\mathrm{}}=\mathrm{\Delta }F/k_BT`$ . Formally, the variance is given by

$$\left(\frac{\sigma _n}{k_BT}\right)^2=_{cn^a}^{\mathrm{}}dy\rho _n(y)\left[\mathrm{log}(1+y/cn^a)\right]^2(f_nf)^2.$$
(16)
Using techniques analogous to those sketched above yields asymptotic expansions for the fluctuations. In the case of finite Boltzmann moments, one finds

$$\left(\sigma _n/k_BT\right)^2(\widehat{\sigma }/\widehat{\mu })^2/n+O\left(n^2\right),$$
(17)
where it should be recalled that the unsubscripted moments refer to the density $`\rho _z`$.

Remarkably, comparison with $`\phi _1`$ for (11) shows that

$$f_nf=(\sigma _n/k_BT)^2/2+O(n^2)$$
(18)
exactly, as $`n\mathrm{}`$, and independent of any parameters of the distribution. This universal law, valid for the case when the second Boltzmann moment is finite, is illustrated in Fig. 2. The gamma distribution of Boltzmann factors was $`\rho _\mathrm{\Gamma }(z,10,2)`$; see (12). The “regulated power law” distribution is defined by $`\rho _{rp}(z)=\alpha /(1+z)^{\alpha +1}`$, and we set $`\alpha =2.5`$.

Conclusions. Motivated by the need to understand the large-$`N`$ asymptotic behavior of free-energy-difference estimates based on a finite amount of data ($`N`$ work values), we have presented a general statistical theory which partially completes the task. Two cases were formally identified, distinguished by whether the second moment of the distribution of *Boltzmann factors* of the required work values is finite. The asymptotic behavior was discussed for both cases, and — for the finite-second-Boltzmann-moment case — an exact solution and a universal law were presented.

Much remains to be done, both in terms of theory and applications. A question of particular practical interest is whether parameter-free extrapolation procedures can be devised, particularly in light of the sensitivity of the asymptotic behavior of $`\mathrm{\Delta }F_n`$ to the width of the distribution of work values.

###### Acknowledgements.

The authors have benefitted greatly from discussions with Chris Jarzynski, Hirsh Nanda, Lawrence Pratt, and David Zuckerman. We gratefully acknowledge funding provided by the NIH (under grant GM54782), the Bard Foundation, and the Department of Physiology. D.M.Z. is the recipient of a National Research Service Award (GM20394) from the NIH.