format: text

text:
# Stochastic Optimal Prediction with Application to Averaged Euler Equations11footnote 1This work was supported in part by the Office of Science, Office of Advanced Scientific Computing Research, Mathematical, Information, and Computational Sciences Division, Applied Mathematical Sciences Subprogram, of the U.S. Department of Energy, under Contract No. DE-AC03-76SF00098, and in part by the National Science Foundation under grant number DMS98-19074.

## 1 Introduction

Many problems in mechanics, in particular problems involving turbulence, cannot be properly resolved because the number of significant degrees of freedom is too large. The problem of making numerical predictions about the behavior of systems that have not been properly resolved has been addressed in ; theoretical results can be found in ; a general introduction to such methods can be found in . When a system is underresolved, nothing much can be said without additional information; in the papers just quoted, it is assumed that the additional information consists of an invariant measure on the space of solutions; this gives rise to an optimal Markovian, deterministic, approximation. An invariant measure constitutes additional information because everything not explicitly known is assumed to be distributed according to the invariant measure; it functions like a prior measure in Bayesian statistics . Unlike what happens in other areas of application of Bayesian statistics, nature often provides a rational choice of prior, invariant measure in the form of a canonical measure. The use of an invariant measure gives rise to approximations that are optimal in a sense that we shall specify below.

In many problems this optimal approximation is still not accurate enough, and a higher-order, stochastic, approximation can be derived. In the present paper we explain these constructions with the help of a simple example, and then apply them to the solution of Averaged Euler equations in two space dimensions. The main difficulty in higher-order prediction lies in finding estimators for stochastic processes whose temporal correlations are determined only empirically.

## 2 Optimal and stochastic prediction for Hamiltonian systems

We present properties associated with general, even infinite dimensional Hamiltonian systems, in the simple case of two oscillators, with position variables $`q_1,q_2`$ and momenta $`p_1,p_2`$, and the Hamiltonian:

$$H=H(q,p)=\frac{1}{2}\left(q_1^2+q_2^2+p_1^2+p_2^2+p_1^2p_2^2\right)$$
(1)
(the “Hald system”). The equations of the motion of the system are:

$$\begin{array}{ccc}\frac{dq_1}{dt}& =& p_1+p_1p_2^2,\\ \frac{dp_1}{dt}& =& q_1,\\ \frac{dq_2}{dt}& =& p_2+p_2p_1^2,\\ \frac{dp_2}{dt}& =& q_2.\end{array}$$
(2)
We pretend that 4 equations in 4 unknowns are too difficult to solve on available computers but that 2 equations in 2 unknowns are accessible (A more realistic situation is one where one has to solve, say, $`10^{20}`$ equations and one can afford only $`10^6`$). Alternately, suppose that for some reason at time $`t=0`$ we only have values for $`q_1,p_1`$ but not for the two other variables. The question is, how does one write equations for $`q_1,p_1`$ without computing $`q_2,p_2`$.

In a standard, Galerkin, approach, one simply sets all the uncomputed variables to zero; this results in the system:

$$\frac{dq_1}{dt}=p_1,\frac{dp_1}{dt}=q_1,$$

which is not a very good approximation.

Suppose however that although the initial conditions for oscillator 2 are unknown, we do know that they are drawn from the canonical distribution:

$`P(x_1q_1<x_1+dx_1,x_2q_2<x_2+dx_2,`$
$`y_1p_1<y_1+dy_1,y_2p_2<y_2+dy_2)=`$
$`Z^1\mathrm{exp}(H(x,y)/T)dxdy,`$ (3)
where $`P`$ is the probability of the event in parentheses, $`Z`$ is a normalization constant that ensures that the sum of all probabilities is $`1`$, $`dxdy=dx_1dx_2dy_1dy_2`$, $`T`$ is a parameter that controls the variance of the samples and is known for physical reasons as the temperature, and $`H(x,y)`$ is the Hamiltonian function (1) with $`q_1`$ replaced by $`x_1`$, $`p_1`$ replaced by $`y_1`$, etc. Equation (3) is often written in the shorter symbolic form

$$P(q,p)=Z^1\mathrm{exp}(H(q,p)/T.$$
(4)
One can readily check that this probability distribution is invariant under the flow defined by (2), i.e., if the initial data are distributed as in (4), then the solutions $`q_{1,2},p_{1,2}`$ have the same distribution at all later times. Nature likes this distribution, and reproduces it often (see any book on statistical mechanics).

Suppose now that the missing initial conditions are drawn from the canonical distribution (4) conditioned by the known information $`q_1,p_1`$, i.e.,

$$P(x_2q_2<x_2+dx_2,y_2p_2<y_2+dy_2)=Z^1\mathrm{exp}(H_{q_1,p_1}(x_2,y_2))dx_2dy_2,$$
(5)
for some $`T`$, where $`H_{q_1,p_1}`$ is $`H`$ where the values of $`x_1,y_1`$ have been given the fixed, known values of these initial conditions. In the language of Bayesian estimation , the canonical distribution (4) is a prior distribution (what we believe the distribution to be before we have any data), and the conditional distribution (5) is a posterior distribution (the prior distribution modified by what we know in a special case). Averages with respect to the conditional distribution (5) are conditional averages, and denoted by $`E[|q_1,p_1]`$ (the information after the vertical line is what we know, and the prior distribution is implied). We now approximate equations (2) by the optimal prediction (OP) equations:

$$\frac{dq_1}{dt}=E[p_1+p_1p_2^2|q_1,p_1]=p_1+p_1E[p_2^2|q_1,p_1],\frac{dp_1}{dt}=q_1;$$
(6)
(clearly $`E[p_1|p_1]=p_1`$). In our particular case,

$`E[p_2^2|q_1,p_1]`$ $`=`$ $`E[p_2^2|p_1]=`$
$`{\displaystyle }{\displaystyle }p_2^2\mathrm{exp}({\displaystyle \frac{1}{2}}(p_1^2+p_2^2+q_1^2+q_2^2+p_1^2p_2^2))dq_2dp_2/`$
$`{\displaystyle \mathrm{exp}((\mathrm{}))𝑑q_2𝑑p_2};`$
(the argument of the exponential is the same in the denominator as in the numerator). After obvious cancellations,

$$E[p_2^2|p_1]=p_2^2\mathrm{exp}(p_2^2/2p_2^2p_1^2/2)𝑑p_2/\mathrm{exp}(\mathrm{})𝑑p_2=\frac{1}{1+p_1^2},$$

(a function of $`p_1`$). A general theorem states that in the mean square with respect to the invariant canonical measure, $`E[p_2^2|p_1]`$ is the best of approximation of $`p_2^2`$ by a function of $`p_1`$. The error in this instance of OP is always smaller than in the Galerkin approximation above, though in this instance not by much (see ). One should think of the system (6) as producing the average of all solutions obtained by having initially values of $`q_1,p_1`$ and sampling the other variables from the conditioned canonical distribution. An important result due to Hald () states that first-order OP for a Hamiltonian system also forms a Hamiltonian system, with a renormalized Hamiltonian which is minus the logarithm of the original Hamiltonian averaged over all the “missing” variables.

First-order OP may be optimal in a mean square sense, but it may not be good enough in many situations, and we wish to do better. In particular, the mean solution of (2) decays, while the solution of the OP equations (6) does not.

This dichotomy can be understood in several equivalent ways. From irreversible statistical mechanics we know that the canonical measure represents thermal equilibrium and that the means of all quantities tend to their equilibrium values even when they are initially conditioned by partial information; the symmetry properties of the Hamiltonian (1) ensure that the asymptotic mean is zero. We will now present a second explanation of the decay of the mean which will motivate our approach to higher-order optimal prediction. Rewrite the OP equations in the form:

$$\frac{dq_1}{dt}=p_1+p_1z(t),\frac{dp_1}{dt}=q_1,$$
(7)
(where $`z(t)`$ is of course $`p_2^2`$). In equations (6) the random function $`z(t)`$ is approximated by its conditional mean. However, in truth $`z(t)`$ varies from realization to realization of the initial data, and we are averaging over systems in which $`z(t)`$ has a mean value and a fluctuation around this mean value. If one thinks of each copy of the system, which conserves energy, as moving on some constant energy surface, the surfaces are slightly different for different copies of the system and the systems move on their surfaces at different rates. The constant energy surfaces are sphere-like, the means of these dispersed systems fall ever closer to the common center of these surfaces, which is the origin in $`qp`$ space. To capture this effect we need to take into account the variability of $`z(t)`$; the system (7) with the initial data $`q_1,p_1`$ can be viewed as a random or stochastic differential equation; the problem is that we have yet to figure out what $`z(t)`$ looks like, as a function of $`t`$ and as a random variable (the randomness coming from the initial conditions). As we now explain, the existence of an invariant measure places constraints on $`z(t)`$ but also helps in modeling it. This modeling has to rely on the specific properties of the system under consideration; some of the elegant generality of first-order optimal prediction will be lost. We shall give below an example of how a term such as $`z(t)`$ can be estimated.

## 3 The Langevin equation and fluctuation / dissipation theorems

Consider a single particle interacting with a thermal sea of other particles, the whole being presumably described in detail by some inaccessibly complicated Hamiltonian system. We wish to describe the evolution of the single particle’s velocity without explicitly calculating the evolution of the sea of particles. The Langevin equation is a standard approximation of this system:

$$\frac{du}{dt}=\gamma u+n(t),$$
(8)
where $`u=u(t)`$ is the unknown particle velocity, $`n(t)`$ is white noise with zero mean and $`\gamma `$ is a constant. The noise $`n(t)`$ repesents the fluctuating force exerted by the sea of particles, while the first term on the right hand side represents the mean force exerted by the sea which opposes the motion of the particle. The Langevin equation is a standard example that shows how an invariant measure constrains the random forcing term in a stochastic or random differential equation. This equation can be solved by elementary means . If one thinks of $`u`$ as the velocity of a particle of mass $`1`$, then one should require that asymptotically, as $`t\mathrm{}`$, the distribution of $`u`$ converge to the canonical distribution with density $`Z^1\mathrm{exp}(u^2/2T)`$; this is achieved if $`\gamma =2/T`$; this is a “fluctuation/dissipation” result . One can understand it as follows: suppose $`\gamma =0`$; then the variance of $`u`$ increases with time $`t`$ (indeed it is proportional to $`t`$). If there is no white noise, the initial variability of $`u`$ decreases because the term in $`\gamma `$ is damping. When the fluctuation/dissipation relation holds, the damping and the fluctuations imposed by the white noise balance asymptotically so that the invariant distribution is reached. Application of the fluctuation/dissipation theorems allows the prediction of the magnitude and direction of the mean force from the invariant measure and the statistical properties of the random force. The formula connecting the damping $`\gamma `$ to the temperature is based on the assumption that the autocorrelation of the noise is a delta function, as it is in white noise; more general fluctuation/dissipation theorems  will not be used in the present paper.

## 4 The Averaged Euler equations

As an example of the application of these ideas, we consider the Averaged Euler equations  in two space dimensions. We choose these equations because an analysis of the three-dimensional Euler equations presupposes an extensive discussion of turbulence, while the statistical mechanics of the usual two-dimensional Euler equations leads to negative temperature states and other unusual phenomena. The two-dimensional Averaged Euler equations describe certain temporal averages of the Euler equations as well as certain viscoelastic flows, and their statistical mechanics is compatible with the machinery we have described.

We introduce the operator

$$A=(1a^2\mathrm{\Delta })^s,$$

where $`\mathrm{\Delta }`$ is the Laplace operator, $`a`$ is a real constant, and $`s`$ is a positive number. If $`s`$ is not an integer, $`A`$ is a pseudo-differential operator. The Averaged Euler equations are:

$$\frac{A𝐮}{t}+(𝐮)A𝐮+(𝐮)^TA𝐮=p,𝐮=0,$$
(9)
where $`𝐮`$ is a vector with components $`u_\alpha ,\alpha =1,2`$. As $`a0`$, these equations formally converge to the Euler equations. We consider a periodic domain, and expand the $`u_\alpha `$ in Fourier series: $`u_\alpha =\widehat{u}_{\alpha ,𝐤}e^{i𝐤𝐱}`$, where $`𝐱=(x_1,x_2),𝐤=(k_1,k_2)`$, and $``$ denotes an inner product. Substitution into equation (9) yields the following equations of motion for the Fourier coefficients:

$$\frac{d}{dt}A(𝐤)\widehat{u}_{\alpha ,𝐤}=i\underset{\beta \gamma 𝐤^{}}{}P_{\alpha \beta }(𝐤)(k_\gamma ^{}A(𝐤)u_{\gamma ,𝐤𝐤^{}}u_{\beta ,𝐤^{}}+k_\beta ^{}A(𝐤𝐤^{})u_{\gamma ,𝐤^{}}u_{\beta ,𝐤𝐤^{}})$$
(10)
where $`P_{\alpha \beta }(𝐤)=\delta _{\alpha \beta }\frac{k_\alpha k_\beta }{k^2},`$ with $`\delta `$=Kronecker delta, $`k^2=k_1^2+k_2^2`$, is the Fourier space projection on the space of divergence-free vectors ($`k_1u_{1,𝐤}+k_2u_{2,𝐤}=0,`$), $`\beta ,\gamma `$ are component indices and $`A(𝐤)=(1+a^2k^2)^s`$ is the Fourier transform of the operator $`A`$ defined above. (This is the straightforward Fourier series of the right-hand side of the projection form of equation 9, with the zero-divergence condition built-in and the pressure eliminated).

Equation (9) conserves an energy, $`(𝐮,A𝐮)`$, where the inner product is the standard $`L_2`$ inner product, and an “enstrophy”, $`(A\xi ,A\xi )`$, where $`\xi `$ is the vorticity $`\xi =\times 𝐮`$. Each of these invariants, as well as any of their linear combinations with positive coefficients, gives rise to an invariant measure with density of the form $`Z^1\mathrm{exp}(C/T)`$, where $`C`$ is a suitable linear combination. However, one can see from general considerations () that the energy is irrelevant: if there is no enstrophy in the expression for the measure the resulting measure is not ergodic, while if the enstrophy is present the energy makes little difference. Thus we consider a measure with density $`Z^1\mathrm{exp}((A\xi ,A\xi )/T)`$, in Fourier variables; (a detailed example of such a construction is given in ). The measure is carried by divergence-free vectors (all vectors whose divergence is not zero have probability 0), and up to the normalizing factor $`Z`$ has the density

$$\mathrm{exp}(k^2(1+a^2k^2)^{2s}|\widehat{𝐮}_𝐤|^2/T),$$
(11)
(For simplicity, we are assuming for the rest of the paper that $`s=1`$ in the operator $`A`$). This expression shows why the negative temperatures of the usual Euler equations do not appear here: they are necessary in the Euler case to keep the enstrophy finite ; here, for large $`k`$, $`\widehat{𝐮}_𝐤k^3`$, which makes the enstrophy $`k^2|\widehat{𝐮}|^2`$ finite even when the temperature $`T`$ is positive.

With this measure, one can see by a simple calculation that if all velocities are divergence-free, then

$$E[\widehat{u}_{\alpha ,𝐤}^{}\widehat{u}_{\beta ,𝐤^{}}]=\frac{T\delta _{𝐤,𝐤^{}}P_{\alpha \beta }(𝐤)}{k^2(1+a^2k^2)^2}.$$
(12)
We want to model the evolution of a small number of Fourier coefficients, those that satisfy $`|𝐤|_{\mathrm{}}=\mathrm{max}(|k_1|,|k_2|)m`$, which we shall call the “resolved modes”. We will refer to the remaining modes as “sampled modes”. We rewrite the evolution equations for the resolved modes as sums of terms that depend only on the resolved modes plus terms that also involve the sampled modes. For this purpose we take advantage of the fact that the right hand side of equation (10) is quadratic in $`𝐮`$ so that its Fourier transform is a convolution of the form

$$\varphi _{\alpha \beta }(𝐤,𝐤^{})\widehat{u}_{\alpha ,𝐤^{}}\widehat{u}_{\beta ,𝐤𝐤^{}},$$
(13)
where the function $`\varphi _{\alpha \beta }`$ consists of expressions that guarantee incompressibility and perform the several differentiations; it depends only on the wave numbers but not on the amplitudes of the Fourier coefficients. The terms in the evolution of the resolved modes can be divided into three groups:

1. Those where both $`|𝐤𝐤^{}|_{\mathrm{}}m`$ and $`|𝐤^{}|_{\mathrm{}}m`$ (i.e, both are in the resolved range); we denote their sum by $`G^1`$.
2. Those where one factor belongs to the resolved range and one does not; the structure of the convolution is such that the factor in the sampled range has a wave number such that $`|𝐤|_{\mathrm{}}2m`$. In other words, resolved modes cannot interact with sampled modes further away. If we can model this subset of the sampled range, we have all the input we need to follow the dynamics of the resolved range. We call their sum $`G^2`$.
3. Those where neither factor belongs to the resolved range; their sum is $`G^3`$.

Thus the Fourier-space evolution equation (10) takes the form:

$$\frac{d}{dt}\widehat{u}_{\alpha ,𝐤}=G_{\alpha ,𝐤}^1+G_{\alpha ,𝐤}^2+G_{\alpha ,𝐤}^3.$$
(14)
In first-order optimal prediction $`G^2+G^3`$ would be approximated by its conditional expectation given the values of the Fourier components in the resolved range. The Averaged Euler equations share with the usual Euler equations the remarkable property that when the measure is Gaussian and based on energy and/or enstrophy, this conditional expectation is zero provided $`\widehat{𝐮}_0=0`$. If $`\widehat{𝐮}_00`$ the conditional expectation takes a simple form, that we shall not specify here because it will not be needed.

## 5 Monte Carlo simulations

The system of equations (10) is simple enough so that we can find its mean solution when we have partial data and the remaining data are drawn from a given distribution, by sampling the distribution $`N_{ensemble}`$ times, evolving the system in time, and averaging. The system is also simple enough that we can determine the statistical properties of the modes, in particular their time covariances. The evolution equations are evolved with a fourth-order Runge-Kutta ODE solver with adaptive time-step control. A typical simulation requires $`10^210^3`$ runs for reasonable accuracy because Monte Carlo methods exhibit errors proportional to $`1/\sqrt{N_{ensemble}}`$.

Our Monte Carlo simulation lead us to the following observations. First, the simulation results are only weakly dependent on the size of the sampled region as long as it includes $`\{𝐤||𝐤|_{\mathrm{}}2m\}`$. The evidence for this observation will be presented elsewhere.

In the rest of this paper the time correlation function

$$C(\alpha ,𝐤_\mathrm{𝟏},t_1;\beta ,𝐤_\mathrm{𝟐},t_2)=<(u_{\alpha ,𝐤_\mathrm{𝟏}}(t_1)\overline{u_{\alpha ,𝐤_\mathrm{𝟏}}(t_1)})^{}(u_{\beta ,𝐤_\mathrm{𝟐}}(t_2)\overline{u_{\beta ,𝐤_\mathrm{𝟐}}(t_2)})>$$
(15)
of the sampled modes plays an important role. Our numerical simulations lead us to the following observations about the correlation function: At time $`t=0`$, when the probability distribution is known, the various Fourier components are independent. We observe that, to a very good approximation, this remains true for $`t>0`$, and the time-correlation function is diagonal in $`𝐤`$. We may therefore restrict consideration to the autocorrelation functions of the sampled modes. The spatial index structure of the autocorrelation is determined by the requirement that the velocity fields be divergence free. We make two further observations about the correlation functions:

* The autocorrelation functions are well approximated in time as Gaussians whose width is a function of $`𝐤`$.
* The peak height of the Gaussian is well approximated by the magnitude of the correlation function in the invariant measure (12). Certainly this must be true asymptotically as $`t\mathrm{}`$. Our numerical calculations indicate that this is approximately true at all times.

These observations allow us to describe the correlation functions with a single parameter $`\sigma (𝐤)`$:

$$C(\alpha ,𝐤,t_1;\beta ,𝐩,t_2)\delta _{𝐤,𝐩}\frac{T}{k^2(1+a^2k^2)^2}P_{\alpha ,\beta }e^{(t_1t_2)^2/\sigma (𝐤)^2}$$
(16)
In general, $`\sigma (𝐤)`$ will be a function of the direction and magnitude of the momentum vector $`𝐤`$, time $`t`$, as well as the initial value of the resolved components. However, in practical application, it will be advantageous to approximate $`\sigma (𝐤)`$ by a function of the magnitude of $`𝐤`$ only. A reasonable approach to this would be to Monte Carlo simulate the correlation function $`C`$ with all $`u_{\alpha ,𝐤}`$ chosen from the invariant measure (i.e. no resolved modes). This produces the correlation functions of the invariant measure, which can only be a function of $`|𝐤|`$. Figure 1 compares $`\sigma (𝐤)`$ at equilibrium, i.e., with all the modes sampled, to $`\sigma (𝐤)`$ calculated when the initial measure has specified resolved modes in the wave-number region $`[5,5]\times [5,5]`$.

Having determined the autocorrelation functions in the prior distribution, we shall now use them to approximate specific initial value problems with prescribed partial data, thus pushing the general methodology of OP to a higher order in the statistics. In particular, $`\sigma (𝐤)`$ will henceforth be modeled as a constant divided by the magnitude of $`𝐤`$, where the constant is determined by the invariant measure, as suggested by Figure 1.

## 6 Approximation of an underresolved system by a stochastic differential equation

Equation (14) shows that the resolved modes interact with the sampled modes only through the terms $`G^2`$ and $`G^3`$. We shall ignore $`G^3`$ since it is higher order in the sampled modes and hence suppressed by factors of $`1/k^3`$. In the following, we shall use $`\widehat{u}_{\alpha ,𝐤}`$ to refer only to the resolved modes, and $`\widehat{v}_{\alpha ,𝐤}`$ to refer to the sampled modes. The sampled modes may be written as a mean plus a fluctuation: $`\widehat{𝐯}_𝐤=\overline{\widehat{𝐯}_𝐤}+\delta \widehat{𝐯}_𝐤`$. Because $`G^2`$ is linear in the sampled modes, we may write it as $`L(\widehat{𝐮})(\overline{\widehat{𝐯}}+\delta \widehat{𝐯})`$. The evolution equation for the resolved modes may now be written as

$$\frac{d}{dt}\widehat{u}_{\alpha ,𝐤}=G_{\alpha ,𝐤}^1(\widehat{𝐮})+L_{\alpha ,𝐤;\beta ,𝐩}\overline{\widehat{v}_{\beta ,𝐩}}+L_{\alpha ,𝐤;\beta ,𝐩}\delta \widehat{v}_{\beta ,𝐩}$$
(17)
We will view the fluctuating part of the sampled Fourier modes, $`\delta \widehat{𝐯}`$, not as a dynamically evolving function of a random initial configuration, but rather as a random variable with specific statistical properties, just like $`z(t)`$ in section 2. By construction $`\delta \widehat{𝐯}`$ has zero mean. The results of the previous section require $`\delta \widehat{𝐯}_𝐤`$ to be uncorrelated with any mode with different $`𝐤`$, and that its autocorrelation have the form of equation (16). We shall not specify any higher order statistics.

The evolution equation above also involves the mean value $`\overline{\widehat{𝐯}}`$. We have not presented any empirical observations about $`\overline{\widehat{𝐯}}`$. However, by analogy to the Langevin equation, we can see that the average effect of the $`\overline{\widehat{𝐯}}`$ term in the evolution equation is to provide a dissipation that counteracts the fluctuating $`\delta \widehat{𝐯}`$ term. We thus write the evolution equation as

$$\frac{d}{dt}\widehat{u}_{\alpha ,𝐤}=G_{\alpha ,𝐤}^1(\widehat{𝐮})\gamma _{\alpha ,𝐤;\beta ,𝐩}\widehat{u}_{\beta ,𝐩}+L_{\alpha ,𝐤;\beta ,𝐩}\delta \widehat{v}_{\beta ,𝐩}$$
(18)
We can use the fluctuation-dissipation theorem to estimate the dissipation matrix $`\gamma `$. Because the matrix $`\gamma `$ is dominated by the diagonal elements, in numerical computations we will only use the diagonal element:

$$\gamma _{𝐤,𝐤}=\underset{𝐤=𝐩+𝐪}{}\frac{T\sqrt{\pi }\sigma (𝐤)(𝐪^{}𝐩)(A(𝐩)p^2A(𝐪)q^2)^2|\widehat{𝐮}_𝐩|^2}{q^4A(𝐪)^2p^2k^2A(𝐤)^2}$$
(19)
where $`q^{}=(q_2,q_1)`$ and q runs over the sampled modes while p runs over the resolved modes.

The new evolution equation is now a stochastic differential equation driven by random inputs $`\delta 𝐯`$. The only Fourier modes that are evolved by the equations are the resolved modes; the sampled modes are not evolved. The statistics of the sampled modes represent the effect of the unresolved modes.

The numerical algorithm described in section 5 can be modified for the stochastic differential equation. First, only the resolved modes are evolved. The Fourier amplitudes $`\widehat{u}_{\alpha ,𝐤}`$ in the sampled region are chosen randomly from a population with the correct autocorrelation function (16). This correctly models the $`L(\widehat{𝐮})\delta \widehat{𝐯}`$ term in equation (18). Third, a dissipation term (19) is introduced on the right hand side to model the $`\gamma u`$ term in (18); this extracts from the noise the main part of its mean. Unlike what we saw in the Langevin equation of section 3 which has purely additive noise, we have no guarantee here that the fluctuation/dissipation formula will take into account all of the dissipation. In the evaluation of the dissipation term, we assume that the autocorrelation of each mode is a constant times a delta function, with the constant equal to the integral of the autocorrelation. What we have done is tranform a problem with random initial data into a stochastic differential equation.

An initial test of the stochastic differential equation approach consists in computing a norm of the average solution as a function of time. We choose as norm the A-enstrophy of the solution:

$$A_enstrophy(𝐮)=\underset{𝐤}{}k^2(1+a^2k^2)^2|\widehat{𝐮}_𝐤|^2.$$
(20)
This norm gives an equal weighting to all the modes. The A-enstrophy is, of course, conserved in time in each realisation of the system, but not on the average ( see the discussion in section 2).

Figure 2 shows the decay of the A-enstrophy of the average solution. One curve shows the true result obtained from many Monte Carlo samples of the full equations. The second curve results from the stochastic differential equation. Clearly the decay characteristics are accurately modeled by the stochastic differential equation. This is an important improvement over first order optimal prediction whose evolution is governed by a renormalized Hamiltonian as noted above, and hence is not dissipative.

## 7 Conclusions

We have used the ideas of optimal prediction to reduce an underresolved problem to a stochastic differential equation. This stochastic differential equation has fewer modes than the full equation, and requires only partial information about the initial state. On the other hand, it does require prior knowledge about the statistics of the solution, as we expect in OP methods. The full power of higher-order, stochastic OP will appear when we create effective variance-reduction techniques for the stochastic differential equation. The real test of the ideas will come when we attempt to solve Euler and Navier-Stokes equations in three space dimensions; the outstanding problem is the formulation of a reasonable invariant measure in those cases. The OP approach transfers the onus of modeling turbulence from trying to guess relations between moments to trying to guess the relevant invariant measures (and awaiting a mathematical derivation of such measures). We also expect that the OP machinery will find important uses in other problems dominated by complexity, for example in molecular dynamics (see e.g.).

Note that Figure 2 illustrates a major principle that is often overlooked: an underresolved conservative system behaves on the average as a dissipative system. The importance of this fact for the understanding of turbulence cannot be overstated.

## 8 Acknowledgements

We would like to thank Prof. G.I. Barenblatt, Prof. O. Hald, Dr. A. Kast, Prof. R. Kupferman, Prof. D. Levy, Mr. P. Okunev, and Prof. B. Turkington, for helpful discussions and for providing much of the mathematical machinery upon which this work is based.


text:
# Modified Wandzura-Wilczek Relation with the Nachtmann Variable

## Abstract

If one retains $`M^2/Q^2`$ terms in the kinematics, the Nachtmann variable $`\xi `$ seems to be more appropriate to describe deep inelastic lepton-nucleon scattering. Up to the first power of $`M^2/Q^2`$, a modified Wandzura-Wilczek relation with respect to $`\xi `$ was derived. Kinematical correction factors are given as functions of $`\xi `$ and $`Q^2`$. A comparison of the modified $`g_2^{WW}(\xi )`$, and original $`g_2^{WW}(x)`$ with the most recent $`g_2`$ data is shown.

preprint:   INPP-UVA-00-02 April 15, 2000 hep-ph/0008104
The standard Bjorken variable $`x_B=Q^2/2Pq`$ is commonly used in the discussion of deep inelastic scattering (DIS). However, if one retains $`M^2/Q^2`$ terms in the kinematics another variable $`\xi `$ introduced by Nachtmann (cf. Greenberg and Bhaumik )

$$\xi =2x/(1+\sqrt{1+4M^2x^2/Q^2})$$
$`(1)`$
seems to be more appropriate to describe DIS processes.

For the deep inelastic polarized lepton-nucleon scattering, the asymmetry depends on two spin structure functions $`g_1`$ and $`g_2`$. The structure function $`g_1(x,Q^2)`$ can be interpreted as a charge-square weighted quark helicity distribution in the parton model, and the EMC  measurement led to a surprising result - the quark spins contribute a very small fraction of the spin of the proton - the so-called “ spin puzzle ”. Since then many theoretical works and experimental measurements have been done to solve this “ puzzle ” . Most recently, the second spin structure function $`g_2`$, which includes both twist-2 and twist-3 contributions, has been measured with relatively higher precision . Using the operator product expansion (OPE) approach, one can obtain a relation between $`g_1`$ and $`g_2`$

$$g_2(x,Q^2)=g_2^{WW}(x,Q^2)+\overline{g}_2(x,Q^2),$$
$`(2)`$
where

$$g_2^{WW}(x,Q^2)g_1(x,Q^2)+_x^1\frac{dy}{y}g_1(y,Q^2).$$
$`(3)`$
If the twist-3 contribution $`\overline{g}_2(x,Q^2)`$ can be neglected, then Eq.(2) reduces to so called the Wandzura-Wilczek relation. An interesting question is that how significant is the twist-3 contribution in $`g_2`$. Model predictions (for instance see ) suggest that the twist-3 contribution is not small compared to $`g_2^{WW}`$. The earlier data given by E143 Collaboration  and most recent data given by E155 Collaboration , however, seem to show that $`g_2(x,Q^2)`$ is close to $`g_2^{WW}(x,Q^2)`$ and the twist-3 part of $`g_2`$ is rather small.

In an earlier unpublished note , we obtained a modified W-W relation (see Eq.(4) below) with respect to the Nachtmann variable $`\xi `$. The result was used by E143 Collaboration . In this brief letter, we present the modified $`g_2^{WW}(\xi )`$ and compare it with most recent data. More discussions on kinematical corrections arising from the target mass effect are given.

The modified Wandzura-Wilczek relation is

$$g_2^{WW}(\xi ,Q^2)=g_1(\xi ,Q^2)+K_2(\xi ,Q^2)_\xi ^1\frac{dy}{y}\left(\frac{g_1(y,Q^2)}{K_1(y,Q^2)}6\frac{M^2y^2}{Q^2}_y^1\frac{dz}{z}\frac{g_1(z,Q^2)}{K_1^{}(z,Q^2)}\right),$$
$`(4)`$
where the kinematic factors $`K_{1,2}`$ and $`K_1^{}`$ are

$$K_1(y,Q^2)=\frac{1M^2y^2/Q^2}{(1+M^2y^2/Q^2)(1+3M^2y^2/Q^2)},$$
$`(5a)`$
$$K_2(\xi ,Q^2)=\frac{1M^2\xi ^2/Q^2}{(1+M^2\xi ^2/Q^2)^2},$$
$`(5b)`$
$$K_1^{}(z,Q^2)=\frac{1M^2z^2/Q^2}{1+M^2z^2/Q^2}.$$
$`(5c)`$
The derivation of (4) with (5a-c) is given in the appendix. Several remarks are in order.

* (i) In the large-$`Q^2`$ limit, all correction factors $`K_1(\xi ,Q^2)`$, $`K_1^{}(\xi ,Q^2)`$ and $`K_2(\xi ,Q^2)`$ given in (5a), (5b), and (5c) approach unity and (4) becomes

$$g_2^{WW}(\xi ,Q^2)=g_1(\xi ,Q^2)+_\xi ^1\frac{dy}{y}\left(g_1(y,Q^2)6\frac{M^2y^2}{Q^2}_y^1\frac{dz}{z}g_1(z,Q^2)\right)$$
$`(6)`$
Considering $`M^2/Q^20`$, and $`\xi x`$, Eq.(6) reduces to the original W-W relation Eq.(3).
* (ii) From Eq.(1), one would expect $`\xi _{min}=0`$ and $`\xi _{max}=2/(1+\sqrt{1+4M^2/Q^2})`$ for $`x=01`$. However, since the true momentum fraction carried by quarks is $`\xi `$ (if quark is massless) rather than $`x`$, hence we should take $`\xi _{max}=1`$ from the beginning. We note that the derivation of Eq.(4) does not depend on the value of $`\xi _{max}`$.
* (iii) To show the correction effect, we first plot $`K_1`$, $`K_2`$ and $`K_1^{}`$ as functions of $`\xi `$ for $`Q^2=3`$ (GeV)2, in Fig.1. One can see that for $`\xi =01`$, $`K_1=10.324`$, $`K_2=10.460`$ and $`K_1^{}=10.581`$. It seems that all correction factors reach the maximum at $`\xi =0`$. However, their combined effect presented in Eq.(4) is not so simple.
* (iv) To show some aspects of the correction effect, we assume that $`g_1`$ in the integral in Eq.(4) is a constant and define a ratio

$$R(\xi ,Q^2)=I(K_1,K_2,\xi ,Q^2)/I(1,1,\xi ,Q^2),$$
$`(7a)`$
where

$$I(K_1,K_2,\xi ,Q^2)K_2(\xi ,Q^2)_\xi ^1\frac{dy}{y}\left(\frac{1}{K_1(y,Q^2)}6\frac{M^2y^2}{Q^2}_y^1\frac{dz}{z}\frac{1}{K_1^{}(z,Q^2)}\right).$$
$`(7b)`$
The ratios $`R(\xi ,Q^2)`$ for $`Q^2=3`$, 5, 10 and 100 (GeV/c)2 as functions of $`\xi `$ are shown in Fig.2. One can see that the ratio is quite large at low $`Q^2`$ and approaches unity when $`Q^2\mathrm{}`$. However, the ratio (7a) only provides an incomplete information of the kinematical target mass correction to $`g_2^{WW}`$. First, the function $`g_1(y,Q^2)`$ is not a constant but function of $`y`$, and secondly, one should take the whole result from Eq.(4), not just the second term.
* (v) As pointed out in  that the original Wandzura-Wilczek relation was derived from the Dirac equation for free massless quarks and no higher twist corrections were included. By using the equation of motion with nonzero quark mass and imposing the gauge invariance, an improved Wandzura-Wilczek relation is obtained in

$$g_2(x)=g_1(x)+_x^1\frac{dy}{y}g_1(y)\frac{m_q}{M}_x^1\frac{dy}{y}\frac{h_T(y)}{y}_x^1\frac{dy}{y}\mathrm{\Gamma }(y),$$
$`(8)`$
where $`h_T(x)`$ is the transverse polarization density and $`\mathrm{\Gamma }(y)`$ is related to the multiparton distribution $`h_T(x,x^{})`$. The quark mass-dependent term ($`m_q/M`$) in (8) is another twist-2 piece in addition to the usual term $`g_2^{WW}(x)`$. The last term in (8) is a twist-3 term which is coming from the quark gluon interactions. Assuming the $`m_q/M`$ term and twist-3 contribution are small, we expect that a modified version of Eq.(8) with the kinematical target mass corrections would be very similar to Eq.(4).
* (vi) Making use of the phenomenologically fitted function to the $`g_1`$ data, the modified $`g_2^{WW}(\xi )`$ in Eq.(4), $`g_2^{WW}(\xi )`$ in Eq.(6), and the original $`g_2^{WW}(x)`$ in Eq.(3) are plotted as functions of $`\xi `$ in Fig.3. The data of $`g_2(x)`$ are taken from E143  and E155 . From Fig.3, one can see that the effect of the kinematical target mass corrections is rather small relative to the experimental errors. All three $`g_2^{WW}`$ curves seems to be consistent with the $`g_2`$ data. However, more precise data are needed for a significant comparison of $`g_2^{WW}`$ and $`g_2`$. Most recently, two papers  published on the same topic - target mass corrections on the Wandzura and Wilczek relation - which found that target mass corrections do not affect the W-W relation (2) if all powers in $`M^2/Q^2`$ are included. We do not know, however, if this conclusion holds for relation (8). Anyway, since our result (4) holds up to the first power of $`M^2/Q^2`$, and target mass corrections have very small effect on the W-W relation, hence there is no contradiction between ours and theirs.
* (vii) It is easy to verify that by changing variable $`\xi `$ to $`x`$ and defining $`a(x,Q^2)\sqrt{1+4M^2x^2/Q^2}1`$, Eq.(4) can be rewritten as

$$g_2^{WW}(x,Q^2)=g_1(x,Q^2)+\frac{1+a(x,Q^2)/2}{(1+a(x,Q^2))^2}_x^1\frac{dy}{y}[\frac{1+2a(y,Q^2)}{1+a(y,Q^2)/2}g_1(y,Q^2)$$

$$\frac{y^2}{(1+a(y))(1+a(y,Q^2)/2)^2}_y^1\frac{dz}{z^3}3a(z)(1+a(z)/2)g_1(z,Q^2)].$$
$`(9)`$
This is the result obtained in .

Acknowledgments

I would like to thank P. K. Kabir for helpful comments. The author thanks L. Stuart for giving Eq.(9) and T. Averett for providing the earlier E143 data when  was completed. The author thanks J. Bluemlein and G. Ridolfi for their useful comments. This work was supported by the Institute of Nuclear and Particle Physics, Department of Physics, University of Virginia, and the Commonwealth of Virginia.

Appendix

From Eqs. (47) and (48) in Wandzura’s paper , for $`n`$=2,4,…, we have

$$_0^{\xi _{max}}𝑑\xi \xi ^n(1+\frac{M^2\xi ^2}{Q^2})[\frac{n}{n+1}(1+\frac{n+2}{n+3}ϵ)g_1+(1+ϵ+\frac{ϵ^2}{n+4})g_2^{WW}]=0,$$
$`(\mathrm{I}.1)`$
where

$$ϵϵ(\xi ,Q^2)(\frac{2M^2\xi ^2}{Q^2})/(1\frac{M^2\xi ^2}{Q^2}).$$

In obtaining Eq.(I.1), $`g_2`$ has been replaced by $`g_2^{WW}`$, or equivalently, $`\overline{g}_2`$ has been neglected.

In the large-$`Q^2`$ limit, $`ϵ0`$, (I.1) becomes

$$_0^1𝑑\xi \xi ^n[\frac{n}{n+1}g_1(\xi ,Q^2)+g_2^{WW}(\xi ,Q^2)]=0(\mathrm{large}Q^2\mathrm{limit}).$$
$`(\mathrm{I}.2\mathrm{a})`$
From (I.2a), one easily obtains

$$g_2^{WW}(\xi ,Q^2)=g_1(\xi ,Q^2)+_\xi ^1\frac{dy}{y}g_1(y,Q^2)(\mathrm{large}Q^2\mathrm{limit}),$$
$`(\mathrm{I}.2\mathrm{b})`$
which is the same form as Eq.(3), but with respect to the variable $`\xi `$. Considering $`\xi x`$ in the large-$`Q^2`$ limit, the relation (I.2b) approaches the original Wandzura-Wilczek relation Eq.(3).

On the other hand, in the large-$`n`$ limit, (I.1) becomes

$$_0^{\xi _{max}}𝑑\xi \xi ^n(1+\frac{M^2\xi ^2}{Q^2})(1+ϵ(\xi ,Q^2))[g_1(\xi ,Q^2)+g_2^{WW}(\xi ,Q^2)]=0(\mathrm{large}n\mathrm{limit}).$$
$`(\mathrm{I}.3\mathrm{a})`$
The main contribution to the integral comes from the large-$`\xi `$ region due to the suppression factor $`\xi ^n`$. It implies that

$$g_2^{WW}(\xi ,Q^2)g_1(\xi ,Q^2)(\xi \xi _{max}).$$
$`(\mathrm{I}.3\mathrm{b})`$
Considering the large-$`Q^2`$ limit (I.2b) and large-$`n`$ limit (I.3b), it is naturally to assume

$$g_2^{WW}(\xi ,Q^2)=g_1(\xi ,Q^2)+K_2(\xi ,Q^2)_\xi ^{\xi _{max}}𝑑yf(y,Q^2)(\mathrm{finite}Q^2),$$
$`(\mathrm{I}.4)`$
where $`K_2(\xi ,Q^2)`$ and $`f(y,Q^2)`$ are two unknown functions to be determined and they have the following behavior in the large-$`Q^2`$ limit

$$K_2(\xi ,Q^2)1,f(y,Q^2)\frac{g_1(y,Q^2)}{y}.$$
$`(\mathrm{I}.5)`$
Contrast $`K_2(\xi ,Q^2)`$ with $`f(y,Q^2)`$, the former is a pure kinematical correction factor and does not depend on $`g_1`$. Substituting (I.4) into (I.1) and neglecting the $`ϵ^2`$ term, we obtain

$$_0^{\xi _{max}}𝑑\xi \xi ^n(1+\frac{M^2\xi ^2}{Q^2})[(1+\frac{2n+3}{n+3}ϵ)\frac{g_1}{n+1}+(1+ϵ)K_2(\xi ,Q^2)_\xi ^{\xi _{max}}𝑑yf(y,Q^2)]=0.$$
$`(\mathrm{I}.6)`$
Unlike the derivation of (I.2b) from (I.2a), we have to use one equation, (I.6), to determine two unknown functions. Since $`K_2`$ is a pure kinematical correction factor satisfies the large-$`Q^2`$ behavior (I.5), we may choose

$$K_2(\xi ,Q^2)(1+ϵ(\xi ,Q^2))^1(1+\frac{M^2\xi ^2}{Q^2})^1,$$
$`(\mathrm{I}.7)`$
and rewrite (I.6) as

$$_0^{\xi _{max}}\frac{d\xi }{n+1}\xi ^{n+1}[f(\xi ,Q^2)(1+\frac{M^2\xi ^2}{Q^2})(1+\frac{2n+3}{n+3}ϵ)\frac{g_1(\xi ,Q^2)}{\xi }]=0,$$
$`(\mathrm{I}.8)`$
where we have exchanged the order of the integrals in the second term of Eq.(I.6).

To determine the second unknown function $`f(\xi ,Q^2)`$, we decompose it into two pieces

$$f(\xi ,Q^2)=f^{(0)}(\xi ,Q^2)+f^{(1)}(\xi ,Q^2),$$
$`(\mathrm{I}.9)`$
where $`f^{(0)}(\xi ,Q^2)O(1)`$ and $`f^{(1)}(\xi ,Q^2)O(M^2/Q^2)`$ is a small term. Since $`f^{(1)}(\xi ,Q^2)0`$ in the limit $`Q^2\mathrm{}`$, the function $`f^{(0)}(\xi ,Q^2)`$ must satisfies large-$`Q^2`$ behavior (I.5). We purposely choose

$$f^{(0)}(\xi ,Q^2)=(1+\frac{M^2\xi ^2}{Q^2})(1+2ϵ)\frac{g_1(\xi ,Q^2)}{\xi }.$$
$`(\mathrm{I}.10)`$
From (I.8), (I.9), and (I.10), we have

$$_0^{\xi _{max}}\frac{d\xi }{n+1}\xi ^{n+1}[f^{(1)}(\xi ,Q^2)+(1+\frac{M^2\xi ^2}{Q^2})\frac{3ϵ(\xi ,Q^2)}{n+3}\frac{g_1(\xi ,Q^2)}{\xi }]=0.$$
$`(\mathrm{I}.11)`$
To determine small unknown function $`f^{(1)}(\xi ,Q^2)`$, we put

$$f^{(1)}(\xi ,Q^2)\xi _\xi ^{\xi _{max}}𝑑y\eta (y,Q^2),$$
$`(\mathrm{I}.12)`$
where $`\eta (y,Q^2)`$ should be the order of $`O(M^2/Q^2)`$. Substituting (I.12) into (I.11), one obtains

$$_0^{\xi _{max}}𝑑\xi \frac{\xi ^{n+3}}{(n+1)(n+3)}[\eta (\xi ,Q^2)+3ϵ(\xi ,Q^2)(1+\frac{M^2\xi ^2}{Q^2})\frac{g_1(\xi ,Q^2)}{\xi ^3}]=0.$$
$`(\mathrm{I}.13)`$
This equation can be satisfied for all $`n`$ ($`n`$=2,4,…) only if the term in the bracket vanishes. Therefore

$$\eta (\xi ,Q^2)=3ϵ(\xi ,Q^2)(1+\frac{M^2\xi ^2}{Q^2})\frac{g_1(\xi ,Q^2)}{\xi ^3},$$
$`(\mathrm{I}.14)`$
which is indeed the order of $`O(M^2/Q^2)`$. Substituting (I.9), (I.10), (I.12), and (I.14) into (I.4), we finally obtain the modified W-W relation Eq.(4) with correction factors (5a), (5b), and (5c).


text:
# Untitled Document

Oscillations of spherical and cylindrical shells

Serkan Erdin,a) and Valery L. Pokrovsky a,b)

a) Department of Physics, Texas A&M University, College Station

TX 77843-4242

b) Landau Institute for Theoretical Physics, Moscow 117940, Russia

## Abstract

We have found the complete spectrum and eigenstates for harmonic oscillations of ideal spherical and cylindrical shells, both being infinetely thin. The spectrum of the cylindrical shell has an infinite number of Goldstone modes corresponding to folding deformations. This infrared catastrophe is overcome by accounting for curvature-dependent part of energy.

INTRODUCTION

The oscillating surfaces have numerous applications in many areas of physics, spanning from biophysics and chemical physics to high energy physics. 1,2,3In this article, we study the oscillation spectra of spherical and cylindrical surfaces, considering them as infinitely thin elastic membranes. We also calculate their heat capacities. Corresponding calculated values of heat capacity are valid at sufficiently low temperature $`T<<\mathrm{\Theta }_D`$, where $`\mathrm{\Theta }_D`$ is Debye temperature, otherwise the effects of discreteness become essential. To our knowledge and big surprise, this problem has not yet been solved, despite of a vast literature on the subject. 4

THE SPECTRUM AND EIGENSTATES OF SPHERICAL AND CYLINDRICAL SHELLS

A. Spherical Membrane

The Lagrangian of a spherical membrane reads :

$``$ $`=`$ $`{\displaystyle \frac{\rho }{2}}{\displaystyle R^2𝑑\mathrm{\Omega }\left[\left(\frac{u_r}{t}\right)^2+\left(\frac{u_\theta }{t}\right)^2+\left(\frac{u_\varphi }{t}\right)^2\right]}`$ (1)
$``$ $`{\displaystyle \frac{\lambda }{2}}{\displaystyle R^2𝑑\mathrm{\Omega }(U_{\theta \theta }+U_{\varphi \varphi })^2}\mu {\displaystyle R^2𝑑\mathrm{\Omega }(U_{\theta \theta }^2+2U_{\theta \varphi }^2+U_{\varphi \varphi }^2)}`$ (2)
where $`R`$ is radius of spherical membrane, $`d\mathrm{\Omega }=\mathrm{sin}\theta d\theta d\varphi `$ and, $`\mu `$ and $`\lambda `$ are the Lame coefficients. $`u_r`$, $`u_\theta `$, $`u_\varphi `$ are components of the displacement vector and $`U_{\theta \theta }`$, $`U_{\varphi \varphi }`$, $`U_{\theta \varphi }`$ are components of the strain tensor. The latter are determined by relative change of distance between two points at a deformation : $`ds^{}_{}{}^{}2=ds^2(1+U_{\alpha \beta }dx_\alpha dx_\beta )`$ and can be expressed in terms of the displacement vector and its derivatives as follows: 5

$$U_{\theta \theta }=\frac{1}{R}\frac{u_\theta }{\theta }+\frac{u_r}{R}$$
(3)
$$U_{\varphi \varphi }=\frac{1}{R\mathrm{sin}\theta }\frac{u_\varphi }{\varphi }+\frac{u_\theta }{R}\mathrm{cot}\theta +\frac{u_r}{R}$$
(4)
$$2U_{\theta \varphi }=\frac{1}{R}\left(\frac{u_\varphi }{\theta }u_\varphi \mathrm{cot}\theta \right)+\frac{1}{R\mathrm{sin}\theta }\frac{u_\theta }{\varphi }$$
(5)
Equations of motions for oscillations with the frequency $`\omega `$ read :

$$u_r=𝒦\left(\frac{u_\theta }{\theta }+u_\theta \mathrm{cot}\theta +\frac{1}{R\mathrm{sin}\theta }\frac{u_\varphi }{\varphi }\right)$$
(6)
$`\left(\lambda +2\mu +2(\lambda +\mu )𝒦\right)\left({\displaystyle \frac{^2u_\theta }{\theta ^2}}+{\displaystyle \frac{(u_\theta \mathrm{cot}\theta )}{\theta }}\right)+\left(\lambda +\mu +2(\lambda +\mu )𝒦\right){\displaystyle \frac{1}{\mathrm{sin}\theta }}{\displaystyle \frac{^2u_\varphi }{\theta \varphi }}`$ (7)
$`\left(\lambda +3\mu +2(\lambda +\mu )𝒦\right){\displaystyle \frac{\mathrm{cot}\theta }{\mathrm{sin}\theta }}{\displaystyle \frac{u_\varphi }{\theta }}+{\displaystyle \frac{\mu }{\mathrm{sin}^2\theta }}{\displaystyle \frac{^2u_\theta }{\varphi ^2}}+(\rho \omega ^2R^2+2\mu )u_\theta =0`$ (8)
$`\left(\lambda +\mu +2(\lambda +\mu )𝒦\right)\left({\displaystyle \frac{1}{\mathrm{sin}\theta }}{\displaystyle \frac{^2u_\theta }{\theta \varphi }}\right)+\left(\lambda +3\mu +2(\lambda +\mu )𝒦\right)\left({\displaystyle \frac{\mathrm{cot}\theta }{\mathrm{sin}\theta }}{\displaystyle \frac{u_\theta }{\varphi }}\right)`$ (9)
$`+\left(\lambda +2\mu +2(\lambda +\mu )𝒦\right)\left({\displaystyle \frac{1}{\mathrm{sin}^2\theta }}{\displaystyle \frac{^2u_\theta }{\varphi ^2}}\right)+\mu \left({\displaystyle \frac{^2u_\varphi }{\theta ^2}}+{\displaystyle \frac{(u_\varphi \mathrm{cot}\theta )}{\theta }}\right)`$ (10)
$`+(\rho \omega ^2R^2+2\mu )u_\varphi =0`$ (11)
where $`𝒦=\frac{2(\lambda +\mu )}{\rho \omega ^2R^24(\lambda +\mu )}`$. There exists a symmetric solution (the ”breathing” mode), $`u_\theta =u_\varphi =0`$ , $`u_r=const0`$. The frequency of this solution is $`\omega _s=\sqrt{\frac{4(\lambda +\mu )}{\rho R^2}}`$. From other solutions, we first consider those independent on $`\varphi `$. There exist three branches of such oscillations. The first two (”longitudinal” modes) correspond to $`u_\varphi =0`$, $`u_\theta 0`$, $`u_r0`$. The frequency of these oscillations are labeled by an integer $`l`$ :

$`\omega _{\pm ,l}^2`$ $`=`$ $`{\displaystyle \frac{[4\lambda +2\mu +(\lambda +2\mu )l(l+1)]}{2\rho R^2}}`$ (12)
$`\pm `$ $`{\displaystyle \frac{\sqrt{(4\lambda +2\mu +(\lambda +2\mu )l(l+1))^216\mu (\lambda +\mu )(l(l+1)2)}}{2\rho R^2}}`$ (13)
In the third branch (a ”shear” mode) $`u_\theta =0`$ , $`u_r=0`$ and $`u_\varphi 0`$. Its frequencies are

$$\omega _{t,l}^2=\frac{\mu }{\rho R^2}(l(l+1)2)$$
(14)
where $`l=2,3,\mathrm{}`$ Explicit expressions for the displacements for the longitudinal branches are: 6

$$u_\varphi =0u_\theta =C\frac{dP_l(\mathrm{cos}\theta )}{d\theta }$$
(15)
$$u_r=\frac{C(\lambda +\mu )}{\rho \omega ^2R^24(\lambda +\mu )}l(l+1)P_l(\mathrm{cos}\theta )$$
(16)
where $`C`$ is constant amplitudes and $`P_l(\mathrm{cos}\theta )`$ are the Legendre polynomials. The displacements in the shear branch solution read :

$$u_r=u_\theta =0u_\varphi =C\frac{dP_l(\mathrm{cos}\theta )}{d\theta }$$
(17)
In order to get eigenstates for $`m0`$, we employ an infinitesimal rotation around $`x`$ axis which turns an arbitrary vector $`\stackrel{}{r}`$ into :

$$\stackrel{}{r}^{^{}}=O(\omega )\stackrel{}{r}=1+\omega g_x\stackrel{}{r}$$
(18)
Here $`ig_x`$ is the generator of the rotation around $`x`$ axis. 7 It generates a transformation of the displacement field $`\stackrel{}{u}(\stackrel{}{r})`$.

$$\stackrel{}{u}^{^{}}(\stackrel{}{r})=O(\omega )\stackrel{}{u}\left(O^1(\omega )\stackrel{}{r}\right)=\stackrel{}{u}(\stackrel{}{r})\omega g_x\stackrel{}{u}(\stackrel{}{r})+\omega (g_x\stackrel{}{r}).u(\stackrel{}{r})\mathrm{}$$
(19)
If $`u_0(\stackrel{}{r})`$ was a solution of the system (5-7), then $`\stackrel{}{u}^{^{}}(\stackrel{}{r})`$ is also a solution. The operation we performed to find new solutions is readily recognizable as the Lee translation introduced in theory of fiber bundles. 8

After a simple transformation, we find an operator $`𝒢`$, which, acting on a solution fixed by two components $`\begin{array}{cccc}& & & \\ & (\mathrm{}& u_\theta & )\mathrm{}\\ & u_\varphi \end{array}`$ , generates another solution with the same frequency.

$`𝒢=\begin{array}{cccc}& & & \\ & (\mathrm{}& \mathrm{sin}\varphi \frac{}{\theta }+\mathrm{cot}\theta \mathrm{sin}\varphi \frac{}{\varphi }& \frac{\mathrm{cos}\varphi }{\mathrm{sin}\theta }& )\mathrm{}\\ & \frac{\mathrm{cos}\varphi }{\mathrm{sin}\theta }& \mathrm{sin}\varphi \frac{}{\theta }+\mathrm{cot}\theta \mathrm{sin}\varphi \frac{}{\varphi }\end{array}`$ (20)
The operator $`𝒢`$ is the analog of the operator of total angular momentum $`𝒥`$ in quantum mechanics, but it is reduced to a subspace of tangential components of the displacement vector $`\stackrel{}{u}`$. The Lagrangian (1) and equations of motion (6-7) are invariant with respect to the transformation (15). Acting by the operator $`𝒢`$ onto solutions with $`m=0`$ (10-12), we generate solutions with all possible values of $`m`$ . They belong to the same eigen-frequencies as the initial solutions with $`m=0`$. The displacements for solutions belonging to the longitudinal branches with $`m0`$ read :

$$u_\theta ^{(m)}(x)=e^{im\varphi }i^m\left[mx(1x^2)^{\frac{m1}{2}}P_l^{(m)}(x)(1x^2)^{\frac{m+1}{2}}P_l^{(m+1)}(x)\right]$$
(21)
$$u_\varphi ^{(m)}(x)=e^{im\varphi }i^{m+1}m(1x^2)^{\frac{m1}{2}}P_l^{(m)}(x)$$
(22)
The displacements for oscillations of the shear branch with $`m0`$ are :

$$u_\theta ^{(m)}(x)=e^{im\varphi }i^{m+1}m(1x^2)^{\frac{m1}{2}}P_l^{(m)}(x)$$
(23)
$$u_\varphi ^{(m)}(x)=e^{im\varphi }i^m\left[mx(1x^2)^{\frac{m1}{2}}P_l^{(m)}(x)(1x^2)^{\frac{m+1}{2}}P_l^{(m+1)}(x)\right]$$
(24)
Given the known frequency spectrum $`\omega _j`$, the thermal energy can be calculated as $`U=_j\frac{\mathrm{}\omega _j}{exp\frac{\mathrm{}\omega _j}{T}1}`$ and the heat capacity can be calculated as $`C=\left(\frac{U}{T}\right)_V`$.

As always, low-lying rotational frequencies $`\omega _k^{(rot)}=\frac{\mathrm{}k(k+1)}{2MR^2}`$ are smaller than low-lying oscillatory frequencies by a factor $`\sqrt{\frac{m}{M}}\frac{a}{R}`$, where $`a`$ is characteristic distances between atoms ( discreetness effect ) and $`m`$ is the mass of electron. That means that at temperature $`T>>T_R=\frac{\mathrm{}^2}{MR^2}`$, the rotational degrees of freedom become substantially classical until $`k\sqrt{\frac{T}{T_R}}`$.Their contribution to the heat capacity is $`C_R=\frac{T}{T_R}\zeta (3)`$. Elastic oscillations remain frozen till a larger temperature $`T_0\frac{\mathrm{}}{R}\sqrt{\frac{\mu }{\rho }}>>T_R`$.

B. Cylindrical shells

In order to find the spectrum and corresponding eigenstates of a cylindrical membrane, we follow the same steps as we did for spherical shell. For cylindrical case, the components of strain tensor in terms of cylindrical coordinates $`r`$ , $`\varphi `$, $`z`$ read: 5

$$U_{\varphi \varphi }=\frac{1}{R}\frac{u_\varphi }{\varphi }+\frac{u_r}{R},U_{zz}=\frac{u_z}{z},2U_{\varphi z}=\frac{1}{R}\frac{u_z}{\varphi }+\frac{u_\varphi }{z}$$
(25)
The Lagrangian reads :

$``$ $`=`$ $`{\displaystyle \frac{\rho }{2}}{\displaystyle R𝑑\varphi 𝑑z\left[\left(\frac{u_r}{t}\right)^2+\left(\frac{u_\varphi }{t}\right)^2+\left(\frac{u_z}{t}\right)^2\right]}`$ (26)
$``$ $`{\displaystyle \frac{\lambda }{2}}{\displaystyle R𝑑\varphi 𝑑z(U_{zz}+U_{\varphi \varphi })^2}\mu {\displaystyle R𝑑\varphi 𝑑z(U_{zz}^2+2U_{z\varphi }^2+U_{\varphi \varphi }^2)}`$ (27)
where $`R`$ is the radius of the cylinder. Equations of motions are :

$$u_r=(\lambda +2\mu )\stackrel{~}{𝒦}\frac{u_\varphi }{\varphi }+\stackrel{~}{𝒦}R\lambda \frac{u_z}{z}$$
(28)
$$u_\varphi =(\lambda +2\mu )\stackrel{~}{𝒦}\frac{^2u_\varphi }{\varphi ^2}(\lambda +\mu )\stackrel{~}{𝒦}\frac{^2u_z}{z\varphi }+\frac{\mu }{\rho \omega ^2R^2}\left((\lambda +2\mu )\stackrel{~}{𝒦}R\frac{^2u_z}{z\varphi }R^2\frac{^2u_\varphi }{z^2}\right)$$
(29)
$`u_z`$ $`=`$ $`(\lambda +\mu )\stackrel{~}{𝒦}R{\displaystyle \frac{^2u_\varphi }{z\varphi }}(\lambda +2\mu )\stackrel{~}{𝒦}R^2{\displaystyle \frac{^2u_z}{z^2}}`$ (30)
$`+`$ $`{\displaystyle \frac{\mu }{\rho \omega ^2R^2}}\left((\lambda +2\mu )\stackrel{~}{𝒦}R{\displaystyle \frac{^2u_\varphi }{z\varphi }}+4(\lambda +\mu )\stackrel{~}{𝒦}R^2{\displaystyle \frac{^2u_z}{z^2}}{\displaystyle \frac{^2u_z}{\varphi ^2}}\right)`$ (31)
where $`\stackrel{~}{𝒦}=\frac{1}{\rho R^2\omega ^2(\lambda +2\mu )}`$. Specifying the $`z`$ dependence of all displacement components as a plane wave $`e^{ikz}`$, and their angular dependence as $`e^{im\varphi }`$, we find a system of linear equations for the amplitudes $`u_r`$, $`u_\theta `$ , $`u_\varphi `$. Its secular equation reads :

$`(\omega ^2\omega _s^2)(\omega ^6(\omega _t^2(R^2k^2+m^2)+\omega _s^2(R^2k^2+m^2+1))\omega ^4`$ (32)
$`+\left(\omega _t^2\omega _s^2\left((R^2k^2+m^2)^2+5R^2k^2+m^2\right)4\omega _t^4R^2k^2\right)\omega ^2`$ (33)
$`+4\omega _t^4(\omega _t^2\omega _s^2)R^4k^4)=0`$ (34)
where $`\omega _s^2=\frac{\lambda +2\mu }{\rho R^2}`$ and $`\omega _t^2=\frac{\mu }{\rho R^2}`$. One of its solutions (symmetric) is $`\omega ^2=\omega _s^2`$. For $`k=0`$, the remaining cubic equation can be solved explicitly for arbitrary $`m`$. The three solutions are :

$$\omega ^2=0,\omega ^2=(m^2+1)\omega _s^2,\omega ^2=m^2\omega _t^2$$
(35)
The first of them implies that there exists infinite number of Goldstone modes corresponding to the folding deformations. For $`m=0`$, the cubic equation gives the following solutions :

$$\omega ^2=k^2R^2\omega _t^2$$
(36)
$$\omega ^2=\frac{\omega _s^2(R^2k^2+1)}{2}\pm \frac{\sqrt{(\omega _s^2(R^2k^2+1))^216\omega _t^2(\omega _s^2\omega _t^2)R^2k^2}}{2}$$
(37)
The physical reason of the peculiarity in the spectrum of cylindrical shell is an opportunity to fold a cylinder conserving its straight-lines and its cross-section contour length. Such a deformation is strainless. The infinite number of folding modes leads to divergence of the heat capacity at any temperature. This difficulty can be avoided by introducing terms with higher derivatives into the Lagrangian. Geometrically, it means that the Gaussian curvature and the mean curvature must be included. In more physical or chemical terms it means that not only the variation of distances between atoms, but also variation of valence angles must be accounted for. Below we write down explicit formulae for the Gaussian curvature $`K`$ and the mean curvature $`M`$ for cylindrical membranes in linear approximation in displacements.

$$K=\frac{1}{R}\frac{^2u_r}{z^2}$$
(38)
$$M=\frac{1}{R}\frac{u_r}{R^2}\frac{1}{R^2}\frac{^2u_r}{\theta ^2}\frac{^2u_r}{z^2}$$
(39)
The term in Lagrangian depending on curvature reads :

$$\mathrm{\Delta }L=\left[\frac{\nu }{2}(KK_0)^2+\frac{\theta }{2}(MM_0)^2\right]d^2S$$
(40)
where $`K_0=0`$ and $`M_0=\frac{1}{R}`$ are the values of $`K`$ and $`M`$ for undeformed cylinder. With this term in Lagrangian we find the gaps $`\mathrm{\Delta }_m`$ in the spectrum corresponding to the former folding modes:

$$\mathrm{\Delta }_m=\sqrt{\frac{\theta }{\rho R^4}}(m^21)$$
(41)
Note that the gaps in spectrum due to curvature effect are propotional to $`R^2`$ and, therefore, are small in comparison to gaps due to strain effects ($`R^1`$). By perturbation method it is possible to find small $`k`$ \- dependent corrections to these frequencies : $`\delta \omega ^2=\frac{\theta }{\rho R^4}[2k^2R^2(m^21)]\frac{(\omega _s^22\omega _t^2)\omega _t^2(kR)^2}{\mathrm{\Delta }_m^2\omega _t^2m^2}`$ or $`\omega _m^2(k)=\mathrm{\Delta }_m^2+V_m^2k^2`$, where $`V_m^2=\frac{\theta }{\rho R^4}[2R^2(m^21)]\frac{(\omega _s^22\omega _t^2)\omega _t^2(R)^2}{\mathrm{\Delta }_m^2\omega _t^2m^2}`$. A clear tendency is that the band becomes more flat with the growth of $`m`$.

For cylindrical shell, we can write the thermal energy in the following form :

$$U=\underset{0}{\overset{\mathrm{}}{}}\frac{dk}{2\pi }\underset{m=2}{\overset{\mathrm{}}{}}\frac{\mathrm{}\omega _m(k)}{exp\frac{\mathrm{}\omega _m(k)}{T}1}$$
(42)
In the temperature interval $`\mathrm{}\mathrm{\Delta }_2<<T<<\mathrm{}\omega _t`$, we can neglect all ” elastic branches” and substitute the summation over the ”valence” branches by integration. It leads to the energy :

$$U=\frac{0.601}{\mathrm{}^2}\sqrt{\frac{1}{\frac{\theta }{\rho R^2}(\omega _s^22\omega _t^2)}}T^3$$
(43)
and the heat capacity :

$$C_v=\frac{1.803}{\mathrm{}^2}\sqrt{\frac{1}{\frac{\theta }{\rho R^2}(\omega _s^22\omega _t^2)}}T^2$$
(44)
For $`T<<\mathrm{}\mathrm{\Delta }_2`$, we can neglect the sum over $`m`$, since only $`\mathrm{\Delta }_2`$ and $`V_2`$ contribute to the energy. The internal energy becomes exponentially small :

$$U=\sqrt{2\pi \mathrm{}T}\frac{\mathrm{\Delta }_2^{\frac{3}{2}}}{V_2}e^{\frac{\mathrm{}\mathrm{\Delta }_2}{T}}$$
(45)
and the heat capacity :

$$C_v=\frac{\sqrt{2\pi \mathrm{}^3}}{T^{\frac{5}{2}}}\frac{\mathrm{\Delta }_2^{\frac{5}{2}}}{V_2}e^{\frac{\mathrm{}\mathrm{\Delta }_2}{T}}$$
(46)
CONCLUSIONS

We have found complete spectrum of oscillations for spherical and cylindrical membranes and corresponding eigenstates. We have showed that cylindrical shell has an infinite number of Goldstone modes due to folding deformations which do not change the elastic energy of shell. We have also showed that this catastrophe can be avoided by introducing curvature-dependent terms to the Hamiltonian. In addition, heat capacities of spherical and cylindrical shells are calculated. They can be measured experimentally as an additional contribution to the heat capacity of a solvent in a diluted solution of long molecules.

ACKNOWLEDGMENTS

This work was inspired by discussion of heat capacity measurements for buckyballs and nanotubes with Vladimir Kopylov.

1 Stanislas Leibler, in Statistical Mechanics of Membranes and Surfaces, edited by D.R.Nelson, T.Piran, S.Weinberg, Vol. 5(World Scientific, Singapore, 1989), p.45

2 J.M. Park and T.C. Lubensky, J.Phys.16, 493 (1996); E.I. Kats, V.V. Lebedev and A.R. Muratov, JETP Letters 63, 216 (1996); J.F. Weather, Nucl. Phys.458, 671 (1996)

3 Science of Fullerenes and Carbon Nanotubes, edited by M.S. Dresselhaus ( Acedemic, New York, 1996 )

4 P.M. Naghdi, in Mechanics of Solids , edited by C. Truesdell, Vol.2 (Springer-Verlag, Berlin, 1972), p.425

5 L.D. Landau and E.M. Lifshitz ,Theory of Elasticity (Pergenon, New York, 1970)

6 W. Soedel , Vibrations of Shells and Plates  (Marcel Dekker, New York, 1993)

7 Group Theory in Physics: Introduction, edited by J.F. Cornwel (Academic, New York, 1997)

8 B.F. Schutz, Geometrical Methods of Mathematical Physics ( Cambridge University Press, Cambridge, 1980 )


text:
# TESTS OF QCD: SUMMARY OF DIS 2000

## 1 Jets in deep-inelastic scattering

Schörner, Pöschl, Chapin, and Caron reported on measurements of inclusive jet (1+1) $`^\mathrm{?}`$, dijet (2+1) $`^{\mathrm{?},\mathrm{?}}`$, and three jet (3+1) $`^\mathrm{?}`$ production by H1 and ZEUS, scanning an enormous region of the kinematic plane and of the jet phase space at HERA. They provide therefore an important and demanding test of the adequacy of the current QCD calculations in next-to-leading order (NLO) of the strong coupling parameter in these regions.

In the photon virtuality $`Q^2`$ the measurements cover the range from 5 GeV2 to 104 GeV2 , in Bjorken-$`x`$ from 10-4 to 0.3, and in $`\xi =x(1+M_{jj}^2/Q^2)`$ from 5$`\times 10^3`$ to 0.1, where $`\xi `$ in leading order (LO) of the strong coupling parameter is the momentum fraction of the gluon or quark emitted by the proton. Jets are defined using the inclusive $`k_T`$-algorithm in the Breit frame. $`^\mathrm{?}`$ In this frame the transverse energy of a jet is dominated by QCD processes. In the jet phase space the measurements span the range in the transverse energy of the jets in the Breit frame from $`5<E_T<100`$ GeV, in $`E_T^2/Q^2`$ from $`0.5`$ to $`100`$, in rapidity from forward to central and backward jets in the HERA frame, and in the three jet invariant mass $`M_{3j}`$ up to $`100`$ GeV.

The comparison of these jet measurements to QCD calculations in NLO allows to draw the following conclusions: 1) Calculations using $`Q^2`$ as renormalization scale can describe almost all of the data at the prize of large uncertainties at low $`Q^2`$, when varying the scale by a factor of 1/4 to 4. 2) When using $`E_T^2`$ as the relevant scale, good agreement with the data is achieved for large jet $`E_T`$ ($`E_T>20`$ GeV) $`^\mathrm{?}`$ and for large $`Q^2`$ ($`Q^2>50`$ GeV2$`^{\mathrm{?},\mathrm{?}}`$. For intermediate values of the ratio $`E_T^2/Q^2`$, between 1 and 50, the calculations are below the data $`^{\mathrm{?},\mathrm{?}}`$. The forward jet cross section shown by H1 $`^\mathrm{?}`$, requiring $`x_{jet}>0.035`$ and $`1.5<\eta _{lab}<2.8`$, but without demanding that $`E_T^2Q^2`$ as done in previous such analyses by H1 and ZEUS $`^\mathrm{?}`$, rises more strongly than the NLO prediction for decreasing $`x`$ down to $`10^4`$. With $`E_T^2`$ as the renormalization scale the uncertainties in the predicted cross sections are however much smaller than for the scale $`Q^2`$. 3) Typically the discrepancies between data and NLO calculations are large, where the corrections from LO to NLO, the hadronisation corrections, and the renormalization scale uncertainties are large. 4) In a large part of the phase space the theoretical uncertainty due to the renormalization scale are much larger than the experimental errors. The latter are dominated by the uncertainty of the hadronic energy scale of the calorimeters of both H1 and ZEUS.

The NLO QCD scale uncertainty estimates for the DIS dijet cross section as a function of $`Q^2`$ and $`<E_T>^2/Q^2`$ are exemplified in Figs. 2 and 2 respectively. In Fig. 2 the scale uncertainty can be seen to grow rapidly for decreasing $`Q^2`$. In the lower $`Q^2`$ region it is much larger than the experimental error. In Fig. 2, in the region where $`E_T^2>Q^2`$ the difference between the NLO predictions using $`\mu _r^2=Q^2`$ and $`\mu _r^2=E_T^2/4`$ is larger than the respective estimated scale uncertainties. Clearly next-to-next-to leading order (NNLO) and also resummed calculations are needed to make further progress.

Forward jet cross sections at very small $`x`$ are supposed to be a clean test of BFKL or CCFM dynamics. The latter should be the way to describe low-$`x`$ final states to leading double-log accuracy. Lönnblad discussed the implementation of the CCFM equation in different Monte Carlo generators. $`^\mathrm{?}`$ Two of the programs implementing CCFM, SMALLX $`^\mathrm{?}`$ and the new complete final state Monte Carlo generator CASCADE $`^\mathrm{?}`$ are able to describe the forward jet cross sections. However they do so only, if the so-called consistency constraint $`k_{ti}^2>z_iq_{ti}^2`$, which has been found to be necessary in BFKL dynamics to account for kinematic effects, is not applied. In addition, the results are found to be sensitive to the treatment of non-singular terms in the gluon splitting function. However, no complete calculation including these terms has been done, and therefore no firm conclusions can be drawn.

## 2 Jets in $`\gamma `$p, $`\gamma ^{}`$p, and $`\gamma \gamma `$

A photon is not just a photon and not just a hadron. In LO QCD it has two components, a direct one, which couples electromagnetically to one of the partons of the other beam, and a resolved one. In the latter it fluctuates into a $`q\overline{q}`$ or even more complicated partonic state. In jet production the probing scale is given by the $`E_T^2`$ of the jets and the ‘size’ of the photon by $`Q^2`$. Therefore, when $`E_T^2Q^2`$ the structure of the photon may be resolved.

Maxfield discussed the assumptions and problems encountered in the measurement of cross sections for low $`E_T`$ jet production by real and virtual photons and the extraction of the effective photon parton densities by H1. $`^\mathrm{?}`$ Measurements at low $`x_\gamma `$, the fractional momentum of the parton in the photon, require measuring low-$`E_T`$ jets. This demands a good experimental understanding and simulation of the pedestal transverse energy due to multiple interactions of the remnant partons of the photon and the proton in addition to the $`E_T`$ of the jets from the hard subprocess.

The photoproduction of dijets as a function of $`E_T`$ has been measured by H1 at high $`E_T`$ up to about 90 GeV and is in good agreement with NLO QCD calculations. $`^\mathrm{?}`$ A ZEUS measurement of the dijet cross section for lower $`E_T`$ as a function of $`x_\gamma `$, also presented by Wing, shows the data overshooting the NLO prediction by about 50% for $`x_\gamma <0.8`$, see Fig 3, which is considered to be due to inadequacies of the current parameterizations of the photon parton densities at scales of about $`E_T^2200`$ GeV2$`^\mathrm{?}`$ Improved photon parton densities should be expected from fits including this new HERA data.

Surrow showed results from OPAL at LEP on dijet production for the case of quasi-real photons from both colliding beams, i.e. both photon virtualities are almost zero. $`^\mathrm{?}`$ They are well described by NLO QCD calculations.

## 3 Measuring $`\alpha _s`$ using jets in deep-inelastic scattering

The extraction of the strong coupling parameter $`\alpha _s`$ with increasing precision is an important task of the HERA experiments. A new result from ZEUS, using measurements of the dijet cross section and of the dijet rate in bins of the photon virtuality $`Q^2`$ was presented by Hadig. $`^\mathrm{?}`$ In their analysis the parton density functions (pdfs) are taken from global fits of other measurements. For the first time, the error on $`\alpha _s`$ due to the uncertainty in the pdfs is estimated including their correlations $`^\mathrm{?}`$. From the dijet cross section they obtain

$`\alpha _s(M_Z)`$ $`=`$ $`0.1186\pm 0.0019(\mathrm{stat}.)_{0.0007}^{+0.0020}(\mathrm{exp}.)_{0.0033}^{+0.0035}(\mathrm{E}\mathrm{scale})`$
$`{}_{0.0038}{}^{+0.0048}(\mathrm{ren}.\mathrm{scale})\pm \mathbf{0.0031}(\mathrm{pdf})\pm 0.0005(\mathrm{hadr}.)`$
and from the dijet rate

$`\alpha _s(M_Z)`$ $`=`$ $`0.1166\pm 0.0019(\mathrm{stat}.)_{0.0005}^{+0.0023}(\mathrm{exp}.)_{0.0034}^{+0.0036}(\mathrm{E}\mathrm{scale})`$
$`{}_{0.0042}{}^{+0.0050}(\mathrm{ren}.\mathrm{scale})\pm \mathbf{0.0012}(\mathrm{pdf})\pm 0.0005(\mathrm{hadr}.).`$
It can be seen that the uncertainties due to the pdfs partially cancel when using the dijet rate. The extracted values have an experimental precision which allows to have an impact on the world average of $`\alpha _s`$. The major errors which need to be reduced in the future are the theory error due to the renormalization scale followed by the experimental error due to the hadronic energy scale of the calorimeter.

## 4 Azimuthal asymmetries in deep-inelastic scattering

A basic prediction of perturbative QCD, which is however difficult to verify experimentally, is the distribution of the azimuthal angle $`\mathrm{\Phi }`$ between the lepton scattering plane and the hadron production plane, defined by the exchanged virtual photon and an outgoing hadron. For parton production this azimuthal dependence has the form

$$d\sigma /d\mathrm{\Phi }=A+B\mathrm{cos}\mathrm{\Phi }+C\mathrm{cos}2\mathrm{\Phi },$$

resulting from the polarisation of the exchanged photon. Brook reported on such a measurement by ZEUS using single hard particles. $`^\mathrm{?}`$ The measured value for $`<\mathrm{cos}\mathrm{\Phi }>`$ is found to be negative, and the value for $`<\mathrm{cos}2\mathrm{\Phi }>`$, measured for the first time, is positive and increases with the transverse momentum of the particles. Both results agree, within large statistical errors, with QCD calculations.

## 5 Search for instantons

Instantons are a fundamental theoretical prediction of QCD. They exist in QCD and weak interactions and are associated with tunneling transitions between topologically non-degenerate vacua. Instanton-induced processes in QCD become calculable and perhaps observable in DIS due to the presence of a hard scale. Cross sections and characteristic event signatures can be calculated using the Monte Carlo generator QCDINS $`^\mathrm{?}`$ and were discussed by Schrempp, together with theoretical uncertainties on the event topology and the rates. $`^\mathrm{?}`$

Mikocki from the H1 collaboration presented first results on a dedicated search for instanton-induced processes in DIS$`^\mathrm{?}`$ They are based on some of the expected characteristics of the hadronic final state of such processes. Choosing simple cuts on three observables, requiring highest background reduction ($``$ 0.1%) while keeping an efficiency of $``$ 10% for instanton-induced processes, 549 events are observed in the data. For the background, i.e. standard DIS, $`435_{22}^{+36}`$ events are estimated using the leading order matrix elements and a parton shower model, which is part of RAPGAP $`^\mathrm{?}`$, and $`363_{26}^{+22}`$ using the colour dipole model, as implemented in ARIADNE $`^\mathrm{?}`$. The errors on the standard DIS background include only the experimental systematic errors and no model uncertainties. An excess in data over standard DIS events is seen which is qualitatively compatible with the expected instanton signal. However, the size and shape of this excess is at the level of the discrepancy between the standard QCD models themselves.

## 6 Status of higher-order calculations

QCD predictions in leading order (LO) of perturbation theory suffer generally from large higher-order corrections and scale uncertainties. They depend on the renormalization scale in the strong coupling constant $`\alpha _s`$ and on the factorization scales in the parton densities and fragmentation functions. At next-to-leading order (NLO), the renormalization and factorization scales appear also explicitly in the hard cross section, thus reducing the scale dependence considerably and making the theoretical predictions more reliable.

NLO QCD calculations for HERA physics have been performed for jet production in DIS, real and virtual photoproduction, heavy flavour production in DIS and photoproduction, prompt photon production in DIS and photoproduction, and inclusive hadron photoproduction (for references see e.g. $`^{\mathrm{?},\mathrm{?}}`$).

The improved performance of HERA, the upcoming luminosity upgrade, and the better understanding of the H1 and ZEUS detectors contribute to reduced statistical and systematic errors on experimental measurements, with the exception of the still substantial energy scale uncertainty. Consequently, comparisons between experimental data and NLO QCD predictions are more and more dominated by theoretical scale uncertainties, limiting the discriminating power of the data for $`\alpha _s`$ and parton density determinations. For many observables it appears thus mandatory to improve the theoretical predictions beyond NLO. The first option consists in calculations in next-to-next-to-leading order (NNLO) of QCD.

Glover reported that the last three years have seen significant progress in this technically difficult field: Double unresolved tree amplitudes and single unresolved one-loop amplitudes have been calculated in the infrared limit, and many massless two-loop integrals have been evaluated, the most difficult ones being the double-box diagrams. $`^\mathrm{?}`$

The remaining challenges are to integrate out the infrared behavior of the unresolved cross sections and to analytically evaluate the two-loop matrix elements. It might still take several years before the first NNLO Monte Carlo program becomes available.

## 7 New results from resummation

An alternate route is to resum QCD corrections from gluon radiation in the soft (small transverse momentum $`q_T`$) and threshold (large $`x_T=2q_T/\sqrt{s}`$) regions of phase space. Resummation techniques rely on the observation that the perturbative series exponentiates. They have been applied to $`e^+e^{}`$ annihilation, Drell-Yan production, and very recently also to DIS, where more work is needed. For photoproduction resummed calculations are not available.

Data for $`q_T`$ distributions in semi-inclusive DIS do not exist. However, the rapidity distribution of the transverse energy flow has been measured by H1 $`^\mathrm{?}`$, and the rapidity $`\eta `$ can be related to $`q_T`$ in the Breit frame if the photon virtuality $`Q`$ is fixed ($`q_T=Qe^\eta `$). Nadolsky showed that the data compare favourably to the resummed predictions. $`^\mathrm{?}`$

Resummed predictions for DIS event shape variables like thrust have been calculated by Dasgupta, but have not yet been confronted with data. $`^\mathrm{?}`$ When compared to NLO thrust distributions, they show a constant difference with respect to DISASTER++ $`^\mathrm{?}`$ as expected, but a thrust-dependent difference with respect to DISENT $`^\mathrm{?}`$.

## 8 Prompt photon production

A physical process where resummation plays a very important role is prompt photon production. Prompt photons could help to constrain the gluon density in the proton and the photon, since initial gluons generally dominate the production cross sections by 80% or more. Unfortunately, this possibility has been put on hold since the transverse momentum distributions of prompt photons in fixed target (e.g. E706 $`^\mathrm{?}`$) and collider (e.g. UA2 $`^\mathrm{?}`$ and CDF $`^\mathrm{?}`$) experiments differ widely from the perturbative NLO expectations. aaaThis is not the case for off-shell photons, which might therefore be better suited to determine the unpolarized and polarized gluon densities in the proton.

While threshold resummation reduces the scale dependence of the theoretical prediction and improves the description of the E706 data at large $`q_T`$, it does not increase the cross section at small $`q_T`$.

It has been speculated that intrinsic $`q_T`$ effects could be responsible for the discrepancy meaning that the partons in the colliding hadrons are not exactly collinear. Phenomenologically this can be motivated by re-expressing the delta-function $`\delta (q_T)`$ as a Gaussian exponential function or using parton densities which are not integrated over $`q_T`$. These considerations must be put on firm theoretical ground: Recoil effects have to be taken into account, and the modified parton distributions have to be related to the usual ones which are used in other scattering processes.

A possible solution, reported on by Vogelsang, consists in the form of a $`q_T`$-profile function, which also permits to resum the low and high $`q_T`$ regions simultaneously. $`^\mathrm{?}`$ Although the numerical result still lies somewhat below the E706 data at low $`q_T`$, the discrepancy is significantly reduced, see Fig. 4.

A drawback of this (as of any resummed) approach is that it is necessary to match the resummed and perturbative predictions at a fixed value of $`q_T`$ which is a priori undetermined.

Prompt photon production has also been measured by ZEUS at HERA in photon-proton scattering and was reported on by Lee. $`^\mathrm{?}`$ The measurements cover the range in transverse momentum $`p_T^\gamma `$ from 5 to 15 GeV and in pseudorapidity $`\eta ^\gamma `$ from -0.7 to 0.9. The $`p_T^\gamma `$ distribution agrees very well with NLO QCD predictions, shedding doubt on the universality of intrinsic $`q_T`$ effects. However, the experimental errors are still dominated by statistics. An update of the analysis is in progress. It is also important to keep in mind that the colliding photon at HERA has a much smaller transverse size than a proton (or even a nucleon) in hadron collisions, which leads to smaller intrinsic $`q_T`$ effects at least on the photon side. The NLO calculations describe the $`\eta ^\gamma `$ distribution only well in forward (proton) direction, but are below the data in the backward region. This would seem to indicate an inadequacy of the current photon parton densities.

## 9 Non-perturbative hadronization corrections

Both NNLO and resummed perturbative calculations are furthermore subject to non-perturbative corrections. These arise, e.g., when partons hadronize into jets, and they can become very important in certain regions of phase space like the backward region of jet photoproduction.

A possible remedy is the combination of perturbative predictions with parton showers and hadronization from existing Monte Carlo models. These models require a well-defined colour flow for the underlying hard scattering matrix element which poses a problem for interference terms. In a first approximation they can be split up and added to the squared matrix elements. Furthermore, beyond LO one has to be careful not to doublecount emission of real partons in the higher order matrix element and in the parton shower. The two can be separated by cutting on the invariant mass of the parton pair.

## 10 Fragmentation functions

The inclusive cross section for the production of single charged hadrons is expressed in terms of universal fragmentation functions (FF) which incorporate the long-distance, non-perturbative physics of the hadronization process. The scale dependence of the fragmentation functions is governed by an evolution equation, similar to the evolution equation for parton densities. The initial condition for the evolution equation, however, is not calculable in perturbation theory and must at present be taken from experiment.

Two new analyses to extract fragmentation functions for charged hadrons from $`e^+e^{}`$ annihilation data have been presented by Pötter and Kretzer. $`^{\mathrm{?},\mathrm{?}}`$ Both analyses have been performed at leading and next-to-leading order in QCD. The fitted fragmentation functions describe all data sets within their respective errors. A comparison between the two independent FF extractions reveals good agreement when the sum over all flavours is taken. However, significant differences are found for the individual flavour fragmentation functions into $`\pi ^\pm `$ and $`K^\pm `$, which are not well constrained by the data. $`^\mathrm{?}`$

Using experimental measurements of inclusive single hadron production at different energies one can extract the strong coupling $`\alpha _s`$ from the scaling violations in the fragmentation functions. A new $`\alpha _s`$ determination with a competitive error and a value consistent with the world average has been presented at this workshop. $`^\mathrm{?}`$

## 11 Inclusive particle production and MLLA

Safonov presented results on jet fragmentation, i.e. inclusive momentum distribution and multiplicities in jets, from CDF. $`^\mathrm{?}`$ He demonstrated that calculations in the modified-leading approximation (MLLA) and assuming local parton hadron duality can provide a simple and compact description of the data. Milstead reported on inclusive particle production in the Breit frame in DIS as studied by H1 and ZEUS. $`^\mathrm{?}`$ Here a consistent description of various moments of the fragmentation function by MLLA has not yet been achieved.

## 12 Event shapes and power corrections

The first new ZEUS measurements of event shapes were reported by Wollmer along with those from H1 by Rabbertz. $`^\mathrm{?}`$ These measurements allow for an extraction of $`\alpha _s`$ and another parameter $`\alpha _0`$ that characterizes the behavior of the strong coupling in both the high energy and low energy (non-perturbative) regimes. Several event shapes are considered, but there is a rather wide spread in results. Biebel reviewed the measurements from LEP $`^\mathrm{?}`$ and Dasgupta the theoretical situation $`^\mathrm{?}`$ of power corrections. The conclusion so far is that there is still a need for deeper understanding of the experimental and theoretical uncertainties associated with these measurements.

## 13 Photoproduction of charm and beauty

The production of heavy flavours provides new opportunities to study the dynamics of perturbative QCD and to extract information on the proton and photon structure. HERA has provided a wide spectrum of $`D^{}`$ measurements and first results on $`D_s`$ mesons by ZEUS were reported by Gladilin. $`^\mathrm{?}`$ It is not guaranteed a priori that charm cross sections can be reliably predicted in perturbation theory. The NLO calculations are indeed plagued by large uncertainties and they tend to underestimate the $`D^{}`$ cross section, in particular in the forward direction and at low $`x_\gamma `$ where resolved photon processes contribute. Potentially large next-to-next-to-leading order corrections and higher-twist contributions may have to be included to improve the theoretical predictions. Other possible explanations for the discrepancies include non-perturbative string effects between the proton remnant and the charm quark or an enhancement of the low-$`x`$ gluon component of the photon, which is currently not well constrained experimentally. A similar picture has emerged in $`\gamma \gamma `$ collisions at LEP. As shown by Andreev, the two-photon $`D^{}`$ cross section is slightly underestimated by NLO theory in the experimentally visible region, in particular at low $`p_t`$ where resolved photon processes are prominent. $`^\mathrm{?}`$

Due to the larger $`b`$ quark mass, theoretical predictions should be under much better control for the beauty cross section. Still, it is well known that Tevatron data for the $`b`$ transverse momentum distribution are systematically larger than NLO QCD predictions. This trend is supported by recent HERA and LEP measurements which show that the beauty photoproduction cross section is significantly higher than the NLO theory. $`^{\mathrm{?},\mathrm{?}}`$ Kuhr presented the H1 and ZEUS results displayed in Fig. 5.

Also the beauty cross section in $`\gamma \gamma `$ collisions at LEP as measured by L3 $`^\mathrm{?}`$

$$\sigma _{tot}(eeeeb\overline{b}X)=(9.9\pm 2.9\pm 3.8)\text{pb}\text{[L3]}$$

is higher than the NLO QCD calculation of $`3`$ pb. The discrepancy between NLO theory predictions and experiment in beauty production should be taken seriously and possible theoretical explanations have been discussed by Ellis in the final plenary session. $`^\mathrm{?}`$ The HERA luminosity upgrade and alternative experimental techniques using microvertex-detection will provide important additional information in the future. $`^\mathrm{?}`$

## 14 Charm production in deep-inelastic scattering

The charm contribution to the total structure function $`F_2`$ at small $`x`$ at HERA is sizeable, up to $`25\%`$. A proper description of the charm contribution to DIS is thus required for a global analysis of structure function data and a precise extraction of the parton densities in the proton.

Considerable theoretical effort has been devoted to including heavy quark effects in DIS. As reviewed by Smith and Tung, different so-called variable flavour number prescriptions have been defined in the literature. $`^{\mathrm{?},\mathrm{?}}`$ These prescriptions incorporate the correct heavy quark threshold behavior and sum terms $`\alpha _s\mathrm{log}Q^2/M_c^2`$ when the hard scale $`Q`$ of the reaction is significantly larger than the heavy-quark mass $`M_c`$. The structure functions calculated using different prescriptions generally differ at finite order in perturbation theory, since the treatment of the matching conditions and the extent to which finite mass effects are retained in the coefficient functions is ambiguous.

Experimentally, charm production in DIS at HERA has not only been studied for inclusive structure functions but differential cross sections for $`D^{}`$ production are available as well. $`^\mathrm{?}`$ The predictions for differential distributions in variable flavour number prescriptions at $`𝒪(\alpha _s)`$ are spoiled by large scale dependences, with the exception of the $`Q^2`$ and $`p_t(D^{})`$ distributions where the agreement with data is good. $`^\mathrm{?}`$ The scale dependence is significantly reduced when higher-order corrections are included. The comparison of HERA data with $`𝒪(\alpha _s^2)`$ calculations $`^\mathrm{?}`$ shows good overall agreement, with a possible slight excess of data in the forward region. Hadronization effects may well be responsible for the deviation, similar to $`D^{}`$ photoproduction. A first measurement of the charm contribution to the photon structure function has been presented at this workshop and satisfactory agreement, within the still sizeable statistical uncertainty, with theory was found. $`^\mathrm{?}`$

## 15 Quarkonium production

Exciting phenomenological developments in quarkonium physics followed from the application of non-relativistic QCD (NRQCD), an effective field theory that includes the so-called colour-octet mechanisms. Although gluon fragmentation into colour-octet charm quark pairs appears as the most plausible explanation of the large direct $`\psi `$ production cross section observed at the Tevatron, the applicability of NRQCD factorization to charmonium production is still not established quantitatively.

Inclusive charmonium production at HERA offers unique possibilities to assess the importance of different quarkonium production mechanisms. Krüger showed that no conclusive evidence for colour-octet processes has been observed in $`J/\psi `$ photoproduction or $`J/\psi `$ production in DIS so far. $`^\mathrm{?}`$ The measurements of $`J/\psi `$ production through resolved photon processes, for which first results have been reported at this workshop $`^\mathrm{?}`$, could be important to clarify the issue once more statistics has been accumulated.

The single most crucial test of the NRQCD approach is the analysis of direct $`J/\psi `$ and $`\psi ^{}`$ polarization at large $`p_t`$ at the Tevatron. NRQCD predicts a substantial fraction of transverse polarization at $`p_t`$ above $`10`$ GeV, as discussed in Lee’s review. $`^\mathrm{?}`$ As shown in Fig. 6 this prediction is not supported by recent Tevatron data. The absence of transverse polarization in $`J/\psi `$ and $`\psi ^{}`$ hadroproduction at large $`p_t`$, if confirmed by the higher statistics data expected at Run II at the Tevatron, would represent a serious problem for the application of NRQCD to charmonium production.

## 16 Discussion session

After 33 experimental and 12 theoretical presentations we ended our working group meetings with a highly attended discussion session on future directions in QCD at HERA. Three ‘provocateurs’, M. Wing, J. Repond, and K. Ellis started off the discussion.

Wing presented his views on the interesting directions in photoproduction. $`^\mathrm{?}`$ Ellis expressed his discomfort with the experimental results of too high cross sections for $`b\overline{b}`$-production compared to current NLO QCD calculations, and he discussed possible theoretical solutions. Repond pointed out that while most of us believe QCD to be the theory of strong interactions, we do lack precision, i.e. $`\alpha _s`$ is only known to about 3%, we want to know where perturbative QCD is good enough, where are higher order or resummation calculations needed and how to do them correctly. To pursue this experimentally, precision data and wide coverage of the available phase space are needed. The luminosity upgrade of HERA will allow to reduce statistical and with much effort by the experimentalists also systematic errors. Also needed are of course precision theory which already now is often trailing the experiments. Repond convincingly pleaded for an increase in the manpower in theory to address, together with the experimentalists, the outstanding QCD-issues at HERA.

We had an interesting discussion and good participation on the issues mentioned above and some others brought up by the audience until we ran out of time. Clearly we did not have enough time, probably we had too many topics and not good enough structuring. The interest in such discussions however could be seen from the lively debate we had and the wishes for more by many people afterwards.

## 17 Executive conclusions

* For many results, on jet cross sections, the determination of $`\alpha _s`$, on charm cross sections, the major uncertainty is theoretical due to the renormalization scale, followed by the experimental error due to the hadronic energy scale of the calorimeter. Both need to be improved in the future.
* Using dijet rates instead of the dijet cross section, the error on $`\alpha _s`$ due to the uncertainties in the parton density functions could be halved. In addition, this error can now be determined, since parton densities with correlated errors have become available.
* Results on dijet production suggest that the parton densities of the photon should be revised and improved by including HERA and LEP data in a global fit.
* A first dedicated search for instantons in deep-inelastic scattering was reported by H1.
* Several important pieces of next-to-next-to-leading order calculations have been computed. The challenge is now to calculate the remainder and to assemble the pieces.
* Resummed calculations for semi-inclusive hadron production and thrust distributions at HERA are now available. Predictions for prompt photon and jet production are still missing.
* Next-to-leading order QCD calculations for charm and beauty photoproduction underestimate experimental data from HERA and LEP. The discrepancy is particularly significant in the case of beauty production where perturbative calculations should be on safe grounds.
* In contrast, the comparison of HERA and LEP data on charm production in DIS shows good overall agreement with $`𝒪(\alpha _s^2)`$ calculations. A proper description of the charm contribution to DIS is important for a precise extraction of the parton densities in the proton.
* HERA measurements so far show no conclusive evidence for colour-octet contributions to charmonium production in $`\gamma p`$ and deep-inelastic reactions. The absence of transverse polarization in $`J/\psi `$ and $`\psi ^{}`$ hadroproduction at large $`p_t`$ represents a serious problem for the application of non-relativistic QCD to charmonium production.
* Finally, one important echo from the discussion session of the last day was that given the high experimental precision which has been reached, more theoretical support is clearly needed.

## Acknowledgments

M. Klasen is supported by Deutsche Forschungsgemeinschaft under Contract KL 1266/1-1. M. Klasen and M. Krämer are supported in part by the European Commission under Contract ERBFMRXCT980194.

## References