text:
# The depinning transition of a driven interface in the random-field Ising model around the upper critical dimension

## I Introduction

Driven Interfaces in quenched disordered systems display with increasing driving force a transition from a pinned interface to a moving interface (see e.g. KARDAR\_2  and references therein). This so-called depinning transition is caused by a competition of the driving force and the quenched disorder. The first one tends to move the interface whereas the latter one hinders the movement. Depinning transitions are observed in a large variety of physical problems, such as fluid invasion in porous materials (see, for instance, Sec. 6.2 in HINRICHSEN\_1  and references therein), depinning of charge density waves FISHER\_1 ; FISHER\_2 , impurity pinning of flux-line in type-II superconductors BLATTER\_1 , contact lines DEGENNES\_1  as well as in field driven ferromagnets, where the interface separates regions of opposite magnetizations BRUINSMA .

A well established model to investigate the depinning transition in disordered ferromagnets is the driven random-field Ising model (RFIM) (see for instance JI\_2 ; DROSSEL\_1 ; BRUINSMA ; AMARAL\_1 ; AMARAL\_2 ). Here, the disorder induces some effective energy barriers which suppress the interface motion. A magnetic driving field $`H`$ reduces these energy barriers but they vanish only if the driving field exceeds the critical value $`H_\text{c}`$. The transition from the pinned to the moving interface can be described as a continuous phase transition and its velocity $`v`$ is interpreted as the order parameter. Without thermal fluctuations ($`T=0`$) the field dependence of the velocity obeys the power-law behavior

$$v(h,T=0)h^\beta $$
(1)
for $`h>0`$, where $`h`$ denotes the reduced driving field $`h=H/H_\text{c}1`$.

The depinning transition is destroyed in the presence of thermal fluctuations ($`T>0`$) which may provide the energy needed to overcome local energy barriers. Although thermal fluctuations drive the system away from criticality the order parameter obeys certain scaling laws and for sufficiently low temperatures the order parameter can be described as a generalized homogenous function ROTERS\_1

$$v(h,T)=\lambda \stackrel{~}{v}(\lambda ^{1/\beta }h,\lambda ^\psi T),$$
(2)
similar to usual equilibrium second order phase transitions. Setting $`\lambda ^{1/\beta }=1`$ one recovers Eq. (1) for zero temperature. Choosing $`\lambda ^\psi T=1`$ one gets for the interface velocity at the critical field $`H_\text{c}`$

$$v(h=0,T)T^{1/\psi }.$$
(3)
This power-law behavior was observed in two- and three-dimensional simulations of the driven RFIM NOWAK\_1 ; ROTERS\_1  as well as in charge density waves in computer simulations and mean-field calculations FISHER\_2 ; MIDDLETON\_1 .

Furthermore, thermal fluctuations cause a creep motion of the interface for small driving fields ($`HH_\text{c}`$) characterized by an Arrhenius like behavior of the velocity IOFFE\_1 ; NATTERMANN\_4 . Recently, this creep motion was observed in experiments considering magnetic domain wall motion in thin films composed of Co and Pt layers LEMERLE\_1 , in renormalization group calculations CHAUVE\_1 ; CHAUVE\_2  regarding the Edwards-Wilkinson equation with quenched disorder (QEW), as well as in numerical simulations of the RFIM ROTERS\_4 .

In equilibrium physics a scaling ansatz according to Eq. (2) usually describes the order parameter as a function of its control parameter and of its conjugated field. Although Eq. (2) can be applied to the depinning transition, i.e. the temperature is a relevant scaling field, $`T`$ is not conjugated to the order parameter. The conjugated field would support the interface motion independent of its strength. But strong thermal fluctuations destroy the interface instead to support the interface motion. Therefore, one has to interpret the value of the thermal exponent $`\psi `$ carefully. For instance, it is not clear wether the obtained values of $`\psi `$ are a characteristic feature of the whole universality class of the depinning transition or just a characteristic feature of the particular considered RFIM. This point could be important for the interpretation of experiments which naturally take place at finite temperatures.

In this paper we reinvestigate the interface dynamics of the driven RFIM and focus our attention to higher dimensions $`d3`$. In particular we consider the scaling behavior at the critical point and determine the exponents $`\beta `$ and $`\psi `$. Our results suggest that the so-called upper critical dimension of the depinning transition of the driven RFIM is $`d_\text{c}=5`$. Above this dimension the scaling behavior is characterized by the mean-field exponents. We compare our results with those obtained from a renormalization group approach of the quenched Edward-Wilkinson equation which is expected to be in the same universality class as the driven RFIM. A summary is given at the end.

## II Model and simulations

We consider the depinning transition RFIM on cubic lattices of linear size $`L`$ in higher dimensions ($`d3`$). The Hamiltonian of the RFIM is given by

$$=\frac{J}{2}\underset{i,j}{}S_iS_jH\underset{i}{}S_i\underset{i}{}h_iS_i\text{ ,}$$
(4)
where the first term characterizes the exchange interaction of neighboring spins ($`S_i=\pm 1`$). The sum is taken over all pairs of neighbored spins. The spins are coupled to a homogenous driving field $`H`$ as well as to a quenched random-field $`h_i`$ with $`h_i=0`$ and $`h_ih_j\delta _{ij}`$. The random field is assumed to be uniformly distributed, i.e., the probability $`p`$ that the random field at site $`i`$ takes some value $`h_i`$ is given by

$$p(h_i)=\{\begin{array}{ccc}(2\mathrm{\Delta })^1& \mathrm{for}& |h_i|<\mathrm{\Delta }\hfill \\ 0& & \text{otherwise.}\hfill \end{array}$$
(5)
Using antiperiodic boundary conditions an interface is induced into the system which can be driven by the field $`H`$ (see ROTERS\_1  for details). A Glauber dynamics with random sequential update and heat-bath transition probabilities is applied to simulate the interface motion (see for instance BINDER\_1 ).

Due to this algorithm not only spins adjacent to the interface but throughout the whole system can flip with a finite probability at temperatures $`T>0`$. In general, this may cause nucleation which always starts with an isolated spin-flip. The minimum energy required for an isolated spin-flip is $`\mathrm{\Delta }E=2(zJH\mathrm{\Delta })`$ with $`z=2^d`$ nearest neighbors on a bcc lattice. Although the corresponding spin-flip probability $`1/[1+\mathrm{exp}(\mathrm{\Delta }E/T)]`$ is small for the considered temperatures ($`T0.2`$), isolated spin-flips are possible and occur. But we observed that these spin-flips are unstable in our simulations, i.e., an isolated spin will flip back in the next update step. Thus the originally induced interfaces is stable during the simulations.

The moving interface corresponds to a magnetization $`M`$ which increases in time $`t`$ (given in Monte Carlo step per spin). The interface velocity, which is the basic quantity in our investigations, is obtained from the time dependence of the magnetization $`v=\text{d}M/\text{d}t`$ where $`\mathrm{}`$ denotes an appropriate disorder average. Starting with a flat interface we performed a sufficiently number of updates until the system reaches after a transient regime the steady state which is characterized by a constant average interface velocity.

As pointed out in previous works NOWAK\_1 ; ROTERS\_1  an appropriate choice of the interface orientation is needed in order to recover that the interface moves for arbitrarily small driving field in the absence of disorder. An appropriate choice is to consider the interface motion along the diagonal direction of a simple cubic lattice. For $`d3`$ an alternative is to examine the interface motion along the $`z`$-axis on a body-centered cubic (bcc) lattice. Since it is much more convenient to implement the latter case in higher dimensions we consider in this work bcc lattices, the more as the lattice structure usually does not affect the universal scaling behavior.

## III d=3

In the case of the three dimensional RFIM we consider bcc lattices of linear size $`L250`$. A snapshot of a moving interface in the steady state is presented in Fig. 1. The obtained values of the interface velocities for $`T=0`$ are plotted in Fig. 2. As one can see $`v`$ tends to zero in the vicinity of $`H1.36`$. Assuming that the scaling behavior of the interface motion is given by Eq. (1) one varies $`H_\text{c}`$ until one gets a straight line in a log-log plot. Convincing results are obtained for $`H_\text{c}=1.357\pm 0.001`$ and the corresponding curve is shown in the inset of Fig. 2. For lower and greater values of $`H_\text{c}`$ we observe significant curvatures in the log-log plot (not shown). In this way we estimate the error-bars in the determination of the critical field. A regression analysis yields the value of the order parameter exponent $`\beta =0.653\pm 0.026`$. This value agrees with $`\beta =0.66\pm 0.04`$ which was obtained from a similar investigation ROTERS\_1  where the interface moves along the diagonal direction of a simple cubic lattice (see inset of Fig. 2). Furthermore, both results are in agreement with $`\beta =0.60\pm 0.11`$ AMARAL\_1 ; AMARAL\_2 , where in $`d=3`$ the influence of helicoidal boundary conditions in one direction and periodic ones in the other direction parallel to the interface was investigated on a simple cubic lattice.

In order to determine the exponent $`\psi `$ we simulated the RFIM around the critical field for $`T=0.025n`$ with $`n\{1,2,3,4,6,8\}`$. The obtained curves are shown in the inset of Fig. 3. According to Eq. (2) the interface velocity scales as

$$v(h,T)=T^{1/\psi }\stackrel{~}{v}(hT^{1/\beta \psi },1).$$
(6)
Plotting $`v(h,T)T^{1/\psi }`$ as a function of $`hT^{1/\beta \psi }`$ one varies $`\beta `$, $`\psi `$, and $`H_\text{c}`$ until one gets a data collapse of the different curves. Convincing data collapses are observed for $`\beta =0.63\pm 0.06`$ $`\psi =2.33\pm 0.2`$ and $`H_\text{c}=1.360\pm 0.01`$ and the corresponding curves are shown in Fig. 3. The obtained values of the order parameter exponent and of the critical field agree within the error-bars with the values of the $`T=0`$ analysis. Furthermore our results are in agreement with similar investigations on a simple cubic lattice ($`\beta =0.63\pm 0.07`$ and $`\psi =2.38\pm 0.2`$, see ROTERS\_1 ).

## IV d=4

In order to determine the order parameter exponent of the four dimensional driven RFIM we measured the interface velocity for bcc lattices of linear sizes $`L140`$. The obtained data for $`T=0`$ are shown in a log-log plot in Fig. 4. After a transient regime which displays a finite curvature we observe an asymptotic power-law behavior for sufficiently small driving field $`H`$. A regression analysis yields $`\beta =0.8\pm 0.06`$ and $`H_\text{c}=1.258\pm 0.002`$.

To determine the exponent $`\psi `$ we simulated the RFIM in the vicinity of the critical field for $`T=0.025n`$ where again $`n\{1,2,3,4,6,8\}`$ was choosen. Similar to the three dimensional case one varies the exponents as well as the critical field until one observes a data collapse. Good results are obtained for $`\beta =0.73\pm 0.13`$, $`\psi =1.72\pm 0.27`$, and $`H_\text{c}=1.256\pm 0.015`$. The corresponding data collapse is shown in Fig. 5. Again, the obtained values of $`\beta `$ and $`H_\text{c}`$ confirm the above presented analysis for $`T=0`$.

## V d=5

In the case of the five-dimensional RFIM system sizes from $`L=10`$ up to $`L=30`$ are simulated. Analyzing the interface motion at $`T=0`$ we observe that the velocity-field dependence can not be described by a pure power-law, i.e., Eq. (1) fails. In Fig. 6 we plotted the logarithmic derivation of the velocity-field dependence

$$\beta _{\text{eff}}=\frac{\mathrm{ln}v}{\mathrm{ln}h}$$
(7)
which can be interpreted as an effective exponent. If the asymptotic scaling behavior obeys Eq. (1) the logarithmic derivative tends to the value of $`\beta `$ for $`HH_\text{c}`$. But as can be seen from Fig.6 no clear saturation takes place for $`d=5`$ as it is observed for the three and four dimensional model. The lack of a clear saturation could be explained by a too large large distance from the critical point, but another reason is possible too.

Significant deviations from a pure power law behavior \[Eq. (1)\] occur for instance at the upper critical dimensional $`d_\text{c}`$ where the scaling behavior is governed by the mean-field exponents modified by logarithmic corrections. The scaling behavior around $`d_\text{c}`$ is well understood within the renormalization group theory (see for instance WEGNER\_2 ; BREZIN\_3 ; PFEUTY\_1 ). For $`d>d_\text{c}`$ the stable fix point of the corresponding renormalization equations is usually a trivial fix point with classical mean-field exponents. This trivial fix point is unstable for $`d<d_\text{c}`$ and a different stable fix point exists with nonclassical exponents. These exponents can be estimated by an $`ϵ`$-expansion, for instance. For $`d=d_\text{c}`$ both fix points are identical and marginally stable. In this case the asymptotic form of the thermodynamic functions is given by the mean-field power-law behavior modified by logarithmic corrections. Applying this approach to the depinning transition, the corresponding ansatz reads

$$v(h,T=0)h^{\beta _{\text{MF}}}|\mathrm{ln}h|^\text{B},$$
(8)
where B denotes an unknown correction exponent. It is worth to mention that in contrast to the values of the critical exponents below the upper critical dimension the above scaling behavior does not rely on approximation schemes like $`ϵ`$\- or $`1/n`$-expansions BREZIN\_2 . Within the renormalization group theory it is an exact result in the limit $`h0`$ (see e.g. WEGNER\_2 ; WEGNER\_1  and references therein for RG investigations and GRIFFIN\_1 ; BRINKMANN\_1  for measurements).

The value $`\beta _{\text{MF}}=1`$ is reported for depinning transitions KARDAR\_2 ; NATTERMANN\_1 ; CHAUVE\_3 . Thus we analyze $`v(h)/h`$ as a function of $`|\mathrm{ln}h|`$ and note again that Eq. (8) describes only the leading order of the scaling behavior, i.e., we expect that the asymptotic behavior of interface velocity obeys

$$[v(h)/h]^{1/\text{B}}=\mathrm{const}|\mathrm{ln}h|.$$
(9)
Therefore, we varied in our analysis the logarithmic correction exponent B and the critical field $`H_\text{c}`$ until we get this expected asymptotic behavior. The best results are obtained for $`\text{B}=0.40\pm 0.09`$ and $`H_\text{c}=1.14235\pm 0.001`$ and the corresponding scaling plot is shown in Fig. 7. The observed asymptotic agreement with Eq. (9) corresponds to a logarithmically $`(1/|\mathrm{ln}h|)`$ convergence of $`\beta _{\text{eff}}`$ to $`\beta _{\text{MF}}=1`$, which explains why no clear saturation of the effective exponent could be observed for $`h0`$ in the five dimensional model.

Similar to the $`T=0`$ scaling behavior one has to modify for $`T>0`$ the scaling ansatz since no data collapse could be obtained by plotting the data according to Eq. (2). Motivated by recently performed investigations of the scaling behavior of an absorbing phase transition around the upper critical dimension LUEB\_22  we assume that the scaling behavior of the order parameter obeys in leading order

$$v(h,T)=T^{1/\psi _{\text{MF}}}|\mathrm{ln}T|^\sigma \stackrel{~}{v}(x,1)$$
(10)
where the scaling argument $`x`$ is given in leading order by

$$x=hT^{1/\beta _{\text{MF}}\psi _{\text{MF}}}|\mathrm{ln}T|^\tau $$
(11)
with $`\beta _{\text{MF}}=1`$. In our analysis we use the value $`\psi _{\mathrm{MF}}=1.49`$ obtained from the analysis of the six dimensional RFIM (see next section). Therefore, we have to vary the exponents $`\sigma `$ and $`\tau `$ in order to observe a data collapse according to Eqs. (10,11). Convincing results are obtained for $`\sigma =0.22\pm 0.16`$, $`\tau =0.19\pm 0.12`$ The corresponding data collapse is shown in Fig. 8.

## VI d=6

Above the upper critical dimension the scaling behavior is characterized by the mean-field exponents, i.e., in leading order the interface velocity is given by

$$v(h,T=0)h.$$
(12)
In Fig. 9 we plot the velocity as a function of the driving field $`H`$ obtained from simulations of system sizes $`L14`$. As one can see the velocity does not display the expected linear behavior. It seems that a linear behavior is only given for small velocities, i.e., our data do not display the pure asymptotic behavior \[Eq. (12)\]. This is confirmed by the behavior of the effective exponent \[Eq.(7)\] which increases fast for $`h0`$ but the actual saturation to $`\beta =1`$ does not take place for the considered values of $`h`$ (see Fig. 6). To observe the asymptotic behavior one has to perform simulations closer to the critical point $`H_\text{c}`$ which requires to simulate larger system sizes. Unfortunately, the limited CPU power makes this impossible.

An alternative is to take the curvature of the function $`v(h)`$ into consideration and to assume that the leading corrections to the asymptotic behavior are of the form

$$v(h,T=0)=v_1h+v_2h^2+𝒪(h^3),$$
(13)
which recovers Eq. (12) for $`h0`$. Fitting our data to this ansatz we get $`H_\text{c}=1.1537\pm 0.003`$, $`v_1=0.2546`$, and $`v_2=0.0895`$. The corresponding curve fits the simulation data quite well as one can see from Fig. 9. In the inset of Fig. 9 we plotted $`v(h)/h`$ as a function of the reduced driving field $`h`$. According to the above ansatz \[Eq. (13)\] one gets a linear behavior, i.e., the deviations from the pure mean-field behavior \[Eq. (12)\] can really be described as quadratic corrections. Thus we get that our numerical data are consistent with the assumption that the six dimensional RFIM depinning transition is characterized by the mean-field exponent $`\beta _{\text{MF}}=1`$.

Again we consider how thermal fluctuations affect the scaling behavior and analyse interface velocities obtained at different temperatures $`T=0.025n`$ with $`n\{1,2,3,4,6,8\}`$. Similar to the situation below the upper critical dimension we assume that the scaling behavior of the interface velocity is given by Eq. (2) where the exponents are given by mean-field values. A convincing data collapse is obtained for $`\psi _{\text{MF}}=1.49\pm 0.15`$ and $`H_\text{c}=1.153\pm 0.02`$ and is plotted in the inset of Fig. 5.

## VII Discussion

A well established realization of interface pinning in a disordered media is the so-called quenched Edwards-Wilkinson (QEW) equation of motion which was intensively investigated in the last decade BRUINSMA ; NARAYAN\_1 ; NATTERMANN\_1 ; CHAUVE\_1 ; CHAUVE\_2 ; CHAUVE\_3 . It is argued that the QEW equation as well as the driven RFIM are characterized by the same critical exponents, i.e., both models belong to the same universality class BRUINSMA ; AMARAL\_1 . Renormalization group analyses of the quenched QEW equation NATTERMANN\_1 ; CHAUVE\_3  predict, in accordance with BRUINSMA , $`d_\text{c}=5`$ and allow to estimate the critical exponents using an $`ϵ`$-expansion. A recently performed two-loop renormalization approach yields CHAUVE\_3

$$\beta _{\text{QEW}}=\mathrm{\hspace{0.33em}1}\frac{1}{9}ϵ0.040123ϵ^2+𝒪(ϵ^3)$$
(14)
where $`ϵ`$ denotes the distance from the upper critical dimension, i.e. $`ϵ=5d`$ (unfortunately, no error-bars can be estimated from an $`ϵ`$-expansion). The corresponding values of the exponents as a function of the dimension are plotted in Fig. 10. The numerically determined exponents $`\beta `$ of the driven RFIM (listed in Table 1) are in a fair agreement with the values of the $`ϵ`$-expansion.

For the QEW equation the temperature exponent $`\psi `$ is not known. Therefore, a direct comparison with the obtained values of the driven RFIM is not possible.

## VIII Conclusions

We studied numerically a field driven interface in the RFIM and determined the order parameter exponent $`\beta `$ as well as the temperature exponent $`\psi `$. Below the upper critical dimension $`d_\text{c}=5`$ the critical exponents depend on the dimension and the values of the exponents correspond to those of a two-loop renormalization group approach of the Edwards-Wilkinson equation CHAUVE\_3 . This suggest that the depinning transition of the RFIM model belongs to the universality class of the quenched Edward-Wilkinson equation. At the upper critical dimension $`d_\text{c}=5`$ the scaling behavior is affected by logarithmical corrections. Above the upper critical dimension we observe that the scaling behavior is characterized in leading order by the corresponding mean-field exponents.

###### Acknowledgements.

This work was supported by the Deutsche Forschungsgemeinschaft via GRK 277 Struktur und Dynamik heterogener Systeme (University of Duisburg) and SFB 491 Magnetische Heteroschichten: Struktur und elektronischer Transport (Duisburg/Bochum). S. Lübeck wishes to thank the Minerva Foundation (Max Planck Gesellschaft) for financial support.


format: text

text:
# Why temperature chaos in spin glasses is hard to observe

## Abstract

The overlap length of a three-dimensional Ising spin glass on a cubic lattice with Gaussian interactions has been estimated numerically by transfer matrix methods and within a Migdal-Kadanoff renormalization group scheme. We find that the overlap length is large, explaining why it has been difficult to observe spin glass chaos in numerical simulations and experiment.

Chaos, rejuvenation, memory, and aging in spin glasses are currently being intensively studied Jonason et al. (1998); Bouchaud et al. (2001); Rizzo (2001); Yoshino et al. (2001); Berthier and Bouchaud (2002); Krzakala and Martin (2002); Jönsson et al. (2002); Takayama and Hukushima (2002); Billoire and Marinari (2000). “Chaos” refers to the property that equilibrium states in the ordered phase of spin glasses are sensitive to arbitrarily small changes in the couplings or in temperature, and is one possible explanation for memory effects in spin glasses Jonason et al. (1998), although other mechanisms are also present Berthier and Bouchaud (2002). There is evidence both for and against temperature chaos in computer simulations: some authors claim to see the effects of chaos Takayama and Hukushima (2002), others have failed to see it Billoire and Marinari (2000, 2002). In this Letter we start from the droplet Fisher and Huse (1986) or scaling Bray and Moore (1987) pictures and calculate numerically the overlap length $`L^{}(T,\mathrm{\Delta }T)`$, i.e. the length scale beyond which spins at temperatures $`T`$ and $`T+\mathrm{\Delta }T`$ become uncorrelated with each other, for a three-dimensional Edwards-Anderson Ising spin glass and also within a Migdal-Kadanoff renormalization group scheme (MKRG), both with a Gaussian distribution of couplings with unit variance. We show that this length scale is large within much of the parameter space, and only just comes down to magnitudes accessible to numerical simulations for certain values of $`T`$ and $`\mathrm{\Delta }T`$. We believe that this is why some workers have been unable to see chaos in their investigations. As a by-product, we will also find the root mean square droplet interface free energy $`F(T)\sqrt{F_{\text{int}}^2(T)}`$ and entropy $`S(T)\sqrt{S_{\text{int}}^2(T)}`$. The angle brackets indicate averaging over realizations of the bond couplings. We show that $`S(T)\sqrt{T}`$ for small $`T`$, contrary to previous arguments Huse and Fisher (1986).

In the MKRG scheme for three dimensions, renormalized bonds after $`n`$ renormalization steps $`J^{(n)}`$ are obtained from the set of bonds $`\{J^{(n1)}\}`$ after $`n1`$ renormalization steps using the relation

$$\frac{J^{(n)}}{T}=\underset{i=1}{\overset{4}{}}\mathrm{tanh}^1\left(\mathrm{tanh}\frac{J_{1i}^{(n1)}}{T}\mathrm{tanh}\frac{J_{2i}^{(n1)}}{T}\right),$$
(1)
where $`J_{1,2i}^{(n1)}`$ are randomly drawn members of the bond pool (“pool method”). For a detailed description of the method in the present context we refer the reader to Banavar and Bray (1987). The renormalized bonds play the role of the interface free energy of a system of linear size $`L=2^n`$, to wit $`\frac{1}{2}F(T)=(J^{(n)})^2^{1/2}`$, the angle brackets indicating the average over the bond pool. Similarly, the interface entropy is given by $`\frac{1}{2}S(T)=lim_{\delta T0}(J^{(n)}J_{}^{}{}_{}{}^{(n)})^2^{1/2}/\delta T`$, where $`J_{}^{}{}_{}{}^{(n)}`$ is the corresponding member of a bond pool which has been evolved at temperature $`T+\delta T`$.

The MKRG applied to chaos in spin glasses has the advantages that large length scales are easily accessible numerically and that it is possible to access and estimate the overlap length directly. Banavar and Bray Banavar and Bray (1987) have shown that the overlap length $`L^{}`$ is related to the ratio of the interface entropy and interface free energy. Introducing the usual generalized stiffness coefficients $`\mathrm{{\rm Y}}(T)`$ and $`\sigma (T)`$ via

$$F(T)=\mathrm{{\rm Y}}(T)L^\theta \text{ and }S(T)=\sigma (T)L^{d_s/2},$$
(2)
(valid for large system sizes $`L`$), a measure of the overlap length for a temperature step $`\mathrm{\Delta }TT`$ is obtained by equating $`F(T)`$ and $`S(T)\mathrm{\Delta }T`$:

$$L^{}=\left(\frac{\mathrm{{\rm Y}}}{\sigma \mathrm{\Delta }T}\right)^{1/\zeta }.$$
(3)
Here, $`\zeta =\frac{d_s}{2}\theta `$ is the chaos exponent, where $`d_s`$ is the fractal dimension of the droplet interface and $`\theta `$ is the droplet excitation energy exponent.

Independently, the overlap length can be obtained as the length scale on which two initially identical bond pools $`\{J\}`$ and $`\{J^{}\}`$, evolved at two different temperatures, decorrelate. A measure for this is $`L^{}=2^{n_0}`$ where $`n_0`$ is such that

$$\frac{\left(J^{(n_0)}J_{}^{}{}_{}{}^{(n_0)}\right)^2}{(J^{(n_0)})^2+(J_{}^{}{}_{}{}^{(n_0)})^2}=\frac{1}{2}.$$
(4)
This method does not require a small temperature difference. When the temperature difference does happen to be small, however, it agrees with the previous method up to a numerical factor close to unity (see below).

Repeating the numerical work from Banavar and Bray (1987) with greater accuracy (pool size 2000000 and ‘long double’ precision) and for a larger set of temperatures, we find $`\mathrm{{\rm Y}}(T)`$ and $`\sigma (T)`$ as shown in Fig. 1.

The data for $`\sigma (T)`$ in this figure has been obtained using a fixed relative temperature change $`ϵ=\delta T/T=10^6`$ at each temperature.

Huse and Fisher Huse and Fisher (1986) argued that a droplet interface can be regarded as a collection of roughly independent two level systems, each of which has on average an entropy proportional to $`T`$, which should give rise to $`S(T)T`$. This is clearly violated by the data in Fig. 1. We will come back to this point later.

The overlap length $`L^{}`$ as obtained from the data for $`\mathrm{{\rm Y}}`$ and $`\sigma `$ using Eq. (3) is shown in Fig. 2.

This requires knowing $`\zeta `$. For the $`d`$-dimensional Migdal-Kadanoff scheme, $`d_s=d1`$ (i.e. $`d_s=2`$ in the present case), and $`\theta `$ was estimated from the numerics to be $`\theta =0.25519\pm 0.00005`$ for $`d=3`$, therefore $`\zeta 0.7448`$. For illustration, the temperature shift has been arbitrarily set to $`\mathrm{\Delta }T=0.01`$ and $`\mathrm{\Delta }T=0.1`$ for this plot. Since for these values of the parameters the temperature shift $`\mathrm{\Delta }T`$ is not much smaller than $`T`$ for small $`T`$, Eq. 3 is not expected to hold and the “true” overlap length differs from the data shown. A direct determination of the overlap length (as the decorrelation length of the two bond pools as explained above) with $`\mathrm{\Delta }T=0.01`$ is shown in Fig. 2 for comparison. A small constant has been added to the latter curve to make the two ways of determining $`L^{}`$ agree at the minimum. (This adjustment reflects the arbitrariness of the constant 1/2 in the definition, through Eq. (4), of $`L^{}`$.) With this correction, the two curves for $`\mathrm{\Delta }T=0.01`$ agree very well in a region around the minimum; they differ for small $`T`$ since $`\mathrm{\Delta }TT`$ and for $`T`$ close to the critical temperature because of increasing influence of critical point controlled decorrelation Nifle and Hilhorst (1992).

Fig. 3 shows a contour plot of the overlap length in the two-temperature plane, obtained from MKRG in the same manner as before.

This figure shows that the overlap length is nowhere less than 40 for $`T_1,T_2<T_c=0.89645\mathrm{}`$. Note that in agreement with Nifle and Hilhorst (1992) the overlap length is nonzero even above $`T_c`$.

Obviously, $`L^{}`$ is extraordinarily large for most values of $`T`$ and never comes into a numerically accessible range, which, by today’s standards, would be around $`L^{}20`$.

While the Migdal-Kadanoff renormalization method allows study of large length scales, it has drawbacks: the fractal dimension of a domain wall $`d_s`$ on a cubic lattice is not well approximated, so the structure of an interface is significantly different on a hierarchical Berker lattice Berker and Ostlund (1979) (for which MKRG is exact) than on a cubic lattice. It is therefore not obvious that the overlap length as obtained above is characteristic of a three-dimensional Edwards-Anderson spin glass as well. In order to test this, we used the transfer matrix method as described in Bray and Moore (1985) to numerically obtain the interface free energy and entropy of small spin glass samples. This method has the disadvantage that due to demands on computer time and memory, only samples of size up to $`L=4`$ are accessible with reasonable statistics. On the positive side, there are no approximations involved whatsoever. This approach is therefore complimentary to the Migdal-Kadanoff scheme. We shall use Eq. (3), which applies generally wherever the droplet picture holds, to provide estimates of $`L^{}`$ from measurements of $`\mathrm{{\rm Y}}`$ and $`\sigma `$.

In order to obtain the interface free energy $`F_{\text{int}}`$ of a spin glass sample of size $`L\times L\times (L+1)`$ at temperature $`T`$, first the free energy of the sample with periodic boundary conditions along the first two dimensions and fixed boundary conditions ($`+1`$ on both sides) along the third dimension is calculated. Next, the free energy of the same sample with boundary condition $`+1`$ on the one side and $`1`$ on the other is calculated. The difference is the interface free energy. The interface entropy $`S_{\text{int}}`$ is obtained by repeating the above for temperature $`T+\delta T`$ and taking the numerical derivative of $`F_{\text{int}}`$. From these data, the root mean square averages $`F(T)`$ and $`S(T)`$ are calculated.

The data for the interface free energy and entropy are shown in Fig. 4 as scaling plots for system sizes $`L=2,3,4`$. The entropy has been obtained with a fixed relative temperature shift $`ϵ=\delta T/T=10^4`$.

The plot for $`F(T)`$ allows for an estimate of $`\theta `$, the one for $`S(T)`$ for an estimate of $`d_s`$ (cf. Eq. (2)). The values obtained are $`\theta 0.18`$ and $`d_s2.38`$, leading to $`\zeta 1.01`$. While these numbers are merely crude estimates (corrections to scaling are expected to be significant for these small system sizes) and therefore no error bars have been supplied, they certainly lie in the expected range and the value of $`d_s`$, in particular, shows that the droplet interface structure is much better captured even for these small systems than in the MKRG. The scaling in Fig. 4 naturally breaks down close to the critical temperature, due to the small system sizes.

The fractal dimension $`d_s`$ can also be estimated at zero temperature from a bond perturbation calculation as in Bray and Moore (1987). The result, shown in Fig. 5, is $`d_s2.7`$, obtained from sample sizes $`L=2,3,4,5`$.

The discrepancy between this result and the one above illustrates the influence of finite size effects. A similar comparison in two dimensions 111Since $`T_c=0`$ in two dimensions, the entropy scaling only works if $`L`$ is much smaller than the correlation length, i.e. only at small temperatures., also shown in Fig. 5, gives the same kind of discrepancy for small $`L`$, which however is removed when going to larger $`L`$. It is found that the value obtained from bond perturbation is more reliable than the one obtained from a scaling plot of the entropy; the bond perturbation data lie almost perfectly on a straight line even down to $`L=2`$. In three dimensions, the same seems to be true since $`d_s2.7`$ compares well with $`d_s2.68`$ from Palassini and Young (1999).

The interface entropy in Fig. 4 again scales as $`S(T)\sqrt{T}`$, showing that this behavior is not special to the MKRG. In order to understand this, we return to the two-level-system argument of Huse and Fisher (1986). The interface entropy is calculated from the change in entropy after a change of boundary conditions, therefore it is necessary to know the excitation energy of each two level system under both boundary conditions. That these will in general be different becomes apparent when considering, e.g., a single spin on the interface: since it is on the interface, at least one but not all of its neighboring spins change sign when changing boundary conditions, thus giving rise to a different effective field for the spin, i.e. a different excitation energy for flipping this spin. Therefore, the interface entropy contribution from a two-level system is in fact $`\mathrm{\Delta }S(\mathrm{\Delta }^+,\mathrm{\Delta }^{})=S_{\text{2-lev}}(\mathrm{\Delta }^+)S_{\text{2-lev}}(\mathrm{\Delta }^{})`$, where $`\mathrm{\Delta }^\pm `$ are the excitation energies under the two boundary conditions. It is easy to check that the entropy of a two-level system is

$$S_{\text{2-lev}}(\mathrm{\Delta })=\mathrm{log}\left(2\mathrm{cosh}\frac{\mathrm{\Delta }}{2T}\right)\frac{\mathrm{\Delta }}{2T}\mathrm{tanh}\frac{\mathrm{\Delta }}{2T}.$$
(5)
The function $`S_{\text{2-lev}}(\mathrm{\Delta })`$ has a maximum at $`\mathrm{\Delta }=0`$ and decays to zero on scale $`T`$, therefore $`\mathrm{\Delta }S(\mathrm{\Delta }^+,\mathrm{\Delta }^{})`$ is non-zero essentially only in two perpendicular strips of width $`T`$ along the $`\mathrm{\Delta }^\pm `$ axes, excluding their overlap region around the origin. Thus only those two-level systems contribute to the interface entropy which have excitation energies $`\mathrm{\Delta }^+<T`$ and $`\mathrm{\Delta }^{}>T`$ (or vice versa). This implies that the second moment, taken with the joint probability distribution $`p(\mathrm{\Delta }^+,\mathrm{\Delta }^{})`$, goes as

$$\mathrm{\Delta }S^2_T^{\mathrm{}}𝑑\mathrm{\Delta }^+_0^T𝑑\mathrm{\Delta }^{}p(\mathrm{\Delta }^+,\mathrm{\Delta }^{})T,T0,$$

provided the marginal distribution $`_0^{\mathrm{}}𝑑\mathrm{\Delta }^+p(\mathrm{\Delta }^+,\mathrm{\Delta }^{})`$ is nonzero for $`\mathrm{\Delta }^{}=0`$. This result is very well supported by the data in Figs. 1 and 4.

The overlap length as estimated from $`F(T)`$, $`S(T)`$, and Eq. (3) is shown in Fig. 6.

While being smaller than for the MKRG, $`L^{}`$ is still very large ($`L^{}>30`$ for $`\mathrm{\Delta }T<0.1`$) and thus (just) out of range of numerical simulations. In fact with the parameters as used in Takayama and Hukushima (2002) ($`T_1=0.7`$, $`T_2=0.4`$), the overlap length is estimated to be $`L^{}20`$, while the authors are using a system of size $`L=24`$. This is in accord with the authors’ claim that they have been able to observe the effects of temperature chaos. Again, we would expect modifications to the small temperature behavior of $`L^{}`$ if a fixed absolute temperature shift is used (cf. Fig. 2), but not in the region of the minimum.

Finally it is worth remarking that the large value of the overlap length will have implications also for experiments on memory and rejuvenation. The effects of chaos are only visible when the length scale $`L(t)`$ within which the spins are well-equilibrated after waiting for a time $`t`$, is larger than the overlap length $`L^{}`$. As $`L(t)`$ only increases very slowly with time, many experiments will not satisfy the criterion $`L(t)>L^{}`$ and it will then be inappropriate to use chaos ideas to explain what is happening.

###### Acknowledgements.

T.A. acknowledges support by the German Academic Exchange Service (DAAD) under a postdoc fellowship.


text:
# Remarks on time-energy uncertainty relations

## 1 Introduction

Time-energy uncertainty relations played an important rôle in the early discussions on the physical interpretation of quantum theory . But contrary to the position-momentum uncertainty relation, their derivation and even precise formulation suffer from the difficulty of assigning an observable (in the sense of selfadjoint operators) of quantum theory with the measurement of time . Meanwhile, as advocated long ago by Ludwig , it is widely accepted that the concept of observables should be generalized, by allowing not only selfadjoint operators (corresponding to projection valued measures) but also positive operator valued measures , and it was quickly realized that one can find such measures which transform covariantly under time translations and fulfil therefore all formal requirements for a time observable . In a recent paper  we gave the first (to our knowledge) general construction of such measures starting from an arbitary positive operator which may be interpreted as the effect whose occurence time is described by the measure; this notion of time is closely related to the concept of time of arrival but in contrast to this our construction always leads to positive operator valued measures (cfr. the book in ). All that is done in a strict quantum language, no classical ideas or generalizations of quantum mechanics are involved (we stress that the construction is valid also in quantum field theory, a subject which we plan to deal with in the future). In the present paper we show that for these time observables the usual time-energy uncertainty relation holds, and that in addition, provided the Hamiltonian is positive, one finds an uncertainty relation for time measurements alone which takes the form

$$\mathrm{\Delta }T\frac{\mathrm{const}}{H}$$

with a universal constant. Arguments for these relations have been given by several authors , mainly in special situations and often only on a heuristic level. It is the aim of the present note to show that these relations can be rigorously derived for an arbitrary time translation covariant positive operator valued measure.

## 2 Covariant Naimark-Stinespring Dilation

Let $``$ be a Hilbert space and $`(U(t))`$ a strongly continuous unitary group on $``$ describing the time evolution. The occurence time of an effect is, following , described by a covariant positive operator valued measure $`F`$, i.e. for any Borel subset $`B`$ of the real line $`𝐑`$ we have a positive contraction $`F(B)`$ on $``$, and the operators $`F(B)`$ satisfy the conditions

$$F(\underset{n}{}B_n)=\underset{n}{}F(B_n)$$

if the sets $`B_n`$ are pairwise disjoint and where the sum on the r.h.s. converges strongly,

$$F(𝐑)=1,F(\mathrm{})=0,$$

$$U(t)F(B)U(t)=F(B+t).$$

We want to construct from these data an extended Hilbert space $`𝒦`$ with a projection $`P`$ onto the subspace $``$, a unitary group $`(V(t))`$ on $`𝒦`$ which reduces on $``$ to the original time translation $`(U(t))`$ and a covariant projection valued measure $`E`$ on $`𝒦`$ such that

$$PE(B)P=F(B)P.$$

Let $`𝒦_0`$ be the space of bounded piecewise continuous functions on $`𝐑`$ with values in $``$. On this space we introduce the positive semidefinite scalar product

$$\mathrm{\Phi },\mathrm{\Psi }=(\mathrm{\Phi }(t),F(dt)\mathrm{\Psi }(t)).$$

$``$ can be isometrically embedded into $`𝒦_0`$ by identifying the elements of $``$ with constant functions. The enlarged Hilbert space $`𝒦`$ is then defined as the completion of the quotient of $`𝒦_0`$ by the null space of the scalar product.

On $`𝒦_0`$ we define

$$(E(B)\mathrm{\Phi })(t)=\{\begin{array}{ccc}\mathrm{\Phi }(t)& ,& tB\\ 0& ,& tB\end{array}$$

$$(V(t)\mathrm{\Phi })(s)=U(t)\mathrm{\Phi }(st),$$

$$P\mathrm{\Phi }=F(dt)\mathrm{\Phi }(t).$$

$`E(B)`$ and $`V(t)`$ map the null space into itself, $`P`$ annihilates it. Therefore they are well defined operators on $`𝒦`$ and form the desired covariant dilation. In particular, $`(E,V)`$ is a system of imprimitivity over $`𝐑`$  and therefore unitarily equivalent to a multiple of the Schrödinger representation.

## 3 The general Time-Energy Uncertainty Relation

Using the covariant dilation described in the previous section we can use the standard uncertainty relation in the Schrödinger representation on $`L^2(𝐑)`$,

$$\mathrm{\Delta }_\mathrm{\Phi }(x)\mathrm{\Delta }_\mathrm{\Phi }(\frac{1}{i}\frac{d}{dx})\frac{1}{2}$$

which holds for all wave functions $`\mathrm{\Phi }`$ in the intersection of the domains of $`x`$ and $`\frac{1}{i}\frac{d}{dx}`$. This follows from the validity of the canonical commutation relations in the sense of quadratic forms,

$$(x\mathrm{\Phi },\frac{1}{i}\frac{d}{dx}\mathrm{\Psi })(\frac{1}{i}\frac{d}{dx}\mathrm{\Phi },x\mathrm{\Psi })=i(\mathrm{\Phi },\mathrm{\Psi })$$

which may be derived from the fact that $`\frac{1}{i}\frac{d}{dx}`$ is the generator of translations $`U(a)`$, and that, by Stone’s Theorem , $`aU(a)\mathrm{\Phi }`$ is strongly differentiable for $`\mathrm{\Phi }`$ in the domain of $`\frac{1}{i}\frac{d}{dx}`$. Namely, we have

$$(x\mathrm{\Phi },\frac{1}{i}\frac{d}{dx}\mathrm{\Psi })(\frac{1}{i}\frac{d}{dx}\mathrm{\Phi },x\mathrm{\Psi })$$

$$=\frac{1}{i}\frac{d}{da}|_{a=0}\left((x\mathrm{\Phi },U(a)\mathrm{\Psi })(U(a)\mathrm{\Phi },x\mathrm{\Psi })\right)$$

$$=\frac{1}{i}\frac{d}{da}|_{a=0}\left((x\mathrm{\Phi },U(a)\mathrm{\Psi })(\mathrm{\Phi },(x+a)U(a)\mathrm{\Psi })\right)$$

$$=\frac{1}{i}\frac{d}{da}|_{a=0}(a)(\mathrm{\Phi },U(a)\mathrm{\Psi })=i(\mathrm{\Phi },\mathrm{\Psi }).$$

We can now state the general time-energy uncertainty relation:

> Let $`F`$ be a time translation covariant positive operator valued measure, and let $`H`$ denote the Hamiltonian. Let $`\mathrm{\Phi }`$ be a unit vector in the domain of the Hamiltonian for which the second moment of the probability measure $`d\mu (t)=(\mathrm{\Phi },F(dt)\mathrm{\Phi })`$ is finite. Then we have the uncertainty relation
>
> $$\mathrm{\Delta }_\mathrm{\Phi }(T_F)\mathrm{\Delta }_\mathrm{\Phi }(H)\frac{1}{2}$$
>
> where $`\mathrm{\Delta }_\mathrm{\Phi }(T_F)`$ is the square root of the variance of $`\mu `$ and $`\mathrm{\Delta }_\mathrm{\Phi }(H)=(H\mathrm{\Phi }^2(\mathrm{\Phi },H\mathrm{\Phi })^2)^{\frac{1}{2}}`$ is the usual energy uncertainty.

Proof: We use the covariant dilation described in Section 2. For $`\mathrm{\Phi }`$ we have

$$(\mathrm{\Phi },E(B)\mathrm{\Phi })=(\mathrm{\Phi },F(B)\mathrm{\Phi }),$$

hence $`\mathrm{\Phi }`$ is in the domain of definition of the selfadjoint operator $`T_E`$ defined by the projection valued measure $`E`$. Moreover, since the dilated time translations $`V(t)`$ restrict on $``$ to the original time translations, $`\mathrm{\Phi }`$ is also in the domain of the generator $`K`$ of $`V`$. But $`K`$ and $`T_E`$ satisfy the canonical commutation relation in the sense of quadratic forms, thus $`\mathrm{\Phi }`$ fulfils the uncertainty relations with respect to $`T_E`$ and $`K`$. The desired time-energy uncertainty relation now simply follows from the equalities

$$\mathrm{\Delta }_\mathrm{\Phi }(T_E)=\mathrm{\Delta }_\mathrm{\Phi }(T_F),\mathrm{\Delta }_\mathrm{\Phi }(K)=\mathrm{\Delta }_\mathrm{\Phi }(H).$$

It is clear that the tricky point was to find a useful representation of the Hilbert space $`𝒦`$ with which we could reduce the computation to the standard position-momentum uncertainty relation. However, we stress that this representation only plays an auxiliary rôle, no physical interpretation has to be associated with it. (An essentially equivalent derivation may already be found in .)

## 4 Uncertainty of Time

The replacement of projections by positive operators in the description of time observables leads to an intrinsic uncertainty. We will assume in this section that the Hamiltonian is positive. Under this condition, we will show that the minimal time uncertainty is inversely proportional to the expectation value of the energy.

Let $`\mathrm{\Phi }`$ be a unit vector for which the time uncertainty $`\mathrm{\Delta }_\mathrm{\Phi }(T_F)`$ and the expectation value of the Hamiltonian are finite. We use the same covariant dilation as before. Because of the positivity of $`H`$, $``$ must be contained in the spectral subspace of $`K`$ corresponding to the positive real axis. We may realize $`𝒦`$ as the space of square integrable functions $`L^2(𝐑,)`$ where $``$ is the Hilbert space which describes the multiplicity of the Schrödinger representation. $`K`$ acts by multiplication and $`T_E`$ as generator of translations. Since $`\mathrm{\Phi }`$ is in the domain of $`T_E`$, it is absolutely continuous, and since $`𝒦_+=L^2(𝐑_+,)`$, $`\mathrm{\Phi }`$ has to vanish at $`x=0`$. Hence $`\mathrm{\Phi }`$ is in the quadratic form domain $`q`$ of the operator $`\frac{d^2}{dx^2}`$ on $`𝒦_+`$ with Dirichlet boundary condition at $`x=0`$ (symbolically $`\frac{d^2}{dx^2}|_𝒟`$).

Since the problem is invariant under time shifts we may assume that the expectation value of $`T_F`$ vanishes, and to determine the infimum (over $`\mathrm{\Phi }`$) of the quantity

$$(\mathrm{\Phi },\frac{d^2}{dx^2}|_𝒟\mathrm{\Phi })(\mathrm{\Phi },x\mathrm{\Phi })^2,$$

it would be sufficient to take it over the set $`𝒮=\{\mathrm{\Phi }𝒦_+,\mathrm{\Phi }=1,\mathrm{\Phi }q(\frac{d^2}{dx^2}|_𝒟)q(x)\}`$.

We use the following relation which is valid for $`a,b>0`$,

$$ab^2=\underset{\lambda >0}{inf}\frac{4}{27\lambda ^2}(a+\lambda b)^3.$$

The relation may be verified by noting that the argument of the infimum assumes the value of the left side for $`\lambda =\frac{2a}{b}`$, hence it suffices to check the inequality

$$ab^2\frac{4}{27\lambda ^2}(a+\lambda b)^3,a,b,\lambda >0.$$

Setting $`c=\lambda b`$, we obtain the equivalent inequality

$$a(\frac{c}{2})^2(\frac{a+c}{3})^3.$$

Taking now the logarithm on both sides we find again an equivalent inequality which is a direct consequence of the concavity of the logarithm.

We therefore obtain the following relation

$$\underset{\mathrm{\Phi }𝒮}{inf}(\mathrm{\Phi },\frac{d^2}{dx^2}|_𝒟\mathrm{\Phi })(\mathrm{\Phi },x\mathrm{\Phi })^2$$

$$=\underset{\mathrm{\Phi }𝒮}{inf}\underset{\lambda >0}{inf}\frac{4}{27\lambda ^2}(\mathrm{\Phi },(\frac{d^2}{dx^2}|_𝒟+\lambda x)\mathrm{\Phi })^3.$$

We may perform on the right hand side first the infimum over $`\mathrm{\Phi }`$. We then can exploit the behaviour of the operator $`\frac{d^2}{dx^2}|_𝒟+x`$ under scale transformations. Namely, let

$$(D(\mu )\mathrm{\Phi })(x)=\mu ^{\frac{1}{2}}\mathrm{\Phi }(\mu x)$$

be the unitary scale transformations on $`𝒦_+`$. Then we have

$$D(\lambda ^{\frac{1}{3}})^1(\frac{d^2}{dx^2}|_𝒟+\lambda x)D(\lambda ^{\frac{1}{3}})=\lambda ^{\frac{2}{3}}(\frac{d^2}{dx^2}|_𝒟+x).$$

Since the set $`𝒮`$ is scale invariant, the infimum over $`\mathrm{\Phi }`$ is independent of $`\lambda `$. We thus obtain

$$\underset{\mathrm{\Phi }𝒮}{inf}(\mathrm{\Phi },\frac{d^2}{dx^2}|_𝒟\mathrm{\Phi })(\mathrm{\Phi },x\mathrm{\Phi })^2=\frac{4}{27}c^3,$$

where $`c`$ is the infimum of the spectrum of $`\frac{d^2}{dx^2}|_𝒟+x`$.

The spectrum of this operator is a pure point spectrum . Its eigenfunctions are

$$\mathrm{\Phi }_n(x)=\mathrm{A}i(x\lambda _n),$$

with eigenvalues $`\lambda _n`$ where $`\mathrm{Ai}`$ is the Airy function and $`\lambda _n`$ are its zeros. The smallest eigenvalue is $`\lambda _1=2.338`$. So we finally arrive at the uncertainty relation

$$\mathrm{\Delta }_\mathrm{\Phi }(T_F)\frac{d}{H_\mathrm{\Phi }}$$

with $`d=1.376`$.

Some comments are in order now:

* The new relation gives a rather large bound if compared to the original time-energy uncertainty, indeed we have

$$\mathrm{\Delta }_\mathrm{\Phi }(T_F)^2H^2_\mathrm{\Phi }=\mathrm{\Delta }_\mathrm{\Phi }(T_F)^2\left(H_\mathrm{\Phi }^2+\mathrm{\Delta }_\mathrm{\Phi }(H)^2\right)d^2+\frac{1}{4},$$

the exact largest lower bound of the left hand side being $`9/4`$. Let us also notice that the bound $`d`$ is universal, i.e., does not depend on the details of the Hamiltonian $`H`$.
* The stated relation is covariant, i.e., energy shifts do not change it. In case the infimum of the Hamiltonian is not zero we may change $`H`$ with $`Hinf(\sigma (H))\mathrm{𝟏}`$, where $`\sigma (A)`$ is the spectrum of the operator $`A`$ and $`\mathrm{𝟏}`$ is the unit operator on the Hilbert space.
* We have an explicit formula for the state with minimum uncertainty, namely the state $`\mathrm{\Phi }_1(x)=\text{Ai}(x\lambda _1)`$. Its shape shows how the energy spectrum has to be distributed in order to have minimal dispersion in time. (Recall that the variable $`x`$ labels the energy of the system.)
* In the light of the last remark one wonders whether it would be possible to prepare such a kind of state in a laboratory and check the relation explicitely.
* The same relation holds for the radial momentum of the system in place of $`T_F`$ and by replacing the Hamiltonian by the radius.


text:
# Competing roughening mechanisms in strained heteroepitaxy: a fast kinetic Monte Carlo study

\[

## Abstract

We study the morphological evolution of strained heteroepitaxial films using kinetic Monte Carlo simulations in two dimensions. A novel Green’s function approach, analogous to boundary integral methods, is used to calculate elastic energies efficiently. We observe island formation at low lattice misfit and high temperature that is consistent with the Asaro-Tiller-Grinfeld instability theory. At high misfit and low temperature, islands or pits form according to the nucleation theory of Tersoff and LeGoues.

\]

Coherent three-dimensional (3D) islands in strained heteroepitaxial films are of great interest because they can self-assemble as quantum dots for possible advanced optoelectronic applications . They are observed in a variety of film-substrate combinations including Ge/Si, InAs/GaAs, InAs/InP, etc. In these systems, island formation follows the Stranski-Krastanov mode. Initially, two-dimensional (2D) layer-by-layer growth leads to a flat wetting layer under stress. Beyond a threshold film thickness, 3D islands emerge on top of the wetting layer partially relieving the stress. The precise island formation mechanism is currently under intensive debate. According to the nucleation theory of Tersoff and LeGoues, the growth of stable islands requires overcoming an energy barrier associated with a critical island size . However, experiments reveal gradual development of ripples  or pre-pyramids  at the initial stage of island formation. These observations are more consistent with the Asaro-Tiller-Grinfeld (ATG) linear instability theory, which predicts that morphological perturbations at sufficiently long wavelengths grow spontaneously and steadily . It was originally proposed for smooth surfaces, but extensions to faceted ones based on non-equilibrium deposition conditions  or finite vicinality of substrates  have been suggested.

To better understand the roughening mechanism of strained layers, we have performed kinetic Monte Carlo simulations using an atomistic model . This approach is computationally very intensive but can reliably account for both lattice discreteness and non-equilibrium conditions. Previous simulations successfully demonstrated island formation in strained layers but only via the nucleation mechanism . In this work we introduce significantly more efficient algorithms. Thus we can explore a much wider range of conditions and observe a rich variety of morphologies in better general agreement with experiments. In particular, by lowering the lattice misfit and raising the temperature, the roughening mechanism crosses over from nucleation to instability controlled.

We adopt the 2D ball and spring model of heteroepitaxy defined on a square lattice first studied by Orr et al  and subsequently by Barabási  and Khor and Das Sarma . Simulations in 3D limited to submonolayer coverage were performed by Meixner et al . Our model parameters are appropriate to the widely studied Si1-xGex/Si system. We assume a substrate lattice constant $`a_s=2.715`$Å so that $`a_s^3`$ gives the correct atomic volume in crystalline silicon. The lattice constant $`a_f`$ of the film material is related to the lattice misfit $`ϵ=(a_fa_s)/a_f`$ which has a compositional dependence $`ϵ=0.04x`$. Nearest and next nearest neighboring atoms are directly connected by elastic springs with force constants $`k_N=13.85eV/a_s^2`$ and $`k_{NN}=k_N/2`$ respectively. This choice gives the correct modulus $`c_{11}`$ of silicon and a shear modulus constant along tangential and diagonal directions, despite a slight anisotropy in the Young’s modulus. The elastic couplings of adatoms with the rest of the system are weak and are completely neglected for better computational efficiency. Solid-on-solid conditions and atomic steps limited to at most two atoms high are assumed. Every topmost atom in the film can hop to a random topmost site $`s`$ columns away where $`s=\pm 1`$, $`\pm 2`$, … or $`\pm s_{max}`$ with equal probability. Previous simulations allowed only nearest neighbor hopping (i.e. $`s_{max}=1`$) . To speed up the simulations, we put $`s_{max}=8`$ or 20 respectively for $`x>0.6`$ or $`x0.6`$. These hopping ranges are much shorter than the dimensions of the relevant structures (islands or pits) on the films and we have checked that decreasing $`s_{max}`$ does not alter our results. The hopping rate $`\mathrm{\Gamma }_m`$ of a topmost atom $`m`$ follows an Arrhenius form:

$$\mathrm{\Gamma }_m=R_0\mathrm{exp}\left[\frac{n_m\gamma \mathrm{\Delta }E_mE_0}{k_BT}\right]$$
(1)
Here, $`n_m`$ is the number of nearest and next nearest neighbors of atom $`m`$. We choose a bond strength $`\gamma =0.4eV`$ which will be explained later. The energy $`\mathrm{\Delta }E_m`$ is the difference in the strain energy $`E_s`$ of the whole lattice at mechanical equilibrium when the site is occupied versus unoccupied. Finally, we put $`E_0=0.53`$eV and $`R_0=2D_0/(\sigma _sa_s)^2`$ with $`D_0=5.2\times 10^{13}\text{Å}^2s^1`$ and $`\sigma _s^2=\frac{1}{6}(s_{max}+1)(2s_{max}+1)`$. This gives the appropriate adatom diffusion coefficient for silicon (100) . Our model follows detailed balance.

The simulations involve intensive computations resulting solely from the long-range nature of elastic interactions. Practically all the CPU time is spent on the repeated calculations of $`E_s`$ which is needed to find $`\mathrm{\Delta }E_m`$ and hence $`\mathrm{\Gamma }_m`$ in Eq. (1). The elastic problem is formulated as follows. First, a flat film is homogeneously strained . This provides a convenient reference position with displacement $`\stackrel{}{u}_i=0`$ for every atom $`i`$. In general, the elastic force on atom $`i`$ by a directly connected neighbor $`j`$ is $`\stackrel{}{f}_{ij}=𝐊_{ij}(\stackrel{}{u}_i\stackrel{}{u}_j)+\stackrel{}{b}_{ij}`$ after linearization where the $`2\times 2`$ symmetric matrix $`𝐊_{ij}=k_{ij}\widehat{n}_{ij}\widehat{n}_{ij}^t`$ is the modulus tensor and $`\stackrel{}{b}_{ij}=(l_{ij}^0l_{ij})𝐊_{ij}\widehat{n}_{ij}`$ arises from the homogeneous stress in flat films. The spring constant $`k_{ij}`$ equals either $`k_N`$ or $`k_{NN}`$ for tangential or diagonal connection respectively. The unit column vector $`\widehat{n}_{ij}`$ points from the unstrained lattice position of atom $`j`$ towards that of atom $`i`$ and $`t`$ denotes transpose. Furthermore, $`l_{ij}^0`$ and $`l_{ij}`$ are respectively the natural and homogeneously strained spring lengths which follow easily from $`a_s`$ and $`ϵ`$. Mechanical equilibrium requires $`_j\stackrel{}{f}_{ij}=0`$ for each atom $`i`$. This leads to a large set of equations coupling the $`\stackrel{}{u}_i`$ of $`all`$ of the atoms. The solution then gives the elastic energy stored in every spring and hence $`E_s`$.

We now introduce a Green’s function approach for calculating $`E_s`$ efficiently requiring the explicit consideration of $`only`$ the surface atoms. It is a lattice analogue of boundary integral methods and is superior to boundary element techniques for our intrinsically discrete problem. We first derive the exact formalism. Figure 1 shows an example of a small lattice of atoms (solid circles). As a mathematical construct, we extend the lattice by adding $`ghost`$ atoms (open circles) with similar elastic properties. Unphysical couplings are hence introduced but can be exactly cancelled by applying external forces $`\stackrel{}{f}_j^e`$ and $`\stackrel{}{f}_j^{}^e`$ to every real surface atom $`j`$ and ghost surface atom $`j^{}`$ respectively with

$`\stackrel{}{f}_j^e`$ $`=`$ $`{\displaystyle \underset{j^{}}{}}(𝐊_{jj^{}}\stackrel{}{u}_j\stackrel{}{b}_{jj^{}})`$ (2)
$`\stackrel{}{f}_j^{}^e`$ $`=`$ $`{\displaystyle \underset{j}{}}𝐊_{jj^{}}\stackrel{}{u}_j`$ (3)
The summation in Eq. (2) is over each ghost atom $`j^{}`$ connected directly to the real atom $`j`$ and it is analogous in Eq. (3). It is easy to see that the real atoms are then exactly decoupled from the ghost atoms which are held precisely at their homogeneously strained positions.

A lattice Green’s function is then applied to express the displacement, $`\stackrel{}{u}_i`$, of every real surface atom $`i`$ under the influence of the external forces:

$$\stackrel{}{u}_i=\underset{j}{}𝐆_{ij}\stackrel{}{f}_j^e+\underset{j^{}}{}𝐆_{ij^{}}\stackrel{}{f}_j^{}^e$$
(4)
Note that the Green’s function $`𝐆`$ is defined for the extended lattice and is independent of the film morphology. It can thus be computed numerically once prior to the start of the simulation. It should not be confused with the half-plane Green’s function which provides simpler but only approximate results . Combining Eqs. (2-4), we arrive at the reduced set of equations

$$\stackrel{}{u}_i=\underset{jj^{}}{}[(𝐆_{ij}𝐆_{ij^{}})𝐊_{jj^{}}\stackrel{}{u}_j𝐆_{ij}\stackrel{}{b}_{jj^{}}]$$
(5)
coupling only the real surface atoms where the sum is over all pairs of directly connected real and ghost surface atoms $`j`$ and $`j^{}`$ respectively. The solution of Eq. (5) gives $`\stackrel{}{u}_i`$.

The elastic energy $`E_s`$ can then be calculated directly from

$$E_s=E_s^0\frac{1}{2}\underset{jj^{}}{}\stackrel{}{b}_{jj^{}}\stackrel{}{u}_j$$
(6)
which is derived from a simple consideration of virtual work. Here, the sum is defined similarly to that in Eq. (5) and $`E_s^0`$ equals the unrelaxed value of $`E_s`$ when $`\stackrel{}{u}_i0`$ which can be straightforwardly computed. The method summarized in Eqs. (5) and (6) is exact and practical for simulations at moderate scales.

We can go further and use a coarse-grained version of our Green’s function approach for a further boost on the computational efficiency. Finding the strain energy $`\mathrm{\Delta }E_m`$ of the atom $`m`$ needed in Eq. (1) requires calculating the strain energy $`E_s`$ of the whole lattice twice with and without the atom $`m`$. Certain fine details of the surface far away are obviously unimportant and can be neglected. Specifically, surface atoms are grouped into sets with the $`I`$th of them denoted by $`\mathrm{\Omega }_I`$. We neglect fluctuations within a set by assuming identical displacement $`\stackrel{}{u}_i\stackrel{}{u}_I`$ for each member $`i\mathrm{\Omega }_I`$. Equation (5) is then approximated as

$$\stackrel{}{u}_I=\underset{J}{}\left[\underset{j\mathrm{\Omega }_J,j^{}}{}(𝐆_{Ij}𝐆_{Ij^{}})𝐊_{jj^{}}\right]\stackrel{}{u}_J\underset{jj^{}}{}𝐆_{Ij}\stackrel{}{b}_{jj^{}}$$
(7)
where $`G_{Ij}=G_{ij}`$ with the lattice point $`i`$ at the centroid of the set $`\mathrm{\Omega }_I`$. Every atom within 3 columns from the atom $`m`$ is not coarsened and constitutes their own single-membered set. Farther away at $`r`$ columns from the atom $`m`$, sets contain atoms in neighborhoods of $`2r/3+1`$ columns wide, a form motivated by simple error analysis. We have checked numerically that a smaller degree of coarsening leads to no noticeable difference in our results.

Surface diffusion can be simulated using the hopping rates in Eq. (1) as $`\mathrm{\Delta }E_m`$ is now readily computable. We adopt an acceptance-rejection algorithm aided by tabulated values of $`\mathrm{\Delta }E_m`$ for $`5^8`$ sample surface configurations. Details will be explained elsewhere. In our main simulations, the lattice is of 1024 atoms wide following periodic boundary conditions. The substrate is 1024 monolayers (ML) thick while the extended lattice on which we compute the Green’s function includes also a film of 80ML. Fixed boundary conditions are applied to the top and bottom layers of the extended lattice.

We first simulate deposition of pure Ge film (i.e. $`x=1`$) with misfit $`ϵ=4`$% at temperature $`T=600K`$. At very high deposition rate $`R=80`$MLs-1 \[Fig. 2(a)\], we observe layer-by-layer growth. At slower deposition rate $`R=8`$MLs-1 \[Fig. 2(b)\], the film is initially flat but pits then develop. A detailed examination of the morphological evolution indicates that the pits appear rather independently and suddenly. Once created, they are immediately bounded by side-walls at an energetically favorable 45 inclination. These features strongly support the nucleation mechanism for their formation, noting that pits are energetically more favorable than islands . This pit nucleation process is similar to that in Fig. 1(d) of Ref. . At $`R=0.8`$MLs-1 \[Fig. 2(c)\], islands with 45 side-walls nucleate at very early stage before the film is sufficiently thick for pit formation. The result is analogous to those in Fig 1(a) of Ref.  and also Refs. . Further decreasing $`R`$ towards realistic values of order 0.01MLs-1 leads to similar but more widely separated islands.

Figure 3 shows results for deposition at misfit $`ϵ=2`$% with $`x=0.5`$ at $`T=1000K`$. Depending on $`R`$, we observe analogous layer-by-layer growth \[Fig. 3(a)\], layer-by-layer growth followed by roughening \[Fig. 3(b)\], and island growth \[Fig. 3(c)\]. However, the islands in Figs. 3(b)-(c) emerge gradually from ripple-like perturbations with local surface inclinations increasing steadily and relatively synchronously in agreement with experiments  and ATG instability theory . This regime has not been reported in previous atomistic simulations mainly due to inaccuracies in accounting for the long-range parts of the elastic interactions or the rather thin substrates used . Instead, it was observed in continuum simulations , which however cannot realize the nucleation mechanism.

The importance of the lattice misfit in deciding the roughening mechanism is particularly easy to understand. The nucleation of islands or pits occurs at a rate $`R_{nucl}\mathrm{exp}(cϵ^4)`$ with $`c`$ being a constant  and becomes very slow at small $`ϵ`$ . The ATG instability with a roughening rate $`R_{inst}ϵ^8`$  then dominates. For Figs. 2(b) and 3(b), we have chosen deposition rates close to the relevant roughening rates. We then observe kinetically limited wetting layers prior to roughening which is characteristic of Stranski-Krastanov growth. However, the threshold thickness depends strongly on $`R`$ contrary to experimental findings . A more realistic model in the future should include other mechanisms such as film-substrate interactions  or nonlinear elasticity  which have been argued to give a more stable wetting layer.

We have also simulated annealing of initially flat films of 30ML at $`T=1000K`$. At this high temperature, roughening is mainly due to the ATG instability, and we observe the development of ripples followed by islands. Figure 4 shows the surface profiles after the islands are fully developed. The cusps on the surfaces  are limited by either the substrate or the local step height limit imposed in our model. We have measured the island size $`l`$ from power spectra of 11 realizations of similar surfaces. We obtain $`lx^{1.8}`$ in reasonable agreement with $`x^2`$ from the instability theory  but slightly different from $`x^1`$ from experiments . The discrepancy with experiments may be improved if the compositional dependence of film properties such as bond energies and diffusion coefficients are properly considered. The island sizes in general lie within the experimental range due to our choice of the bond strength $`\gamma =0.4eV`$, which in fact is a reasonable value.

We have already discussed possible extensions of the model to include film-substrate chemical interactions, nonlinear elasticity, and composition-dependent material properties. To further enhance the morphological resemblance with experiments, one can explore more sophisticated forms of bond energies favoring 2D analogs of (105) and (113) facets .

We should note that the validity of ATG instability theory in 3D is not clear due to a singular form of the equilibrium surface energy in the presence of facets . Critical tests of the theory by 3D simulations are important, and may be feasible using our method. For example, our 2D simulation leading to Fig. 2(b) involves about $`6\times 10^7`$ hops and ran for 18 hours on a pentium 2GHz computer. Thus, extensions to 3D should be practical.

In conclusion, using accelerated algorithms which properly and efficiently account for long-range elastic interactions, we have simulated deposition and annealing of strained heteroepitaxial layers in 2D. At low misfit and high temperature, we observe ripples and subsequently gradual island formation in consistent with the ATG instability theory. At high misfit and low temperature, islands or pits are generated via the nucleation pathway. These suggest a non-trivial competition between roughening mechanisms, although reliable quantitative determination of the crossover conditions is beyond the scope of our model. The ATG instability is the most promising description of island formation in Si1-xGex films at low, and probably also at high lattice misfit . However, the nucleation mechanism applied to high misfit regimes in certain experimental situations has not been ruled out . Thus the competition of mechanisms can be important for interpreting experimental results. In our simulations, for deposition rates close to the relevant roughening rate, kinetically limited wetting layers develop prior to roughening. At lower but more realistic deposition rates, islands form at an early stage and are more widely separated. This may be related to the great variation in how closely the islands are packed under various conditions in experiments .

We thank B.G. Orr for helpful comments. CHL is supported by PolyU grant no. 5309/01P.


text:
# Bicocca-FT-02-17July 2002 HEAVY FLAVOUR PRODUCTION11footnote 1Invited talk at the 6-th International Workshop on Heavy Quarks and Leptons, May 27/June 1, 2002, Vietri s/m, Salerno, Italy.

## 1 Introduction

The phenomenology of heavy flavour production has been so far a mixed success. Although the order of magnitude of the total cross sections, and the shape of differential distribution is reasonably predicted, large discrepancy are present, especially for $`b`$ production.

The theoretical framework for the description of heavy flavour production is the QCD improved parton model. Besides the well-established NLO corrections to the inclusive production of heavy quark in hadron-hadron , hadron-photon , and photon-photon collisions , much theoretical work has been done in the resummation of contributions enhanced in certain regions of phase space: the Sudakov region, the large transverse momentum region and the small-$`x`$ region.

The theoretical effort involved is justified by the large variety of applications that heavy quark production physics has, in top, bottom and charm production. Besides the need of modeling these processes, heavy quark production is an important benchmark for testing QCD and parton model ideas, due to the relative complexity of the production process, the large range of masses available, and the existence of different production environments, like $`e^+e^{}`$ annihilation, hadron, photon-hadron and photon-photon collisions.

## 2 Total cross section for top and bottom

Top production  has been a most remarkable success of the theoretical model. The measured cross section has been found in good agreement with theoretical calculations, as shown in fig. 1.

Resummation of soft gluon effects  reduces the theoretical uncertainty in the cross section, pushing it toward the high side of the theoretical band. It remains, however, inside the theoretical band of the fixed order calculation, thus showing consistency with the estimated error.

Recently, the HERA-B experiment has measured the $`b\overline{b}`$ total cross section . This experiment is sensitive to the moderate transverse momentum region, where the bulk of the total cross section is concentrated. Since the production is (in a certain sense) close to threshold, resummation of Sudakov effects is useful also in this case.

A comparison of the HERA-b measurement with a theoretical calculation is shown in fig. 2. The HERA-b result is compatible with the central value prediction, with the $`b`$ pole mass around $`4.75`$ GeV. Higher precision may constrain further the $`b`$ quark mass.

## 3 Differential distributions

CDF has a longstanding disagreement with QCD in $`b`$ production. A very recent publication of the $`B^+`$ differential cross section  has quantified the disagreement as a factor of $`2.9\pm 0.2\pm 0.4`$ in the ratio of the measured cross section over the theoretical prediction. This discrepancy has been present since a long time, and it has been observed both in CDF and D0. Some authors  have argued that the discrepancy could be interpreted as a signal for Supersimmetry.

Because of the large theoretical uncertainties, this discrepancy has been often downplayed. In fact, several effects may conspire to raise the $`b`$ cross section to an appropriate value. It was early recognized that small-$`x`$ effects may be important in $`b`$ production at the Tevatron . Resummation of these effects  leads to an increase of the cross section of the order of 30%. Threshold effects are small in this case , being below 15%, but they are nevertheless positive. Resummation of large $`\mathrm{log}p_\mathrm{T}`$ yields an increase in the $`p_\mathrm{T}`$ spectrum of the order of 20% in the intermediate $`p_\mathrm{T}`$ region . The presence of many possible enhancements should not lead, however, to excessive optimism. First of all, it is not clear whether these effect can be added up without overcounting. Furthermore, they are all higher order effects, and thus should not push the cross section too far out of the theoretical band, which includes estimates of unknown higher order effects.

It has been observed since some time that an improper understanding of fragmentation effects may be one of the causes of the Tevatron discrepancy. This possibility stems from the fact that $`b`$ quark cross sections are in reasonable agreement with the Tevatron measurements of the $`B`$ meson spectrum, while the cross section obtained by applying a standard fragmentation function of the Peterson form  with $`ϵ=0.006`$ to the quark cross section yields too soft a spectrum, as can be seen in fig. 3.

In ref.  it was suggested to study $`b`$ quark jets rather than $`B`$ meson’s distributions. In fact, the jet momentum should be less sensitive to fragmentation effects than the hadron momentum. A D0 study  has demonstrated that by considering $`b`$ jets instead of $`B`$ hadrons the agreement between theory and data improves considerably (see fig. 4).

## 4 Fragmentation

If fragmentation effects are so important, they should be assessed very carefully before claiming a discrepancy with data. Non-perturbative fragmentation is usually introduced by writing the hadron-level cross section for $`B`$ mesons as

$$\frac{d\sigma ^B}{dp_\mathrm{T}}=𝑑\widehat{p}_\mathrm{T}𝑑z\frac{d\sigma ^b}{d\widehat{p}_\mathrm{T}}D(z)\delta (p_\mathrm{T}z\widehat{p}_\mathrm{T}),$$
(1)
where $`p_\mathrm{T}`$ is the hadron, and $`\widehat{p}_\mathrm{T}`$ is the quark transverse momentum. Traditionally, $`D(z)`$ is given by the Peterson form

$$D(z)\frac{z(1z)^2}{((1z)^2+ϵz)^2}.$$
(2)
Strictly speaking, this procedure is justifiable only in the limit of large transverse momenta (i.e. $`p_Tm_b`$). It can however be consider a rough model of the effect of fragmentation at moderate and small transverse momenta, as long as the variable $`z`$ is referred to the momentum or the kinetic energy, rather than the energy or the $`+`$ component, of the quark momentum.

Eq. (2) is often used with $`ϵ=0.006`$ for bottom and $`ϵ=0.06`$ for charm, in conjunction with shower Montecarlo programs. These values of $`ϵ`$ were obtained from fits to inclusive heavy flavour production distributions in $`e^+e^{}`$ collisions. The shower Montecarlo accounts for hard fragmentation effects (i.e. hard gluon radiation), where hard means from the typical transverse momentum of the process down to some cutoff, of the order of the hadron mass. This accounts for scaling violation in the inclusive cross section for the production heavy quarks222The inclusive cross section for single particle production in $`e^+e^{}`$ collisions is sometimes also called fragmentation function, which is often a source of confusion.

It has become common practice to adopt the value $`ϵ=0.006`$ also in conjunction with NLO calculations of $`b`$ production cross sections. This practice is not totally justified, since the procedure used to obtain this value of $`ϵ`$ does not match the accuracy of the NLO calculation. A better approach is possible; one can use NLL resummation of transverse momentum logarithms in the context of $`e^+e^{}`$ data in order to extract an appropriate non-perturbative fragmentation function, and supplement NLO heavy flavour production calculations with NLL resummation, in order to correctly account for scaling violation. The theoretical tools to do this are in fact available in the literature:

* Single inclusive particle production in hadronic collisions . Single hadron production are described in term of NLO single parton cross section convoluted with a NLL fragmentation function;
* Heavy quark Fragmentation Function ; a method for the computation of the heavy quark fragmentation function at all orders in perturbation theory is developed, and applied at NLL. Several applications in $`e^+e^{}`$ physics have appeared .
* Single inclusive heavy quark production at large $`p_T`$ ; item 1 and 2 are combined to give a NLL resummation of transverse momentum logarithms in heavy quark production;
* FONLL calculation of single inclusive heavy quark production; item 3 is merged without overcounting with standard NLO calculations. This procedure has been implemented both in hadroproduction  and in photoproduction .

When using these ingredients, it was found that the CDF discrepancy is considerably reduced . It can be quantified as a factor of $`1.7\pm 0.5\text{(theory)}\pm 0.5\text{(expt)}`$.

Reduction in the discrepancy is due to four basic points:

* FONLL calculation brings about a 20% increase in the intermediate $`p_T`$ region.
* The fragmentation step for going from a “perturbative” $`b`$ quark to a $`B`$ hadron seems to be too strong at small $`p_T`$ in the CDF implementation. We get a 20% increase in the small $`p_T`$ region.
* A Peterson fragmentation function with $`ϵ=0.006`$ is too soft (present $`e^+e^{}`$ data favour values around $`0.002`$); this is a 20% effect.
* A more accurate use of $`e^+e^{}`$ input data on the fragmentation function, described in the following section, brings about another 20% effect.

The net effect is $`1.2^3=1.7`$ and $`2.9/1.7=1.7`$, so that the 2.9 CDF discrepancy becomes our $`1.7`$.

## 5 A better use of $`e^+e^{}`$ data on $`b`$ fragmentation

$`D(z)`$ is extracted from fits to $`e^+e^{}`$ data. One has

$$\underset{\text{measured}}{\underset{}{\frac{d\sigma (e^+e^{}B+X)}{dp}}}=𝑑\widehat{p}𝑑z\underset{\text{computed}}{\underset{}{\frac{d\sigma (e^+e^{}b+X)}{d\widehat{p}}}}\underset{\text{fitted}}{\underset{}{D(z)}}\delta (pz\widehat{p})$$
(3)
From the above equation it is easy to show that

$$x_B^{N1}=x_b^{N1}D_N$$
(4)
where

$$x_B=\frac{p}{p_{\mathrm{max}}},x_b=\frac{\widehat{p}}{\widehat{p}_{\mathrm{max}}},D_N=𝑑zz^{N1}D(z).$$
(5)
It turns out that for the computation of the hadronic cross sections, only the first few moments of the fragmentation function are important. This follows from the fact that the heavy quark hadroproduction cross section is a steeply falling function of the transverse momentum, following roughly a power law. Assuming the form

$$\frac{d\sigma (H_1H_2b+X)}{d\widehat{p}_\mathrm{T}}A\widehat{p}_\mathrm{T}^n$$
(6)
one finds

$$\frac{d\sigma (H_1H_2b+X)}{dp_\mathrm{T}}𝑑z𝑑\widehat{p}_\mathrm{T}D(z)\frac{A}{\widehat{p}_\mathrm{T}^n}\delta (p_\mathrm{T}z\widehat{p}_\mathrm{T})=\frac{A}{p_\mathrm{T}^n}D_n=\frac{d\sigma }{d\widehat{p}_\mathrm{T}}D_n.$$
(7)
This simple fact was noticed long ago . The exponent $`n`$ ranges from 3 to 5 in hadronic collisions. In ref.  it was shown that eq. (7) is also an excellent approximation to the exact convolution.

When fitting the non-perturbative fragmentation function from $`e^+e^{}`$ data, the low moments relevant for hadroproduction are not well fitted. The Peterson form with $`ϵ=0.002`$ gives a good fit to the shape of the inclusive cross section, but not to the moments. For the purpose of predicting hadronic cross section it is instead better to make sure that the first few moments are accurately fitted. In fig. 7, it is shown that a fit to the second moment alone with the simple form $`x(1x)^\beta `$, yields a good fit to all moments in the interesting region (i.e. between 3 and 5). We also notice that the second moment fit ($`N=2`$ fit from now on), the $`ϵ=0.002`$ and the $`ϵ=0.006`$ fits differ by 20% around the $`N=4`$ moment.

The reason why $`x`$ fits do not yield good fits to moments is easy to explain. The large $`x`$ region in the inclusive production of heavy flavoured hadrons in $`e^+e^{}`$ annihilation is a very difficult region to model. It is in fact affected by large perturbative (i.e. Sudakov) and non-perturbative effects. It is unlikely that the simple forms used commonly for the non-perturbative fragmentation function can accurately model these complex effects. Thus, the large $`x`$ region is usually excluded from the fits. On the other hand, this region is the most important one for moments.

From fig. 7, another important fact becomes apparent. When using NLL calculations for the perturbative part of the fragmentation function, the importance of the non-perturbative effects on moments is greatly reduced. This can be taken as an indication of the fact that the heavy quark fragmentation function is, to a considerable extent, perturbatively calculable.

## 6 Other $`B`$ production data

The procedure of ref.  may affect other $`b`$ production data at hadron colliders. Here I consider, as an example, the inclusive muon cross section in the central and forward region, where the muons come from the semileptonic decays of the $`b`$ flavoured hadrons. These distributions were measured by the D0 experiment, and found to be higher than theoretical prediction by up to a factor of 4 in the forward region. The result of the calculation performed with the method of ref.  is reported in fig. 8.

As one can see, the agreement with data is now quite remarkable.

## 7 Conclusions

The theory of heavy flavour production seems to give a good qualitative description of the available data. In the case of top production, the comparison between theory and experiment is satisfactory also at a quantitative level.

Recent progress has taken place in the field of $`b`$ hadroproduction. The HERA-b experiment has provided a cross section for $`b`$ production at relatively low CM energy. Some progress in understanding the role of fragmentation has considerably reduced the longstanding problem of the $`b`$ momentum spectrum at the Tevatron.

Major problems do remain in the (perhaps less developed) areas of bottom production in $`\gamma \gamma `$ and $`\gamma p`$ collisions. A discussion of these problems is given in ref. .