text:
# The scaling region of the lattice 𝑂(𝑁) sigma model at finite temperature

## 1 Introduction

The behavior of symmetries at finite temperature is one of the most outstanding and relevant problems in current areas of particle physics such as cosmology, relativistic heavy-ion collisions and the quark-gluon plasma. It was shown some time ago that spontaneously broken symmetries in quantum field theories get restored at sufficiently high temperature . In the high temperature limit, however, perturbation theory becomes unreliable when infrared divergences are not properly summed. Various groups developed temperature dependent renormalization group techniques to study the critical properties of the $`O(N)`$ scalar field theories (for a few examples see  and references therein).

In recent years considerable work has been done on the physics of the finite temperature chiral phase transition in QCD and related models. A compelling idea based on universality and dimensional reduction first put forward in , is that in QCD with two massless quarks the physics near the transition can be described by the three-dimensional $`O(4)`$-symmetric sigma model. Although, there is little disagreement that the chiral phase transition in QCD is second order, no quantitative work or simulations have been done that decisively determine its universality class, because simulations are performed away from the continuum limit on small volumes and show substantial discretization and finite size effects. The knowledge of the scaling relations are crucial for a successful extrapolation to small mass and large volume near $`T_c`$.

The reasoning behind the dimensional reduction scenario is based on the dominance of the light degrees of freedom near $`T_c`$. In the imaginary time formalism of finite temperature field theory a given system is defined on a manifold $`S^1\times R^d`$ with the inverse temperature $`\beta T^1`$ being the circumference of $`S^1`$. The boundary condition in the time direction is periodic for bosons and anti-periodic for fermions. At temperatures larger than the physical scales of the system, the nonzero Matsubara modes acquire a heavy mass and non-static configurations are strongly suppressed in the Boltzmann sum. If the original four-dimensional theory has fermions, then because of the anti-periodicity in the temporal direction, all fermionic modes become massive and can be integrated over. This is not to say that non-static modes have no effect. Rather they generate counter-terms to the three-dimensional theory. The physical picture of the dimensional reduction scenario can be described nicely by using a slightly different language presented in . As one approaches the critical point, the fluctuations in space become correlated over a distance $`\xi `$ which diverges at the critical point. In the critical region, not only do the fluctuations extend over large spatial scales, they also become very slow and there is a divergent relaxation time $`\tau `$; the effect is known as critical slowing down. The appropriate quantum of energy is then $`\mathrm{}\omega _c`$ with $`\omega _c1/\tau `$. Hence if the transition is at a finite temperature $`T_c`$, $`\mathrm{}\omega _ck_BT`$ near the critical point, which means that long wavelength fluctuations are classical thermal fluctuations and quantum mechanics is irrelevant. Equivalently, since the time interval $`\beta `$ is very short compared to the relaxation time $`\tau `$, typical time histories will consist of a single configuration which is the same at each time slice. The dynamics therefore drops out of the problem and a Boltzmann weight $`\mathrm{exp}(\beta H_{\mathrm{classical}})`$ is recovered from the path integral.

In this work we study numerically the finite temperature transition both in the $`(3+1)d`$ and $`(2+1)d`$ $`\varphi ^4`$ model. This is equivalent to studying the phenomenon of dimensional crossover in a $`d`$-dimensional layer $`O(N)`$ spin model of infinite extent in $`d1`$ dimensions and finite thickness in the remaining dimension. In other words temperature effects become predominant as a consequence of a change in dimension. The formal derivations of the dimensional reduction scenario require that the scalar mass is small compared to the temperature. We check whether this condition is necessary and we study the effects of the infrared fluctuations on the width of the critical region under different circumstances. We performed simulations of the model with different $`O(N)`$ (for $`N=1,2`$ and $`3`$) symmetries in $`(3+1)d`$ and study the dependence of the width of the $`3d`$ scaling region on $`N`$. In the case of the $`(3+1)d`$ $`O(2)`$ model we compare the widths of the scaling region in the low temperature broken phase and the high temperature symmetric phase. We also study the scaling region of the $`(2+1)d`$ $`O(1)Z_2`$ model. In addition we are interested to see how finite size effects can distort the results in such a simple model that can be simulated on lattices with relatively large temporal and spatial sizes. We extract critical exponents by measuring the order parameter, the susceptibility, and the correlation length.

## 2 The model

In the continuum limit the Lagrangian of the model reads:

$$=\frac{1}{2}(\stackrel{}{\varphi })^2+\frac{m_0^2}{2}\stackrel{}{\varphi }^2+\frac{g}{4!}(\stackrel{}{\varphi }^2)^2.$$
(1)
In order to study the problem of symmetry restoration in a non-perturbative way, we regard this model as an effective low energy theory by introducing a lattice regularization with a lattice constant $`a=1/\mathrm{\Lambda }`$, where $`\mathrm{\Lambda }`$ is the ultraviolet cut-off. Summing the bare Lagrangian over all lattice sites and transforming the continuum parameters into suitably defined lattice parameters according to a standard convention

$$a\stackrel{}{\varphi }(x)=\kappa ^{1/2}\stackrel{}{\mathrm{\Phi }}_x,a^2m_0^2=\frac{24\lambda }{\kappa }8,g=\frac{24\lambda }{\kappa ^2},$$
(2)
we are led to the following action:

$$S[\stackrel{}{\mathrm{\Phi }}]=\underset{x}{}\left[\kappa \underset{\mu }{}\mathrm{\Phi }_x^\alpha \mathrm{\Phi }_{x+\mu }^\alpha +\mathrm{\Phi }_x^\alpha \mathrm{\Phi }_x^\alpha +\lambda (\mathrm{\Phi }_x^\alpha \mathrm{\Phi }_x^\alpha 1)^2\right],$$
(3)
where the index $`\mu `$ runs over the nearest neighbors, ($`\mu =1,\mathrm{},d`$), and the index $`\alpha `$ runs over the $`N`$ components of the field $`\stackrel{}{\mathrm{\Phi }}`$. For convenience, in our numerical work we concentrated on the limit of an infinitely strong bare coupling where the field $`\stackrel{}{\mathrm{\Phi }}`$ is constrained to unit length: $`\stackrel{}{\mathrm{\Phi }}\stackrel{}{\mathrm{\Phi }}=1`$. In this case the action $`S`$ reduces to the form:

$$S[\mathrm{\Phi }]=\kappa \underset{x,\mu }{}\mathrm{\Phi }_x^\alpha \mathrm{\Phi }_{x+\mu }^\alpha ,$$
(4)
which is equivalent to the energy of a configuration of the classical ferromagnetic $`O(N)`$ spin model.

Finite physical temperature on the lattice means a finite size of the temporal extent $`L_t`$. The physical temperature is given by $`T=1/(L_ta)`$. The temperature of the system may be altered, by varying either the coupling $`\kappa `$, which amounts to varying the lattice spacing or by varying $`L_t`$. To minimize the effects of lattice artifacts the condition $`\mathrm{\Lambda }T`$ must be satisfied which requires a lattice temporal extent $`L_t1`$. This condition not being met implies that the sum over the Matsubara modes is truncated at a small value. It is well-known that the $`(3+1)d`$ renormalized $`\varphi ^4`$ theory is a free theory. However, the regularized theory may well describe the correct physics of a broad spectrum of processes below the cut-off scale.

As discussed in  the regularized model can be interpreted in terms of physical terms by moving along the curves of constant physics (CCPs) which are characterized by $`g_R=\mathrm{const}.`$ Therefore, when approaching the $`T=0`$ second order critical line the mass of the sigma meson $`M_\sigma `$ decreases and in order to keep $`T/M_\sigma `$ constant one has to increase $`L_t`$. Near the $`T=0`$ critical line the ratio $`T/M_\sigma `$ should not change any more along a CCP. It was shown within a specific regularization scheme to leading order in $`1/N`$  that $`T_c/M_\sigma =\sqrt{36/(Ng_R)}`$ and by solving the gap equation on the lattice it was shown  that the corrections due to a finite $`L_t`$ are negligible for $`L_t4`$. In order to suppress possible discretization effects arising from $`1/N`$ corrections we performed our simulations with $`L_t=12`$. The $`(2+1)d`$ $`\varphi ^4`$ theory has an interacting continuum limit and it has been shown that the non-linear sigma model is renormalizable in the $`1/N`$ expansion . The phenomenon of symmetry restoration at finite temperature in the three-dimensional theory was investigated by means of lattice simulations and finite size scaling techniques by Bimonte et al. . The authors of  located the intersections of critical lines corresponding to different values of $`L_t`$ with various CCPs lying in the ordered phase of the $`T=0`$ theory and showed that for $`L_t3`$ the critical temperatures approach the finite physical temperature of the continuum theory.

Furthermore, the condition $`L_s\xi `$ must be satisfied in order to obtain results that approximate the thermodynamic limit. For this reason at certain values of coupling $`\kappa `$ we performed simulations on lattices with different spatial extent in order to detect and control finite size effects.

In our simulations we used the single-cluster algorithm , because it reduces dramatically the correlations between successive configurations near a critical point. The cluster update is as follows: First choose randomly a reflection axis in the plane, denote the component of the spin $`\stackrel{}{\mathrm{\Phi }}_i`$ that is parallel to this reflection axis by $`\mathrm{\Phi }_i^{||}`$ and that which is orthogonal by $`\mathrm{\Phi }_i^{}`$. Then choose randomly a site $`i`$ of the lattice as a starting point for the cluster $`𝒞`$. Visit all neighbor sites $`j`$ of $`i`$ and allow them to join the cluster with probability

$$p(i,j)=1\mathrm{exp}[\kappa (\mathrm{\Phi }_i^{}\mathrm{\Phi }_j^{}+|\mathrm{\Phi }_i^{}\mathrm{\Phi }_j^{}|)].$$
(5)
After this is done, visit all neighbors of the new sites in the cluster and add them to the cluster with probability $`p(i,j)`$. Iterate this step until no new sites enter the cluster. Finally flip the signs of all $`\mathrm{\Phi }_i^{}`$ contained in the cluster.

## 3 The observables

In this section we present the various observables measured in the broken and symmetric phases of the model in our lattice simulations. The magnetization $`M^\alpha `$ of a given configuration is defined as:

$$M^\alpha \frac{1}{V}\underset{x}{}\mathrm{\Phi }_x^\alpha ,$$
(6)
where $`VL_s^3L_t`$ is the space-time volume. However, in a finite volume and without the benefit of an explicit symmetry breaking interaction with an external magnetic field, the direction of symmetry breaking changes over the course of the run so that $`M^\alpha `$ averages to zero over the ensemble. It is in this way that the absence of spontaneous symmetry breaking on a finite lattice is enforced. Another option is to introduce a field operator

$$\mathrm{\Phi }_{\sigma ,x}\frac{\mathrm{\Phi }_x^\alpha M^\alpha }{|M|},$$
(7)
which is a projection of $`\mathrm{\Phi }_x^\alpha `$ to the direction of $`M^\alpha `$ separately for each configuration . The effective order parameter

$$\mathrm{\Sigma }\frac{1}{V}\underset{x}{}\mathrm{\Phi }_{\sigma ,x}=|M|,$$
(8)
differs from the true order parameter extrapolated to the zero magnetic field limit by a factor $`V^{\frac{1}{2}}`$ . One can easily show that the susceptibility of the order parameter can be expressed in terms of moments of the magnetization as follows:

$$\chi =V(|M|^2|M|^2).$$
(9)
In the symmetric phase one gets $`\chi =V|M|^2`$. Cluster algorithms enable us to reduce the variance of the expectation values in the symmetric phase by using improved estimators . The improved estimator for the susceptibility is given by:

$$\chi _{\mathrm{imp}}=\frac{N}{|C|}\left(\underset{i𝒞}{}\mathrm{\Phi }_i^{}\right)^2,$$
(10)
where $`|C|`$ denotes the number of spins in the cluster $`𝒞`$.

We determined the temporal exponential correlation length $`\xi _{\mathrm{exp}}^t`$ of the sigma meson from the asymptotic exponential decay of the zero momentum connected two-point correlation function constructed from the field $`\mathrm{\Phi }_{\sigma ,x}`$, i.e.,

$`O(t)`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{x}}{}}\mathrm{\Phi }_{\sigma ,x}`$ (11)
$`G(t)`$ $``$ $`O_0O_t\mathrm{\Sigma }^2\stackrel{t1}{=}A\left[\mathrm{exp}\left({\displaystyle \frac{t}{\xi _{\mathrm{exp}}}}\right)+\mathrm{exp}\left({\displaystyle \frac{(L_t1)}{\xi _{\mathrm{exp}}}}\right)\right].`$ (12)
Similarly, we measured the spatial correlation length $`\xi _{\mathrm{exp}}^s`$ from the exponential decay of $`G(x)`$ in the spatial direction. The asymmetry of the lattice at nonzero temperature implies that $`\xi _{\mathrm{exp}}^s\xi _{\mathrm{exp}}^t`$, because temperature effects on the mesonic self-energy become direction dependent. The breaking of Lorentz invariance implies that the spectral functions and the dispersion relations are modified in a non-trivial way by the thermal statistical ensemble. However, both $`\xi _{\mathrm{exp}}^s`$ and $`\xi _{\mathrm{exp}}^t`$ are expected to diverge as $`TT_c`$ according to the scaling law $`(T_cT)^\nu `$, where $`\nu `$ is a critical exponent of the corresponding $`O(N)`$ $`3d`$ spin model. One can construct correlation functions from various other operators $`O(t)`$ that couple to the scalar particle $`\sigma `$, but according to  the operators we defined in eq.(11) give a better signal than other options. In the symmetric phase the two-point correlation function is $`G(t)=(_i\stackrel{}{\mathrm{\Phi }}(\stackrel{}{x}_i,t))_j(\stackrel{}{\mathrm{\Phi }}(\stackrel{}{x}_j,0))`$ and its improved estimator is given by

$$G(t)_{\mathrm{imp}}=\frac{N}{|𝒞|}\delta _{ij}(𝒞)\mathrm{\Phi }_i^{}\mathrm{\Phi }_j^{},$$
(13)
where $`\delta _{ij}(𝒞)=1`$ if $`i`$ and $`j`$ belong to the same cluster otherwise $`\delta _{ij}(𝒞)=0`$ .

We also measured the so-called second moment correlation length $`\xi _{2\mathrm{n}\mathrm{d}}`$, defined by

$$\xi _{2\mathrm{n}\mathrm{d}}^2\frac{1}{2d}\frac{_x|x|^2G(x)}{_xG(x)}=\left(\overline{G}^1(k^2)\frac{d}{dk^2}\overline{G}(k^2)\right)_{k^2=0},$$
(14)
where $`\overline{G}(k)_j\mathrm{exp}(ikx_j)\stackrel{}{\mathrm{\Phi }}_0\stackrel{}{\mathrm{\Phi }}_j`$ is the Fourier transform of the two-point correlation function. On the lattice $`\xi _{2\mathrm{n}\mathrm{d}}`$ is expressed as

$$\xi _{2\mathrm{n}\mathrm{d}}=\left(\frac{(\chi /F)1}{4\mathrm{sin}^2(\pi /L)}\right)^{1/2},$$
(15)
where $`L=L_t`$ or $`L_s`$, $`F\overline{G}(k)|_{|k|=2\pi /L}`$. The improved $`\xi _{2\mathrm{n}\mathrm{d}}`$ requires $`\overline{G}(k)_{\mathrm{imp}}`$ given by

$$\overline{G}(k)_{\mathrm{imp}}=\frac{N}{|𝒞|}\left[\left(\underset{i𝒞}{}\mathrm{\Phi }_i^{}\mathrm{cos}(kx_i)\right)^2+\left(\underset{i𝒞}{}\mathrm{\Phi }_i^{}\mathrm{sin}(kx_i)\right)^2\right].$$
(16)
This estimator of the correlation length is very popular, because its numerical evaluation is simpler than that of the exponential correlation length which requires an identification of an asymptotic exponential decay. However, the two definitions of the correlation length do not coincide, since in $`\xi _{\mathrm{exp}}`$ only the ground state enters, while in the case of $`\xi _{2\mathrm{n}\mathrm{d}}`$ according to eq.(14) a mixture of the full spectrum is taken into account. As mentioned in  the difference from one of the ratio $`\xi _{\mathrm{exp}}/\xi _{2\mathrm{n}\mathrm{d}}`$ gives an idea of the density of the lowest states of the spectrum. It is well-known however, that in the high temperature phase the asymptotic exponential decay behavior of $`G(t)`$ is reached at much smaller values of $`t`$ than in the broken phase (see e.g. ). This behavior of $`G(t)`$ implies that the difference between $`\xi _{\mathrm{exp}}`$ and $`\xi _{2\mathrm{n}\mathrm{d}}`$ in the symmetric phase is smaller than in the broken phase. In our simulations we measured $`\xi _{2\mathrm{n}\mathrm{d}}`$ only in the symmetric phase and we observed that the difference between $`\xi _{\mathrm{exp}}`$ and $`\xi _{2\mathrm{n}\mathrm{d}}`$ vanishes within error.

## 4 Results

In this section we present results from numerical simulations of the four-dimensional $`O(N)`$ ($`N=1,2`$ and $`3`$) non-linear sigma model at $`T0`$. We also present results from simulations preformed in the low temperature phase of the three-dimensional Ising model. By measuring the correlation length $`\xi `$ at the crossover into the $`3d`$ scaling region we study the dependence of the width of the critical region on the symmetry group $`O(N)`$ and on whether the symmetry is manifest ($`T>T_c`$) or broken ($`T<T_c`$). We define the crossover coupling $`\kappa _{\mathrm{cross}}`$ as the value of $`\kappa `$ that if it is included in the fitting window, the extracted exponent $`\beta _{\mathrm{mag}}`$ does not agree within statistical error with the exponent of the $`3d`$ model $`\beta _{\mathrm{mag}}^{3\mathrm{d}}`$. As we show in the next paragraphs, the change in the value of the exponent is accompanied by a substantial increase in the $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}.`$ of the fit, implying that the quality of the fit drops rapidly once we enter the crossover region. We managed to identify the value of $`\kappa _{\mathrm{cross}}`$ with relatively good accuracy, because as we will see in the following paragraphs the value of $`\kappa `$ was varied by small increments at the crossover into the $`3d`$ scaling region. Throughout this work we used a large lattice temporal extent $`L_t=12`$ in order to minimize discretization effects.

First, we present the results for the $`(3+1)d`$ Ising model. For certain values of $`\kappa `$ we monitored finite size effects on the order parameter by simulating the model on lattices with two different values of spatial extent $`L_s`$. For the final analysis we chose the data sets for which the systematic errors on $`\mathrm{\Sigma }`$ are less than $`0.5\%`$. More precisely for the range of couplings $`0.15040\kappa 0.15100`$ we chose $`L_s=48`$, for $`0.15010\kappa 0.15030`$ $`L_s=72`$ and for $`0.14993\kappa 0.150075`$ $`L_s=96`$. For each value of $`\kappa `$ we typically generated $`10^6`$ sweeps of the algorithm. The magnetization data was fitted to the form

$$\mathrm{\Sigma }=A(\kappa \kappa _c)^{\beta _{\mathrm{mag}}}.$$
(17)
The results for different ranges of $`\kappa `$ are presented in Table 1. The simple power law fit is very good for data close to the transition, i.e. for $`\kappa 0.15005`$. For these fits the exponent $`\beta _{\mathrm{mag}}`$ is consistent within error with a recent estimate of the critical index of the $`3d`$ Ising model $`\beta _{\mathrm{mag}}^{3\mathrm{d}}=0.3269(6)`$ . As shown in Table 1, the value of $`\beta _{\mathrm{mag}}`$ deviates from the $`3d`$ values and the quality of the fit gets worse as we include more data points at larger values of $`\kappa `$. Therefore, we conclude that the crossover to the $`3d`$ universality class occurs at $`\kappa _{\mathrm{cross}}0.150075`$. The data and the fitting curve to the data with $`\kappa 0.150025`$ are shown in Fig. 1. The statistical errors are smaller than the size of the symbols and for this reason we don’t plot them in the figure. The data for the susceptibility and the correlation length for these values of the coupling were very noisy and did not allow us to measure the exponents $`\gamma `$ and $`\nu `$.

We measured the values of the spatial and temporal correlation lengths at the crossover into the $`3d`$ Ising scaling region by performing a simulation at $`\kappa =0.150075`$ on a $`12\times 48^3`$ lattice. For this measurement we chose a lattice with smaller spatial size $`L_s`$, because the correlation function was very noisy and we had to generate $`15`$ million configurations in order to measure the correlation lengths $`\xi ^t`$ and $`\xi ^s`$ with statistical errors less than $`7\%`$. By fitting the temporal and spatial correlation functions to the single-pole ansatz we extracted $`\xi ^t=6.5(4)`$ and $`\xi ^s=8.8(5)`$. The overall error was estimated from the observation of the fit sensitivity to the variation of the parameters and to the range of points used in the fit. We conclude that in order to enter the $`3d`$ scaling region it was not necessary to work so close to the transition itself that the correlation lengths would be large compared to the temporal extent of the lattice. The importance of this comment lies in the fact that the formal derivations of the dimensional reduction scenario require that the zero mode dominates the Matsubara sum and this occurs when the scalar mass is small compared to the temperature. Apparently, our result shows that the higher Matsubara frequencies do not affect the critical singularities but they do affect non-universal aspects of the transition such as the width of the scaling region.

As the temperature approaches $`T_c`$ the mass of the sigma meson decreases to zero and this results in an increase of the mesonic number density. A consequence of this is that the spectral function may become broader due to collisional broadening. One can study in detail the properties of the mesons at $`T0`$ by using the maximum entropy method to extract their spectral functions from high statistics correlation functions measured on anisotropic lattices with temporal sizes larger than $`L_t=12`$. However, for our purposes the $`L_t=12`$ temporal extent is large enough for an estimate of $`\xi ^t`$ at $`\kappa _{\mathrm{cross}}`$. To the extent that we recognized particle-like excitations at $`T0`$ and $`\xi ^t\xi ^s`$, implies that at $`\kappa _{\mathrm{cross}}`$ the breaking of Lorentz invariance caused a significant “renormalization” of the speed of light in the mesonic dispersion relation.

In this paragraph we discuss the numerical results in the low temperature phase of the O(2) symmetric sigma model. We generated approximately $`80,000`$ configurations for each value of the coupling on the $`12\times 24^3`$ lattice and approximately $`40,000`$ configurations for each value of $`\kappa `$ on the $`12\times 48^3`$ lattice. The results for the fits of $`\mathrm{\Sigma }`$ to the standard scaling law for different values of $`\kappa `$ are shown in Table 2. By monitoring the change in the quality of the fits and the deviation in the value of $`\beta _{\mathrm{mag}}`$ from recent estimates of $`\beta _{\mathrm{mag}}^{3\mathrm{d}}=0.3485(2)`$  we conclude that $`\kappa _{\mathrm{cross}}0.3325`$. An interesting observation is that finite size effects are negligible even when the asymmetry ratio is two. In Fig. 2 we plot the order parameter versus coupling. The continuous line represents the power law fit to the $`12\times 48^3`$ data for $`\kappa 0.3275`$ .

The values of the spatial and temporal correlation lengths measured at the crossover coupling $`\kappa _{\mathrm{cross}}=0.3325`$ are $`\xi _{\mathrm{cross}}^s=0.94(2)`$ and $`\xi _{\mathrm{cross}}^t=0.97(3)`$. We conclude that in the $`O(2)`$ model the crossover to the $`3d`$ universality class occurs at a much smaller value of $`\xi /L_t`$ than in the $`O(1)`$ model. This explains why the critical region in this case is more easily accessible on lattices with a smaller asymmetry ratio than in the case of the Ising model. The fact that the value of $`\xi _{\mathrm{cross}}^t`$ is very close to the value of $`\xi _{\mathrm{cross}}^s`$ implies that the thermal fluctuations are relatively small at the crossover into the $`3d`$ scaling region. The diverging correlation length of the Goldstone boson is a source of additional finite size effects in a system. In this case however, the values of the order parameter in the $`3d`$ scaling region measured on small lattices are close to its thermodynamic limit values, because the Goldstone boson had the effect of increasing the width of the critical region.

In Fig. 3 we present the results for the susceptibility $`\chi `$ of the order parameter for the $`O(2)`$ model. As expected finite size effects and statistical errors are much bigger for $`\chi `$ than for $`\mathrm{\Sigma }`$ and for this reason we are unable to extract a meaningful value for the exponent $`\gamma `$.

We also performed simulations of the $`O(3)`$ model on lattices with sizes $`12\times 48^3`$ and $`12\times 72^3`$ and we observed a similar scenario to the one described above for the $`O(2)`$ model. Again, we generated approximately $`30,000`$ to $`50,000`$ configurations for each value of $`\kappa `$ on both lattice sizes. The values of $`\beta _{\mathrm{mag}}`$ extracted from fits near the transition are very close to recent numerical results which predicted $`\beta _{\mathrm{mag}}^{3\mathrm{d}}=0.3689(3)`$ . The discrepancy in the values of $`\beta _{\mathrm{mag}}`$ extracted from the simulations on the two different lattice sizes and the difference in the value of $`\chi ^2/\mathrm{d}.\mathrm{o}.\mathrm{f}.`$ could be attributed to finite size effects on the smaller lattice. As expected the finite size effects in the $`O(3)`$ model are larger than in the $`O(2)`$ model, because it has twice as many Goldstone bosons. Again by monitoring the quality of the fit and the deviation of $`\beta _{\mathrm{mag}}`$ from $`\beta _{\mathrm{mag}}^{3\mathrm{d}}`$ we found that the crossover into the $`3d`$ scaling region occurs at $`\kappa _{\mathrm{cross}}0.4975`$. The values of the spatial and temporal correlation lengths measured on a $`12\times 24^3`$ lattice at the crossover coupling $`\kappa _{\mathrm{cross}}=0.4975`$ are $`\xi _{\mathrm{cross}}^s=0.92(3)`$ and $`\xi _{\mathrm{cross}}^t=0.90(2)`$. Given the uncertainty in the measurement of the $`\kappa _{\mathrm{cross}}`$, we conclude that the crossover in the $`O(2)`$ and $`O(3)`$ models occurs at approximately the same value of the correlation length.

We also performed simulations of the $`(2+1)d`$ Ising model. The results for the order parameter are shown in Fig. 5. In this case, because the dimensionality is lower, the infrared fluctuations are stronger than in the $`(3+1)d`$ model with the same symmetry. Therefore, in order to suppress finite size effects we simulated the model on lattices with larger asymmetry ratios than in the $`(3+1)d`$ case. The values of $`L_s`$ and the amount of statistics generated at various values of $`\kappa `$ are shown in Table 4. By fitting the data with the four smallest values of $`\kappa `$ to the standard scaling relation we extracted $`\kappa _c=0.22394(1)`$ and $`\beta _{\mathrm{mag}}=0.138(8)`$, which is close to the $`2d`$ Ising exponent $`\beta _{\mathrm{mag}}^{2\mathrm{d}}=0.125`$. We conclude that the crossover into the $`2d`$ scaling region is at $`\kappa _{\mathrm{cross}}0.22417`$. The value of the spatial correlation length at $`\kappa _{\mathrm{cross}}`$ measured on a $`12\times 240^2`$ lattice by generating $`5.5`$ million configurations is $`\xi _{\mathrm{cross}}^s=25.0(2.0)`$. Unfortunately, the large statistical fluctuations in the data set did not allow us to measure the temporal correlation length.

Next, we present results in the symmetric phase of the $`O(2)`$ model. We measured the improved estimators for the susceptibility (see eq. (10)), for the exponential correlation length $`\xi _{\mathrm{exp}}`$ extracted from the asymptotic decay of $`C(t)_{\mathrm{imp}}`$ (see eq. (13)) and for the second moment correlation length $`\xi _{2\mathrm{n}\mathrm{d}}`$ (see eq. (15)). It is well-known that in the high temperature phase of this model $`L_s/\xi 7`$ is sufficient to give thermodynamic limit results. In order to detect finite size effects, we performed simulations at $`\kappa =0.2930`$ on lattices with $`L_s=24`$ and $`36`$ and we found that in both cases the values of the susceptibility and the correlation lengths agree within errors, implying that finite-size effects are under control. In most cases we kept $`L_s/\xi >9`$ by increasing $`L_s`$ from $`24`$ to $`96`$ as we approached the critical coupling with constant $`L_t=12`$. For each value of the coupling we generated approximately $`10^7`$ configurations. We also checked that the deviations between $`\xi _{\mathrm{exp}}`$ and $`\xi _{2\mathrm{n}\mathrm{d}}`$ are less than $`0.5\%`$ for all the values of $`\kappa `$. An interesting observation is that $`\xi ^t`$ and $`\xi ^s`$ remain equal for all values of $`\kappa `$ studied in this work. The value of $`\kappa `$ which is closest to $`\kappa _c`$ is $`0.3015`$ and $`\xi =14.18(7)`$ at this value of the coupling. This observation implies that in the symmetric phase the quantum fluctuations remain more important than the thermal fluctuations up to relatively large values of $`\xi /L_t`$.

We also measured the exponents $`\gamma `$ and $`\nu `$ in the $`O(2)`$-symmetric phase by fitting our data to the standard scaling laws

$$\chi =B(\kappa _c\kappa )^\nu ;\xi =C(\kappa _c\kappa )^\gamma .$$
(18)
The results for the fits in different scaling windows are shown in Table 5. The data and the fitting curves are shown in Figs. 6 and 7. The exponents measured are significantly different from recent predictions of the three-dimensional critical exponents $`\gamma ^{3\mathrm{d}}=1.3177(5)`$ and $`\nu ^{3\mathrm{d}}=0.67155(27)`$ . Our results are consistent however, with the observation mentioned above that in the symmetric phase the thermal fluctuations do not influence decisively the critical properties of the model even up to values of the correlation length that are much larger than the value of $`\xi _{\mathrm{cross}}`$ in the broken phase.

## 5 Conclusions

We presented results from numerical simulations of the $`O(N)`$ sigma model at finite temperature for $`N=1,2,3`$ in four dimensions and $`N=1`$ in three dimensions. We investigated the phenomenon of dimensional reduction and studied the dependence of the width of the critical region on the number of field components $`N`$ and on whether the $`O(N)`$ symmetry is spontaneously broken or restored. Contrary to the naive geometric scenario which requires that the condition $`\xi L_t`$ must be satisfied in order to reach the $`d1`$ scaling region, we showed that in the low temperature phase the correlation length of the sigma meson doesn’t have to be much larger than $`L_t`$ in order to reach the dimensionally reduced critical region. Our results imply that the nonzero Matsubara modes do not affect the critical singularities but they can affect non-universal aspects of the transition such as the width of the scaling region. We measured the correlation length of the sigma meson in the temporal and spatial directions at the crossover into the $`d1`$ scaling region and we found that for $`\frac{\mathrm{\Lambda }}{T}=L_t=12`$ in the case of the $`O(2)`$ model $`\xi _{\mathrm{cross}}^t=0.94(2)`$ and $`\xi _{\mathrm{cross}}^s=0.97(3)`$ and for the $`O(3)`$ model $`\xi _{\mathrm{cross}}^t=0.90(2)`$ and $`\xi _{\mathrm{cross}}^s=0.92(3)`$. We should note however, that the deviations from the $`3d`$ scaling laws at these relatively small correlation lengths should also be affected by discretization effects. In the case of the Ising model we measured $`\xi _{\mathrm{cross}}^t=6.5(4)`$ and $`\xi _{\mathrm{cross}}^s=8.8(5)`$, implying that the critical region is much wider in the presence of the infrared fluctuations due to the Goldstone bosons. The critical regions for $`N=2`$ and $`3`$ have approximately the same width and therefore we don’t expect a significant change in the width of the scaling region for $`N>3`$. The fact that for the continuous symmetries $`\xi _{\mathrm{cross}}^s\xi _{\mathrm{cross}}^t`$ implies that the thermal fluctuations are still very small at the crossover temperature.

It is not yet known whether the dimensional reduction scenario is the correct description of the QCD thermal phase transition. If it is, then the thermal transition with massless quarks may belong to the $`O(4)`$ universality class, whereas when the quarks are massive the tricritical point in the $`(\mu ,T)`$ plane becomes critical and it must belong to the $`3d`$ Ising universality class. According to our results the latter must have a much thinner scaling region than the massless quarks transition.

We also showed that the critical region of the $`(3+1)d`$ $`O(2)`$ model in the symmetric phase is much thinner than in the broken phase. Despite the fact that we simulated the model very close to the critical point and went up to correlation lengths $`14`$, the values of the measured exponents are not consistent with the $`3d`$ exponents; they are rather effective exponents with values between mean field exponents and $`3d`$ exponents.

We also studied the symmetry restoration transition in the $`(2+1)d`$ Ising model. This model has a well-defined interacting continuum limit and its behavior near $`T_c`$ is particularly hard to explore analytically, because infrared divergences get worsened in three dimensions. We showed that in order to reach the $`2d`$ scaling region on lattices with temporal extent $`L_t=12`$ the spatial correlation length $`\xi ^s=25.0(2.0)`$ and in order to minimize finite size effects the simulations near $`T_c`$ were performed on lattices with large spatial extents that are larger by almost an order of magnitude than in the four-dimensional case. In this lower dimensionality we studied only $`N=1`$, because according to the Coleman-Mermin-Wagner theorem  for any continuous symmetry in two spatial dimensions the infrared divergent fluctuations wash out the order parameter and the symmetry is manifest for $`T>0`$.

In the future we would like to study the width of the scaling region of the $`O(N)`$ linear sigma model and compare with the results from the nonlinear sigma model presented in this paper. It was shown in the past that the dimensional reduction scenario is the correct description of the chiral phase transition in four-fermion models , because the infrared region is dominated by the zero Matsubara mode of the bosonic field and the nonzero bosonic and fermionic modes decouple. Simulations of the $`(3+1)d`$ four-fermion model with $`Z_2`$ and $`SU(2)\times SU(2)`$ chiral symmetries are underway in order to study the dependence of the width of the critical region on the degree of chiral symmetry of the model.

###### Acknowledgments.

Discussions with Simon Hands and John Kogut are greatly appreciated. Costas Strouthos was supported by a Leverhulme Trust grant. Ioannis Tziligakis was partially supported by the Research Board at the University of Illinois, RES BRD KOGUT J 1-1-28212. Funds provided by this grant were used to purchase a dual processor Pentium IV,2.2 Ghz computer, on which the simulations were carried out. Ioannis Tziligakis would also like to thank Costas Papanikolas and the Institute for Accelerator Systems and Applications (IASA) at the University of Athens, Greece for their hospitality in the fall of 2002.


format: text

text:
# Boson Normal Ordering Problem and Generalized Bell Numbers

## Abstract

For any function $`F(x)`$ having a Taylor expansion we solve the boson normal ordering problem for $`F\left[(a^{})^ra^s\right]`$, with $`r,s`$ positive integers, $`[a,a^{}]=1`$, i.e. we provide exact and explicit expressions for its normal form $`𝒩\left\{F\left[(a^{})^ra^s\right]\right\}=F\left[(a^{})^ra^s\right]`$, where in $`𝒩\left(F\right)`$ all $`a`$’s are to the right. The solution involves integer sequences of numbers which, for $`r,s1`$, are generalizations of the conventional Bell and Stirling numbers whose values they assume for $`r=s=1`$. A complete theory of such generalized combinatorial numbers is given including closed-form expressions (extended Dobinski - type formulas), recursion relations and generating functions. These last are special expectation values in boson coherent states.

Consider a function $`F(x)`$ having a Taylor expansion around $`x=0`$, i.e. $`F(x)=_{k=0}^{\mathrm{}}\frac{F^{(k)}(0)}{k!}x^k`$. In this note we will collect the formulas concerning our solution of the normal ordering problem for $`F\left[(a^{})^ra^s\right]`$, where $`a,a^{}`$ are the boson annihilation and creation operators, $`[a,a^{}]=1`$, and $`r`$ and $`s`$ are positive integers. The normally ordered form of the operator $`F\left[(a^{})^ra^s\right]`$ is denoted by $`𝒩\left\{F\left[(a^{})^ra^s\right]\right\}`$, where in $`𝒩(F)`$ all the $`a`$’s are to the right. It satisfies the operator identity:

$`𝒩\left\{F\left[(a^{})^ra^s\right]\right\}=F\left[(a^{})^ra^s\right].`$ (1)
Furthermore, an auxiliary symbol $`:O(a,a^{}):`$ is used, which means expand $`O`$ in powers of $`a`$ and $`a^{}`$ and order normally *without* taking into account the commutation rule ,.

For the moment we restrict ourselves to the case $`rs`$, the other alternative being treated later. We do not give the proofs of the formulas here; they will be given elsewhere . The case $`r=s=1`$ is known , and some features of the $`r>1,s=1`$ case have been published . When needed we shall refer to , and  where particular cases of our more general formulas are discussed.

We define the set of positive integers $`S_{r,s}(n,k)`$ entering the expansion:

$`\left[(a^{})^ra^s\right]^n=𝒩\left\{\left[(a^{})^ra^s\right]^n\right\}=(a^{})^{n(rs)}\left[{\displaystyle \underset{k=s}{\overset{ns}{}}}S_{r,s}(n,k)(a^{})^ka^k\right],`$ (2)
for $`n=1,2,\mathrm{}`$ . Once the $`S_{r,s}(n,k)`$ are known, the normal ordering of $`\left[(a^{})^ra^s\right]^n`$ is achieved. The same applies to any operator $`F\left[(a^{})^ra^s\right]`$ if $`F(x)`$ has a Taylor expansion around $`x=0`$. The row sums of the triangle $`S_{r,s}(n,k)`$ are given by:

$`B_{r,s}(n)={\displaystyle \underset{k=s}{\overset{ns}{}}}S_{r,s}(n,k),`$ (3)
which for any $`t`$ extend to the polynomials of order $`ns`$ defined by:

$`B_{r,s}(n,t)={\displaystyle \underset{k=s}{\overset{ns}{}}}S_{r,s}(n,k)t^k.`$ (4)
All the subsequent formulas are consequences of the following relation, linking the polynomial of Eq.(4) to a certain infinite series:

$`e^t{\displaystyle \underset{k=s}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k!}}{\displaystyle \underset{j=1}{\overset{n}{}}}[(k+(j1)(rs))(k+(j1)(rs)1)`$ (5)
$`\mathrm{}(k+(j1)(rs)s+1)]t^k=`$ (6)
$`={\displaystyle \underset{k=s}{\overset{ns}{}}}S_{r,s}(n,k)t^k,n=1,2\mathrm{}`$ (7)
which for $`r,s1`$ and $`t=1`$ is an analogue of the celebrated Dobinski relation, which expresses the combinatorial Bell numbers $`B_{1,1}(n)`$ as a sum of an infinite series ,,:

$`B_{1,1}(n)={\displaystyle \frac{1}{e}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{k^n}{k!}}.`$ (8)
For a review of various characteristics of Bell numbers $`B_{1,1}(n)`$ see . It seems that Bell and Stirling numbers are now beginning to appear in textbooks of mathematical physics .

By setting $`t=1`$ in Eq.(5) we obtain equivalent forms of expressions for the generalized Bell numbers $`B_{r,s}(n),n=1,2\mathrm{}`$ :

$`r>s:`$ $`B_{r,s}(n)`$ $`=B_{r,s}(n,1)=`$ (14)
$`={\displaystyle \frac{1}{e}}{\displaystyle \underset{k=s}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k!}}{\displaystyle \underset{j=1}{\overset{n}{}}}[(k+(j1)(rs))(k+(j1)(rs)1)`$
$`\mathrm{}(k+(j1)(rs)s+1)]`$
$`={\displaystyle \frac{1}{e}}{\displaystyle \underset{k=s}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k!}}{\displaystyle \underset{j=1}{\overset{n}{}}}\left[k+(j1)(rs)\right]^{\underset{¯}{s}}`$
$`={\displaystyle \frac{1}{e}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k!}}{\displaystyle \underset{j=1}{\overset{n1}{}}}{\displaystyle \frac{(k+jr(j1)s)!}{(k+j(rs))!}}`$
$`={\displaystyle \frac{(rs)^{s(n1)}}{e}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k!}}\left[{\displaystyle \underset{j=1}{\overset{s}{}}}{\displaystyle \frac{\mathrm{\Gamma }(n+\frac{k+j}{rs})}{\mathrm{\Gamma }(1+\frac{k+j}{rs})}}\right],`$
$`r=s:`$ $`B_{r,r}(n)`$ $`=B_{r,r}(n,1)={\displaystyle \frac{1}{e}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k!}}\left[{\displaystyle \frac{(k+r)!}{k!}}\right]^{n1}.`$ (16)
For all $`r,s`$ we set $`B_{r,s}(0)=1`$ by convention. In Eq.(14) we have employed the notation $`m^{\underset{¯}{s}}:=m(m1)(m2)\mathrm{}(ms+1)`$ for the falling factorial , and in Eq.(14) $`\mathrm{\Gamma }(y)`$ is the Euler gamma function.

The numbers $`S_{r,s}(n,k)`$ are non-zero for $`skns`$, and satisfy by convention $`S_{r,s}(n,0)=\delta _{n,0}`$. We shall refer to them as generalized Stirling numbers of the second kind. Their exact expressions in the form of a finite sum are:

$`S_{r,s}(n,k)`$ $`=`$ $`{\displaystyle \frac{(1)^k}{k!}}{\displaystyle \underset{p=s}{\overset{k}{}}}(1)^p\left(\begin{array}{c}k\\ p\end{array}\right){\displaystyle \underset{j=1}{\overset{n}{}}}[(p+(j1)(rs))(p+(j1)(rs)1)`$ (20)
$`\mathrm{}(p+(j1)(rs)s+1)]`$
$`=`$ $`{\displaystyle \frac{(1)^k}{k!}}{\displaystyle \underset{p=s}{\overset{k}{}}}(1)^p\left(\begin{array}{c}k\\ p\end{array}\right){\displaystyle \underset{j=1}{\overset{n}{}}}\left(p+(j1)(rs)\right)^{\underset{¯}{s}}`$ (23)
$`=`$ $`{\displaystyle \frac{(1)^k}{k!}}\left\{\left(x^r{\displaystyle \frac{d^s}{dx^s}}\right)^n\left[(1x)^k{\displaystyle \underset{p=0}{\overset{s1}{}}}\left(\begin{array}{c}k\\ p\end{array}\right)(x)^p\right]\right\}_{x=1},`$ (26)
from which for $`r=s=1`$ one obtains the standard form of the classical Stirling numbers of the second kind :

$`S_{1,1}(n,k)={\displaystyle \frac{(1)^k}{k!}}{\displaystyle \underset{p=1}{\overset{k}{}}}(1)^p\left(\begin{array}{c}k\\ p\end{array}\right)p^n,`$ (29)
and for $`r=2,s=1`$,

$`S_{2,1}(n,k)={\displaystyle \frac{n!}{k!}}\left(\begin{array}{c}n1\\ k1\end{array}\right),`$ (32)
which are the so-called unsigned Lah numbers ,.

That the numbers $`S_{r,s}(n,k)`$ are natural extensions of $`S_{1,1}(n,k)`$ can be neatly seen by observing their action on the space of polynomials generated by falling factorials. Let $`x`$ be any real number. Then,

$`{\displaystyle \underset{j=1}{\overset{n}{}}}(x+(j1)(rs))(x+(j1)(rs)1)\mathrm{}(x+(j1)(rs)s+1)=`$ (33)
$`={\displaystyle \underset{k=s}{\overset{ns}{}}}S_{r,s}(n,k)x(x1)\mathrm{}(xk+1),`$ (34)
which, when rewritten with $`m^{\underset{¯}{s}}`$ symbols become:

$`{\displaystyle \underset{j=1}{\overset{n}{}}}(x+(j1)(rs))^{\underset{¯}{s}}={\displaystyle \underset{k=s}{\overset{ns}{}}}S_{r,s}(n,k)x^{\underset{¯}{k}}.`$ (35)
This last equation when specified to $`r=s`$ gives a particularly transparent interpretation of $`S_{r,r}(n,k)`$ as connection coefficients:

$`[x^{\underset{¯}{r}}]^n={\displaystyle \underset{k=r}{\overset{nr}{}}}S_{r,r}(n,k)x^{\underset{¯}{k}}.`$ (36)
By choosing $`r=s=1`$ it boils down to the defining equations for $`S_{1,1}(n,k)`$ in terms of $`x^{\underset{¯}{k}}`$, used by J. Stirling himself :

$`x^n={\displaystyle \underset{k=1}{\overset{n}{}}}S_{1,1}(n,k)x^{\underset{¯}{k}}.`$ (37)
The relations for $`S_{r,s}(n,k)`$ become particularly appealing for $`r=s`$: Eq.(23) can be reformulated as:

$`S_{r,r}(n,k)`$ $`=`$ $`{\displaystyle \frac{(1)^k}{k!}}{\displaystyle \underset{p=r}{\overset{k}{}}}(1)^p\left(\begin{array}{c}k\\ p\end{array}\right)\left[p(p1)\mathrm{}(pr+1)\right]^n`$ (40)
$`=`$ $`{\displaystyle \frac{(1)^k}{k!}}{\displaystyle \underset{p=r}{\overset{k}{}}}(1)^p\left(\begin{array}{c}k\\ p\end{array}\right)\left[p^{\underset{¯}{r}}\right]^n,`$ (43)
which differs from Eq.(29) in that the sum in Eq.(43) starts from $`p=r`$ and $`p^n`$ becomes $`[p^{\underset{¯}{r}}]^n`$. The recurence relations for $`S_{r,r}(n,k)`$ are the following:

$`S_{r,r}(1,r)=1,S_{r,r}(n,k)=0,k<r,nr<k(n+1)r,`$ (44)
$`S_{r,r}(n+1,k)={\displaystyle \underset{p=0}{\overset{r}{}}}\left(\begin{array}{c}k+pr\\ p\end{array}\right)r^{\underset{¯}{p}}S_{r,r}(n,k+pr),rknr,n>1,`$ (47)
where the definition $`r^{\underset{¯}{0}}:=1`$ is made. Note that for $`r=s=1`$ we get the known relation for conventional Stirling numbers i.e.: $`S_{1,1}(n+1,k)=kS_{1,1}(n,k)+S_{1,1}(n,k1)`$ with appropriate initial conditions .

The “non-diagonal” generalized Bell numbers $`B_{r,s}(n)`$ can always be expressed as special values of generalized hypergeometric functions $`{}_{p}{}^{}F_{q}^{}`$. Algebraic manipulation of Eqs.(14)-(16) yields the following examples:

$`r>1,s=1`$:

$`B_{r,1}(n)`$ is a combination of $`r1`$ different hypergeometric functions of type $`{}_{1}{}^{}F_{r1}^{}(\mathrm{};x)`$, each of them evaluated at the same value of argument $`x=(r1)^{1r}`$; here are some lowest order cases:

$`B_{2,1}(n)`$ $`=`$ $`{\displaystyle \frac{n!}{e}}{}_{1}{}^{}F_{1}^{}(n+1;2;1)=(n1)!L_{n1}^{(1)}(1),`$ (48)
$`B_{3,1}(n)`$ $`=`$ $`{\displaystyle \frac{2^{n1}}{e}}\left({\displaystyle \frac{2\mathrm{\Gamma }(n+\frac{1}{2})}{\sqrt{\pi }}}{}_{1}{}^{}F_{2}^{}(n+{\displaystyle \frac{1}{2}};{\displaystyle \frac{1}{2}},{\displaystyle \frac{3}{2}};{\displaystyle \frac{1}{4}})+n!{}_{1}{}^{}F_{2}^{}(n+1;{\displaystyle \frac{3}{2}},2;{\displaystyle \frac{1}{4}})\right),`$ (49)
$`B_{4,1}(n)`$ $`=`$ $`{\displaystyle \frac{3^{n1}}{2e}}({\displaystyle \frac{3^{3/2}\mathrm{\Gamma }(\frac{2}{3})\mathrm{\Gamma }(n+\frac{1}{3})}{\pi }}{}_{1}{}^{}F_{3}^{}(n+{\displaystyle \frac{1}{3}};{\displaystyle \frac{1}{3}},{\displaystyle \frac{2}{3}},{\displaystyle \frac{4}{3}};{\displaystyle \frac{1}{27}})+`$ (51)
$`{\displaystyle \frac{3\mathrm{\Gamma }(n+\frac{2}{3})}{\mathrm{\Gamma }(\frac{2}{3})}}{}_{1}{}^{}F_{3}^{}(n+{\displaystyle \frac{2}{3}};{\displaystyle \frac{2}{3}},{\displaystyle \frac{4}{3}},{\displaystyle \frac{5}{3}};{\displaystyle \frac{1}{27}})+n!{}_{1}{}^{}F_{3}^{}(n+1;{\displaystyle \frac{4}{3}},{\displaystyle \frac{5}{3}},2;{\displaystyle \frac{1}{27}})),`$
etc. In Eq.(48) $`L_m^{(\alpha )}(y)`$ is the associated Laguerre polynomial.

Similarly the series $`B_{2r,r}(n)`$ can be written down in a compact form using the confluent hypergeometric function of Kummer:

$`B_{2r,r}(n)={\displaystyle \frac{(rn)!}{er!}}{}_{1}{}^{}F_{1}^{}(rn+1,r+1;1).`$ (52)
A still more general family of sequences arising from Eqs.(14)-(16) has the form ($`p,r=1,2\mathrm{}`$):

$`B_{pr+p,pr}(n)`$ $`=`$ $`{\displaystyle \frac{1}{e}}\left[{\displaystyle \underset{j=1}{\overset{r}{}}}{\displaystyle \frac{(p(n1)+j)!}{(pj)!}}\right]`$ (55)
$`{}_{r}{}^{}F_{r}^{}(pn+1,pn+1+p,\mathrm{},pn+1+p(r1);1+p,1+2p\mathrm{},1+rp;1),`$
etc.

In contrast, the ”diagonal” numbers $`B_{r,r}(n)`$ of Eq.(16), which also can be rewritten as:

$`B_{r,r}(n)={\displaystyle \frac{1}{e}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(k+r1)!}}\left[k(k+1)\mathrm{}(k+r1)\right]^n,n=1,2,\mathrm{}`$ (56)
cannot be expressed through hypergeometric functions. However, the sequences $`B_{r,r}(n)`$ possess a particularity that they can be always expressed in terms of conventional Bell numbers with $`r`$-nomial (binomial, trinomial…) coefficients. For example:

$`B_{2,2}(n)={\displaystyle \underset{k=0}{\overset{n1}{}}}\left(\begin{array}{c}n1\\ k\end{array}\right)B_{1,1}(n+k).`$ (59)
Several low-order triangles of $`S_{r,s}(n,k)`$ and their associated $`B_{r,s}(n)`$ are presented in Table 1.

It turns out that various generating functions for $`B_{r,s}(n)`$ and $`S_{r,s}(n,k)`$ can be related to special quantum states, called coherent states , defined as linear combinations of the eigenstates of the harmonic oscillator, $`H=a^{}a,H|n=n|n,n|n^{}=\delta _{n,n^{}}`$ and defined as:

$`|z=e^{\frac{|z|^2}{2}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{z^n}{\sqrt{n!}}}|n,`$ (60)
(with $`z|z=1`$), for $`z`$ complex. The states $`|z`$ satisfy:

$`a|z=z|z.`$ (61)
It has been noticed in  that,

$`z|(a^{}a)^n|z\stackrel{|z|=1}{=}B_{1,1}(n)`$ (62)
and the exponential generating function (egf) of the numbers $`B_{1,1}(n)`$, i.e. $`_{n=0}^{\mathrm{}}B_{1,1}(n)\frac{\lambda ^n}{n!}`$ satisfies

$`z|e^{\lambda a^{}a}|s=z|:e^{a^{}a(e^\lambda 1)}:|z\stackrel{|z|=1}{=}e^{e^\lambda 1}={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}B_{1,1}(n){\displaystyle \frac{\lambda ^n}{n!}}`$ (63)
which is a restatement of the known fact that ,,:

$`𝒩\left(e^{\lambda a^{}a}\right)=:e^{a^{}a(e^\lambda 1)}:.`$ (64)
How does it generalize to $`r,s1`$ ?

Eqs.(2) and (3) give directly:

$`z|\left[(a^{})^ra^s\right]^n|z\stackrel{z=1}{=}B_{r,s}(n),`$ (65)
valid for all $`r,s`$.

We first treat the case $`r=s=1,2\mathrm{}`$. Note in this context that a Hermitian Hamiltonian $`(a^{})^ra^r`$ is of great importance in quantum optics, as for $`r=2,3,\mathrm{}`$ it provides a description of non-linear Kerr-type media. Using Eq.(43) we can find the following egf of $`S_{r,r}(n,k)`$:

$`{\displaystyle \underset{n=k/r}{}}{\displaystyle \frac{x^n}{n!}}S_{r,r}(n,k)={\displaystyle \frac{(1)^k}{k!}}{\displaystyle \underset{p=r}{\overset{k}{}}}(1)^p\left(\begin{array}{c}k\\ p\end{array}\right)\left(e^{xp(p1)\mathrm{}(pr+1)}1\right),`$ (68)
where $`y`$ is the ceiling function , defined as the nearest integer greater or equal to $`y`$. This yields via exchange of order of summation:

$`𝒩\left(e^{\lambda (a^{})^ra^r}\right)=1+:{\displaystyle \underset{k=r}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^k}{k!}}\left[{\displaystyle \underset{p=r}{\overset{k}{}}}(1)^p\left(\begin{array}{c}k\\ p\end{array}\right)\left(e^{\lambda p(p1)\mathrm{}(pr+1)}1\right)\right](a^{}a)^k:`$ (71)
and

$`z|e^{\lambda (a^{})^ra^r}|z\stackrel{|z|=1}{=}1+{\displaystyle \underset{k=r}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(1)^k}{k!}}\left[{\displaystyle \underset{p=r}{\overset{k}{}}}(1)^p\left(\begin{array}{c}k\\ p\end{array}\right)\left(e^{\lambda p^{\underset{¯}{r}}}1\right)\right].`$ (74)
Another case which can be written in a closed form is $`r>1`$, $`s=1`$, for which the egf for $`S_{r,1}(n,k)`$ reads:

$`{\displaystyle \underset{n=k/r}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^n}{n!}}S_{r,1}(n,k)={\displaystyle \frac{1}{k!}}\left\{\left(1(r1)x\right)^{\frac{1}{r1}}1\right\}.`$ (75)
One then obtains ,:

$`𝒩\left(e^{\lambda (a^{})^ra}\right)=:exp\left\{[(1\lambda (a^{})^{r1}(r1))^{\frac{1}{r1}}1]a^{}a\right\}:`$ (76)
For arbitrary $`r>s`$ the following formula is still valid:

$`𝒩\left(e^{\lambda (a^{})^ra^s}\right)=1+:{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\lambda ^n}{n!}}(a^{})^{n(rs)}\left({\displaystyle \underset{k=s}{\overset{ns}{}}}S_{r,s}(n,k)(a^{}a)^k\right):`$ (77)
where the explicit form of $`S_{r,s}(n,k)`$ may be used, see Eqs.(20)-(26).

A close look at Eqs.(74) and (65) reveals that

$`z|e^{\lambda (a^{})^ra^r}|z\stackrel{|z|=1}{=}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}B_{r,r}(n){\displaystyle \frac{\lambda ^n}{n!}}.`$ (78)
Comparing Eqs.(78),(76) and (74) we conclude that for $`r>1,s=1`$ and $`r=s`$ the egf of respective generalized Bell numbers are special matrix elements of $`e^{\lambda (a^{})^ra^s}`$ in coherent states.

However, for arbitrary $`r>s>1`$ examination of Eqs.(14)-(14) confirms that the numbers $`B_{r,s}(n)`$ increase so rapidly with $`n`$ that one cannot define their egf’s meaningfully. In particular, this is reflected by Eq.(77) which is true in its operator form, but one needs to exercise care when calculating its matrix elements as the convergence of the result may require limitations on $`\lambda `$.

A well-defined and convergent procedure for such sequences is to consider what we call *hypergeometric generating functions* (hgf), i.e. the egf for the ratios $`B_{r,s}/(n!)^t`$, where $`t`$ is an appropriately chosen integer. A case in point is the series $`B_{3,2}(n)`$ which may be written explicitly from Eq.(55) as:

$`B_{3,2}(n)={\displaystyle \frac{1}{e}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(n+k)!(n+k+1)!}{k!(k+1)!(k+2)!}}.`$ (79)
Its hgf $`G_{3,2}(\lambda )`$ is then:

$`G_{3,2}(\lambda )={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left[{\displaystyle \frac{B_{3,2}(n)}{n!}}\right]{\displaystyle \frac{\lambda ^n}{n!}}={\displaystyle \frac{1}{e}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(k+2)!}}{}_{2}{}^{}F_{1}^{}(k+2,k+1;1;\lambda ).`$ (80)
Similarly for $`G_{4,2}(\lambda )`$ one obtains:

$`G_{4,2}(\lambda )={\displaystyle \frac{1}{e}}{\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(k+2)!}}{}_{2}{}^{}F_{1}^{}({\displaystyle \frac{k+2}{2}},{\displaystyle \frac{k}{2}}+1;1;4\lambda ).`$ (81)
Eq.(52) implies more generally:

$`G_{2r,r}(\lambda )={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left[{\displaystyle \frac{B_{2r,r}(n)}{(n!)^{r1}}}\right]{\displaystyle \frac{\lambda ^n}{n!}}`$ (82)
$`=\{\begin{array}{cc}\frac{1}{1!e}_{k=0}^{\mathrm{}}\frac{1}{(k+1)!}{}_{2}{}^{}F_{1}^{}(\frac{k+1}{2},\frac{k}{2}+1;1;4\lambda ),\hfill & r=2,\hfill \\ \frac{1}{2!e}_{k=0}^{\mathrm{}}\frac{1}{(k+2)!}{}_{3}{}^{}F_{2}^{}(\frac{k+1}{3},\frac{k+2}{3},\frac{k+3}{3};1,1;27\lambda ),\hfill & r=3,\hfill \\ \frac{1}{3!e}_{k=0}^{\mathrm{}}\frac{1}{(k+3)!}{}_{4}{}^{}F_{3}^{}(\frac{k+1}{4},\mathrm{},\frac{k+4}{4};1,1,1;256\lambda ),\hfill & r=4\mathrm{}\text{etc.}\hfill \end{array}`$ (86)
See  for other instances where this type of hgf appear.

All the previous considerations also apply to the case $`r<s`$. In this case we define the generalized Stirling numbers by (note the difference with Eq.(2)):

$`rs:\left[(a^{})^ra^s\right]^n=𝒩\left\{\left[(a^{})^ra^s\right]^n\right\}=\left[{\displaystyle \underset{k=r}{\overset{nr}{}}}S_{r,s}(n,k)(a^{})^ka^k\right]a^{n(sr)},`$ (87)
By taking the Hermitian conjugate of Eq.(2), with the change $`rs`$, we obtain

$`rs:\left[(a^{})^ra^s\right]^n=\left[{\displaystyle \underset{k=r}{\overset{nr}{}}}S_{s,r}(n,k)(a^{})^ka^k\right]a^{n(sr)},`$ (88)
whence the symmetry relation:

$`S_{r,s}(n,k)=S_{s,r}(n,k),skns,(rs).`$ (89)
The normal ordering of $`F[a^s(a^{})^r]`$ involves the antinormal to normal order, in the terminology of . With this in mind we define the *anti*-Stirling numbers of the second kind $`\stackrel{~}{S}_{r,s}(n,k)`$ as:

$`rs:[a^s(a^{})^r]^n=(a^{})^{n(rs)}\left[{\displaystyle \underset{k=0}{\overset{ns}{}}}\stackrel{~}{S}_{s,r}(n,k)(a^{})^ka^k\right].`$ (90)
It can be demonstrated  that they satisfy the following symmetry properties:

$`\stackrel{~}{S}_{s,r}(n,k)=\stackrel{~}{S}_{r,s}(n,k),0kns,(rs),`$ (91)
$`\stackrel{~}{S}_{r,s}(n,k)=S_{r,s}(n+1,k+s),0kns,(rs).`$ (92)
We have obtained recursion relations for $`B_{r,s}(n)`$ for certain values of $`r,s`$. They include the $`r1,s=1`$ case:

$`B_{r,1}(n+1)={\displaystyle \underset{k=0}{\overset{n}{}}}\left(\begin{array}{c}n\\ k\end{array}\right)\left[{\displaystyle \underset{j=0}{\overset{nk}{}}}(r+(j1)(r1))\right]B_{r,1}(k),`$ (95)
which for $`r=2,s=1`$ may be written explicitly as:

$`B_{2,1}(n+1)={\displaystyle \underset{k=0}{\overset{n}{}}}\left(\begin{array}{c}n\\ k\end{array}\right)(nk+1)!B_{2,1}(k),`$ (98)
and for $`r=s=1`$ reduces to the known relation ,, which is the binomial transform:

$`B_{1,1}(n+1)={\displaystyle \underset{k=0}{\overset{n}{}}}\left(\begin{array}{c}n\\ k\end{array}\right)B_{1,1}(k).`$ (101)
If in Eq.(1) $`F(x)`$ has a Taylor expansion of $`F(x)`$ around $`x_00`$, then

$`𝒩\left\{F\left[(a^{})^ra^s\right]\right\}=𝒩\left({\displaystyle \underset{k=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{F^{(k)}(x_o)}{k!}}\left((a^{})^ra^sx_o\right)^k\right),`$ (102)
and our above results should be supplemented by the appropriate binomial expansion coefficients.

All the $`B_{r,s}(n)`$ described here are the $`n`$-th moments of positive functions on the positive half axis; some solutions of the associated Stieltjes moment problem have been recently discussed at some length . In addition we have found that even larger classes of combinatorial sequences are solutions of the moment problem and they can be used for the construction of new families of coherent states .

We are considering the possible combinatorial interpretation of the above results.
Acknowledgments

We thank G. Duchamp, L. Haddad, A. Horzela, M. Mendez and J.Y. Thibon for enlightening discussions.
N.J.A. Sloane’s Encyclopedia of Integer Sequences (http://www.research.att.com/~njas/sequences) was an essential help in sorting out the properties of sequences we have encountered in this work.


text:
# 1 Twist-3 DVCS amplitude

## 1 Twist-3 DVCS amplitude

Deeply Virtual Compton Scattering (DVCS)  is a subject of rather intensive investigations during the last few years (for a detailed recent review of DVCS see ). It is a two-photon process:

$$\gamma ^{}(q)+N(P)\gamma ^{}(q^{})+N^{}(P^{}),$$
(1)
in which the initial virtual photon $`\gamma ^{}`$, being produced by a lepton, interacts with a parton of the initial nucleon $`N`$, the final real photon $`\gamma ^{}`$ is radiated by the parton, and the nucleon $`N^{}`$ in the final state is formed by the active parton and the rest of the initial nucleon $`N`$. The reaction is considered in the Bjorken limit, when the virtuality of the initial photon $`Q^2=q^2P^2`$ is large, Bjorken variable $`x_{\text{Bj}}=Q^2/2(Pq)`$ is fixed and the momentum transfer squared is small: $`(P^{}P)^2Q^2`$. The first experimental results on DVCS can be found in Refs. .

We will work in the infinite momentum frame defined by the light-cone vectors $`n^\mu =1/(\overline{P}_+\sqrt{2})(1,0,0,1)`$ and $`p^\mu =\overline{P}_+/\sqrt{2}(1,0,0,1)`$, with $`\overline{P}^\mu =(P+P^{})^\mu /2`$ being the average hadron momentum (note, that $`n^2=p^2=0`$ and $`np=n\overline{P}=1`$). Let us define the momentum transfer as $`\mathrm{\Delta }^\mu =(P^{}P)^\mu `$ and its square as $`t=(P^{}P)^2`$. In this paper we shall use a kinematical variable $`\xi `$, analogous to the Bjorken one, defined as $`\xi =(n\mathrm{\Delta })/2`$. A more detailed description of the DVCS kinematics can be found, e.g., in .

To leading order in the electromagnetic coupling constant, the amplitude of DVCS is proportional to the hadronic tensor $`H_{\mu \nu }`$,

$$H_{\mu \nu }=id^4xe^{iqx}P^{}|TJ_\mu (x)J_\nu (0)|P.$$
(2)
To order $`(1/Q)^0`$, it can be expressed through the non-diagonal hadron matrix elements $`_\alpha `$ and $`\stackrel{~}{}_\alpha `$ of gauge-invariant light-cone bilocal operators:

$`H_{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(g_{\mu \nu }^{}\right){\displaystyle _1^1}𝑑xC_+(x,\xi )n^\alpha _\alpha (x,\xi ,t)`$ (3)
$`{\displaystyle \frac{i}{2}}\epsilon _{\mu \nu }^{}{\displaystyle _1^1}𝑑xC_{}(x,\xi )n^\alpha \stackrel{~}{}_\alpha (x,\xi ,t)+𝒪\left({\displaystyle \frac{1}{Q}}\right),`$
where $`_\alpha `$ and $`\stackrel{~}{}_\alpha `$ are defined as

$`_\alpha (x,\xi ,t)`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{d\lambda }{2\pi }}e^{i\lambda x}P^{},S^{}\left|\overline{\psi }\left({\displaystyle \frac{\lambda }{2}}n\right)\gamma _\alpha \psi \left({\displaystyle \frac{\lambda }{2}}n\right)\right|P,S,`$ (4)
$`\stackrel{~}{}_\alpha (x,\xi ,t)`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{d\lambda }{2\pi }}e^{i\lambda x}P^{},S^{}\left|\overline{\psi }\left({\displaystyle \frac{\lambda }{2}}n\right)\gamma _\alpha \gamma _5\psi \left({\displaystyle \frac{\lambda }{2}}n\right)\right|P,S`$ (5)
(let us note, that here we do not write explicitly the phase factors which make the bilocal operators to be those gauge-invariant). The coefficient functions $`C_\pm (x,\xi )`$ are

$$C_\pm (x,\xi )=\frac{1}{x\xi +i\epsilon }\pm \frac{1}{x+\xi i\epsilon },$$
(6)
and the transverse tensors $`g_{\mu \nu }^{}`$ and $`\epsilon _{\mu \nu }^{}`$ are

$$g_{\mu \nu }^{}=g_{\mu \nu }+n_\mu p_\nu +n_\nu p_\mu ,\epsilon _{\mu \nu }^{}=\epsilon _{\mu \nu \lambda \sigma }n^\lambda p^\sigma .$$
(7)
In general, the functions $`_\alpha `$ and $`\stackrel{~}{}_\alpha `$ contain both the purely quark and mixed quark-gluon contributions. However, in the considered case, the quark-gluon contribution vanishes due to contraction with the light-cone vector $`n^\alpha `$. Therefore the scalar products $`n^\alpha _\alpha `$ and $`n^\alpha \stackrel{~}{}_\alpha `$ can be parameterized in terms of four functions $`H_q(x,\xi ,t)`$, $`E_q(x,\xi ,t)`$ and $`\stackrel{~}{H}_q(x,\xi ,t)`$, $`\stackrel{~}{E}_q(x,\xi ,t)`$ in the following way :

$`n^\alpha _\alpha `$ $`=`$ $`H_q(x,\xi ,t)\overline{N}(P^{},S^{})\widehat{n}N(P,S)+E_q(x,\xi ,t)\overline{N}(P^{},S^{}){\displaystyle \frac{i\sigma ^{\mu \nu }n_\mu \mathrm{\Delta }_\nu }{2M}}N(P,S),`$ (8)
$`n^\alpha \stackrel{~}{}_\alpha `$ $`=`$ $`\stackrel{~}{H}_q(x,\xi ,t)\overline{N}(P^{},S^{})\widehat{n}\gamma _5N(P,S)+\stackrel{~}{E}_q(x,\xi ,t)\overline{N}(P^{},S^{}){\displaystyle \frac{\gamma _5\mathrm{\Delta }n}{2M}}N(P,S).`$ (9)
Here $`N(P,S)`$ is the Dirac bispinor of a hadron (a proton in our case), $`M`$ is the hadron mass and we use the notation $`\widehat{x}=x^\mu \gamma _\mu `$ throughout the paper. The functions $`H`$, $`E`$ and $`\stackrel{~}{H}`$, $`\stackrel{~}{E}`$ are the Generalized Parton Distributions (GPD) (for a review see Refs. ). In the forward limit $`P^{}=P`$ the distributions $`H_q`$ and $`\stackrel{~}{H}_q`$ coincide with the usual parton densities $`q(x)`$ and $`\mathrm{\Delta }q(x)`$. The distributions $`E_q`$ and $`\stackrel{~}{E}_q`$ are completely new ones. Thus the investigation of DVCS allows us to get an additional information about the structure of hadrons. In particular, the GPDs contain an information about the hadron spin distribution among the partons of the hadron .

The amplitude of DVCS has corrections of two types: firstly, loop corrections in powers of the electromagnetic and strong coupling constants and, secondly, higher twist corrections in powers of $`1/Q`$. We will not consider the loop corrections in this paper. Making allowance for the first correction in powers of $`1/Q`$, the hadronic tensor $`H_{\mu \nu }`$ takes the following form  (for the pion target see ):

$`H_{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\left(g_{\mu \nu }^{}\right){\displaystyle \frac{\mathrm{\Delta }_\mu ^{}\overline{P}_\nu }{\left(\overline{P}q^{}\right)}}\right]{\displaystyle _1^1}𝑑xC_+(x,\xi )n^\alpha _\alpha (x,\xi ,t)`$ (10)
$`+`$ $`{\displaystyle \frac{i}{2}}\epsilon _\mu ^\lambda \left[\left(g_{\lambda \nu }^{}\right){\displaystyle \frac{\mathrm{\Delta }_\lambda ^{}\overline{P}_\nu }{\left(\overline{P}q^{}\right)}}\right]{\displaystyle _1^1}𝑑xC_{}(x,\xi )n^\alpha \stackrel{~}{}_\alpha (x,\xi ,t)`$
$``$ $`{\displaystyle \frac{(q+4\xi \overline{P})_\mu }{2\left(\overline{P}q^{}\right)}}[(g_{\lambda \nu }^{}){\displaystyle \frac{\mathrm{\Delta }_\lambda ^{}\overline{P}_\nu }{\left(\overline{P}q^{}\right)}}]{\displaystyle _1^1}dx\{_{}^\lambda (x,\xi ,t)C_+(x,\xi )`$
$`i\epsilon _{}^{\lambda \sigma }\stackrel{~}{}_\sigma ^{}(x,\xi ,t)C_{}(x,\xi )\}`$
$`+`$ $`𝒪\left({\displaystyle \frac{1}{Q^2}}\right).`$
Note that $`H_{\mu \nu }`$ satisfies the transversity conditions

$`q^\mu H_{\mu \nu }=0,H_{\mu \nu }(q^{})^\nu =0.`$ (11)
We see that the correction contains a term proportional to the transverse part of the vector functions $`_\alpha `$ and $`\stackrel{~}{}_\alpha `$. This means that, in comparison to the leading order, to order $`1/Q`$ we deal with corrections of two types: those “kinematical” associated with purely quark operators and “genuine twist-3” ones originated from the quark-gluon operators, which have the form $`\overline{\psi }G\psi `$, where $`\psi `$ is the operator of the quark field and $`G`$ denotes the operator of the gluon field strength tensor. The “kinematical” contribution can be expressed through the twist-2 quark GPDs $`H`$, $`E`$, $`\stackrel{~}{H}`$, $`\stackrel{~}{E}`$ .

In the case of polarized DIS the quark-gluon correction is supposed to be small relative to the “kinematical” one. This assumption is known as Wadzura-Wilczek approximation . The experimental data  point at its reliability. It is assumed the same approximation to be valid in the case of DVCS. The aim of this paper is to test this hypothesis in the model of the instanton vacuum .

A method based on the model of the instanton vacuum, which allows to calculate hadronic matrix elements of the quark-gluon operators, was suggested . In the case of DIS it was shown  that the operators of twist-3 are of order $`(\rho /R)^4\mathrm{ln}(\rho /R)`$, where $`\rho `$ is an average size of instantons and $`R`$ is an average distance between instantons in the instanton medium. Thus, the matrix elements are parametrically small due to the packing fraction $`\rho /R1/3`$. In this paper we apply this analysis to the nondiagonal hadronic matrix elements of the quark-gluon operators in the framework of DVCS. We show that in the model of the instanton vacuum they are parametrically suppressed by the packing fraction of the instanton vacuum. In the forward limit we reproduce the results for the diagonal case, that is for polarized DIS.

## 2 “Genuine twist-3” effects

Let us consider the following vector and axial vector bilocal operators

$`V_\alpha (x,x)`$ $`=`$ $`\overline{\psi }(x)\gamma _\alpha [x,x]\psi (x),`$ (12)
$`A_\alpha (x,x)`$ $`=`$ $`\overline{\psi }(x)\gamma _\alpha \gamma _5[x,x]\psi (x).`$ (13)
Here

$$[x,y]=P\mathrm{exp}\left\{ig_0^1𝑑u(xy)^\nu A_\nu (ux+(1u)y)\right\}$$
(14)
is the phase factor ordered along a straight line connecting the points $`x`$ and $`y`$. In what follows we will not write it explicitly assuming its presence in all non-local gauge invariant operators. In oder to separate the “genuine twist-3” contributions to the operators (12), (13) from the “kinematical” ones we shall use a technique elaborated in . According to the result obtained by Balitsky and Braun  in the case of massless fermion operators there are the following operator identities ($`\epsilon _{0123}=+1`$):

$`\overline{\psi }(x)\gamma _\alpha \psi (x)`$ $`=`$ $`{\displaystyle \frac{}{x^\alpha }}{\displaystyle _0^1}𝑑u\overline{\psi }(ux)\widehat{x}\psi (ux)`$ (15)
$``$ $`i\epsilon _{\alpha \beta \lambda \sigma }{\displaystyle _0^1}𝑑uux^\beta 𝔇^\lambda \left[\overline{\psi }(ux)\gamma ^\sigma \gamma _5\psi (ux)\right]+{\displaystyle _0^1}𝑑u𝔊_\alpha (ux),`$
$`\overline{\psi }(x)\gamma _\alpha \gamma _5\psi (x)`$ $`=`$ $`{\displaystyle \frac{}{x^\alpha }}{\displaystyle _0^1}𝑑u\overline{\psi }(ux)\widehat{x}\gamma _5\psi (ux)`$ (16)
$``$ $`i\epsilon _{\alpha \beta \lambda \sigma }{\displaystyle _0^1}𝑑uux^\beta 𝔇^\lambda \left[\overline{\psi }(ux)\gamma ^\sigma \psi (ux)\right]+{\displaystyle _0^1}𝑑u𝔊_{5\alpha }(ux),`$
where

$`𝔊_\alpha (x)`$ $`=`$ $`g{\displaystyle _1^1}𝑑v\overline{\psi }(x)\left[\epsilon _{\alpha \mu \nu \tau }G^{\mu \rho }(vx)x_\rho x^\nu \gamma ^\tau \gamma _5ivG_{\alpha \rho }(vx)x^\rho \widehat{x}\right]\psi (x),`$ (17)
$`𝔊_{5\alpha }(x)`$ $`=`$ $`g{\displaystyle _1^1}𝑑v\overline{\psi }(x)\left[\epsilon _{\alpha \mu \nu \tau }G^{\mu \rho }(vx)x_\rho x^\nu \gamma ^\tau \gamma _5ivG_{\alpha \rho }(vx)x^\rho \widehat{x}\right]\gamma _5\psi (x).`$ (18)
In these equations the derivative $`𝔇^\lambda `$ is defined for any combination $`\mathrm{\Gamma }`$ of Dirac $`\gamma `$-matrices as

$$𝔇^\lambda \left\{\overline{\psi }(x)\mathrm{\Gamma }[x,x]\psi (x)\right\}=\underset{y0}{lim}\frac{}{y_\lambda }\left\{\overline{\psi }(x+y)\mathrm{\Gamma }[x+y,x+y]\psi (x+y)\right\}.$$
(19)
Later we will use the fact that the non-diagonal matrix element of (19) is proportional to the momentum transfer:

$$P^{},S^{}\left|𝔇^\lambda \left\{\overline{\psi }(x)\mathrm{\Gamma }[x,x]\psi (x)\right\}\right|P,S=i(P^{}P)^\lambda P^{},S^{}\left|\overline{\psi }(x)\mathrm{\Gamma }[x,x]\psi (x)\right|P,S.$$
(20)
If we took the fermion mass into account, there would appear an additional mass term in the expression (16) for the axial vector bilocal operator

$$+2im_0^1𝑑uu\overline{\psi }(ux)i\sigma _{\alpha \beta }x^\beta \gamma _5\psi (ux).$$
(21)
We would now like to express the operators (12), (13) through the symmetrical ones:

$`V_{sym}(x,x)`$ $`=`$ $`\overline{\psi }(x)\widehat{x}[x,x]\psi (x),`$ (22)
$`A_{sym}(x,x)`$ $`=`$ $`\overline{\psi }(x)\widehat{x}\gamma _5[x,x]\psi (x).`$ (23)
These operators are “symmetrical” in a sense that on light cone ($`x^2=0`$) they expand in series of symmetrical traceless local operators. In order to do so we should solve the system of ordinary differential equations (15), (16) with respect to $`V_\alpha (x,x)`$ and $`A_\alpha (x,x)`$ (see Appendix of ref.  for details). The vector solution can be presented in the following form:

$$V_\alpha (x,x)=V_\alpha ^{WW}(x,x)+V_\alpha ^{tw3}(x,x),$$
(24)
where

$$V_\alpha ^{WW}(x,x)=\frac{1}{2}_0^1𝑑u\left\{e^{(1u)[(x𝔇)^2x^2𝔇^2]^{1/2}}+e^{(1u)[(x𝔇)^2x^2𝔇^2]^{1/2}}\right\}$$

$$\{[x_\alpha 𝔇^2(x𝔇)𝔇_\alpha ]_0^udvv\overline{\psi }(vx)\widehat{x}[vx,vx]\psi (vx)+\frac{}{x^\alpha }\left[\overline{\psi }(ux)\widehat{x}[ux,ux]\psi (ux)\right]$$

$$i\epsilon _{\alpha \beta \lambda \sigma }x^\beta 𝔇^\lambda \frac{}{x_\sigma }_0^udv\overline{\psi }(vx)\widehat{x}\gamma _5[vx,vx]\psi (vx)\},$$
(25)
$$V_\alpha ^{tw3}(x,x)=\frac{1}{2}_0^1𝑑u\left\{e^{(1u)[(x𝔇)^2x^2𝔇^2]^{1/2}}+e^{(1u)[(x𝔇)^2x^2𝔇^2]^{1/2}}\right\}$$

$$\{i[x^2𝔇_\alpha x_\alpha (x𝔇)]_0^udvv_v^vdt\overline{\psi }(vx)[vx,tx]gG_{\mu \nu }(tx)x^\mu \gamma ^\nu [tx,vx]\psi (vx)$$

$$+𝔊_\alpha (ux)i\epsilon _{\alpha \beta \lambda \sigma }x^\beta 𝔇^\lambda _0^udv𝔊_5^\sigma (vx)\}.$$
(26)
Here $`V_\alpha ^{WW}`$ is the “kinematical” (or Wandzura-Wilczek) contribution to the bilocal vector operator (12). The WW-contribution for the spin-$`\frac{1}{2}`$ target was discussed in details in Refs. . The second term $`V_\alpha ^{tw3}`$ in (24) is the “genuine twist-3” (quark-gluon) contribution to the vector operator. In the case $`x^2=0`$ the formulae (24)-(26) reproduce the result of Belitsky and Müller . It is of necessity to note that the “genuine twist-3” contribution to the scalar product $`x^\alpha V_\alpha (x,x)`$ vanishes, with the first term of the “kinematical” contribution only being survived. In the case when $`x^\alpha =\lambda n^\alpha /2`$ this means that the quark-gluon contribution to the function $`n^\alpha _\alpha `$ vanishes, giving rise to the decomposition (8).

The corresponding expressions for the axial vector operator (13) can be obtained by the following substitution in Eqs. (24)-(26):

$`\mathrm{\Gamma }\psi `$ $``$ $`\mathrm{\Gamma }\gamma _5\psi ,`$
$`𝔊_\alpha `$ $``$ $`𝔊_{5\alpha }.`$
## 3 Diagonal nucleon matrix elements

Let us consider polarized structure functions defined as integrals over diagonal matrix elements of string operators (for a review see  and references therein):

$`{\displaystyle \frac{1}{2M}}{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{d\lambda }{2\pi }}e^{i\lambda x}P,S\left|\overline{\psi }\left({\displaystyle \frac{\lambda n}{2}}\right)\gamma _\mu \gamma _5\psi \left({\displaystyle \frac{\lambda n}{2}}\right)\right|P,S`$ $`=`$ $`g_1(x)(Sn)p_\mu +g_T(x)S_\mu `$ (27)
$`+`$ $`g_3(x)M^2(Sn)n_\mu ,`$
$`{\displaystyle \frac{1}{2M}}{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}{\displaystyle \frac{d\lambda }{2\pi }}e^{i\lambda x}P,S\left|\overline{\psi }\left({\displaystyle \frac{\lambda n}{2}}\right)i\sigma _{\mu \nu }\gamma _5\psi \left({\displaystyle \frac{\lambda n}{2}}\right)\right|P,S`$ $`=`$ $`h_1(x){\displaystyle \frac{1}{M}}(S_\mu p_\nu S_\nu p_\mu )`$ (28)
$`+`$ $`h_L(x)M(Sn)(p_\mu n_\nu p_\nu n_\mu )`$
$`+`$ $`h_3(x)M(S_\mu n_\nu S_\nu n_\mu ),`$
where the spin 4-vector is normalized as $`S^2=1`$. Taking the matrix element of (16) with the mass term (21) we arrive at formulas known as Wandzura-Wilczek relations (). For the region $`0<x<1`$ we have:

$$S_\alpha \left\{g_T(x)_x^1\frac{dy}{y}g_1(y)\frac{m}{M}\left[\frac{h_1(x)}{x}_x^1\frac{dy}{y^2}h_1(y)\right]\right\}=$$

$$=\frac{1}{2M}_x^1\frac{dy}{y}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{d\lambda }{2\pi }e^{i\lambda y}P,S_{}\left|𝔊_{5\alpha }\left(\frac{\lambda n}{2}\right)\right|P,S_{}.$$
(29)
An analogous equation for the region $`1<x<0`$ reads:

$$S_\alpha \left\{g_T(x)+_1^x\frac{dy}{y}g_1(y)\frac{m}{M}\left[\frac{h_1(x)}{x}+_1^x\frac{dy}{y^2}h_1(y)\right]\right\}=$$

$$=\frac{1}{2M}_1^x\frac{dy}{y}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{d\lambda }{2\pi }e^{i\lambda y}P,S_{}\left|𝔊_{5\alpha }\left(\frac{\lambda n}{2}\right)\right|P,S_{}.$$
(30)
It is easy to get the Wanzsura-Wilczek relations for the Mellin moments. For example, in the massless case they have the form:

$$S_\alpha _1^1𝑑xx^m\left[g_T(x)\frac{1}{m+1}g_1(x)\right]=\frac{(i)^m}{2M(m+1)}\frac{^m}{\lambda ^m}P,S_{}\left|𝔊_{5\alpha }\left(\frac{\lambda n}{2}\right)\right|P,S_{}|_{\lambda =0}.$$
(31)
In particular, for $`m=0`$ we have:

$$_1^1𝑑x\left[g_T(x)g_1(x)\right]=0;$$
(32)
for $`m=1`$ we get:

$$_1^1𝑑xx\left[g_T(x)\frac{1}{2}g_1(x)\right]=0;$$
(33)
and for $`m=2`$ we arrive at:

$$S_\alpha _1^1𝑑xx^2\left[g_T(x)\frac{1}{3}g_1(x)\right]=$$

$$=\frac{1}{6M}P,S_{}\left|g\overline{\psi }\left[\stackrel{~}{G}_{\alpha \beta }n^\beta \widehat{n}n_\alpha \stackrel{~}{G}_{\lambda \beta }n^\beta \gamma ^\lambda \right]\psi \right|P,S_{}.$$
(34)
The second term in (34) is purely longitudinal in $`\alpha `$ and does not give any contribution, therefore the transverse part of (34) takes the form

$$S_\alpha _1^1𝑑xx^2\left[g_T(x)\frac{1}{3}g_1(x)\right]=\frac{1}{6M}P,S_{}\left|g\overline{\psi }\stackrel{~}{G}_{\alpha \beta }n^\beta \widehat{n}\psi \right|P,S_{},\alpha =1,2.$$
(35)
The validity of (32)-(34) follows from the observation that $`𝔊_5(\lambda n/2)\lambda ^2`$ and, consequently, the first non-vanishing moment in (31) is the $`x^2`$-moment. The eqs. (32), (33) known as Burkhardt-Cottingham  and Efremov-Leader-Teryaev  sum rules correspondingly. Using the definition of the constant $`d^{(2)}`$, which measures the effect of genuine twist-3 contributions:

$$_1^1𝑑xx^2\left[g_T(x)\frac{1}{3}g_1(x)\right]=\frac{2}{3}d^{(2)},$$
(36)
we come to the identity

$$S_{}^\alpha P,S_{}\left|g\overline{\psi }\stackrel{~}{G}_{\alpha \beta }n^\beta \widehat{n}\psi \right|P,S_{}=4Md^{(2)},$$
(37)
which illustrates that the constant $`d^{(2)}`$ characterizes the correlations of quarks and gluons in a hadron.

## 4 Nondiagonal nucleon matrix elements

In this section we would like to concentrate on gluon part $`V_\alpha ^{tw3}(x,x)`$ of eqs. (24)-(26) with $`x=\frac{1}{2}\lambda n`$. Fourier transformation of its nondiagonal nucleon matrix element

$$_\alpha ^{tw3}(P^{},S^{};P,S;x)=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{d\lambda }{2\pi }e^{i\lambda x}P^{},S^{}\left|V_\alpha ^{tw3}(\frac{\lambda n}{2},\frac{\lambda n}{2})\right|P,S$$
(38)
is the gluon part of the vector GPD, where the spin 4-vectors are normalized as $`S^2=S^2=1`$. Mellin moments of $`_\alpha ^{tw3}`$ are given by

$$_1^1𝑑xx^m_\alpha ^{tw3}(P^{},S^{};P,S;x)=(i)^m\frac{^m}{\lambda ^m}P^{},S^{}\left|V_\alpha ^{tw3}(\frac{\lambda n}{2},\frac{\lambda n}{2})\right|P,S|_{\lambda =0}.$$
(39)
It is not difficult to predict the results of the calculation of the first moments in (39). The $`𝔊`$-term in (26) is proportional to $`x^2`$ (17) and, consequently, to $`\lambda ^2`$. The first and the third terms in (26) are proportional to $`x^3`$ and, consequently, to $`\lambda ^3`$. This means that $`x^0`$\- and $`x^1`$-moments should vanish; the $`𝔊`$-term defines the $`x^2`$-moment, with the other terms of (26) being given their contributions to the highest moments only. Thus one can get by direct calculation that

$$_1^1𝑑x_\alpha ^{tw3}(x;P^{},S^{};P,S)=_1^1𝑑xx_\alpha ^{tw3}(x;P^{},S^{};P,S)=0,$$
(40)
These results mean that the first two moments of the GPDs do not get gluon contributions, the WW approximation turns out to be exact for them.

The first nontrivial Mellin moment sensitive to the genuine twist-3 contributions is the $`x^2`$-moment. It has the following form:

$`{\displaystyle _1^1}𝑑xx^2_\alpha ^{tw3}(x;P^{},S^{};P,S)`$ $`=`$ $`{\displaystyle \frac{1}{3}}P^{},S^{}|g\overline{\psi }[\stackrel{~}{G}_{\alpha \beta }n^\beta \widehat{n}`$ (41)
$`n_\alpha \stackrel{~}{G}_{\lambda \beta }n^\beta \gamma ^\lambda ]\gamma _5\psi |P,S.`$
The transverse part of eq. (41) has the form

$$_1^1𝑑xx^2_\alpha ^{tw3}(x;P^{},S^{};P,S)=\frac{1}{3}P^{},S^{}\left|g\overline{\psi }\stackrel{~}{G}_{\alpha \beta }n^\beta \widehat{n}\gamma _5\psi \right|P,S,\alpha =1,2.$$
(42)
The expressions for the axial moments can be obtained using the substitution: $`\mathrm{\Gamma }\psi \mathrm{\Gamma }\gamma _5\psi `$.

## 5 Estimate of the nucleon matrix elements

In order to estimate the magnitude of the second moments (35), (42), we should analyze the nondiagonal hadron matrix elements of the quark-gluon operators:

$`O_{}`$ $`=`$ $`g\overline{\psi }\stackrel{~}{G}_\lambda n^\lambda \widehat{n}\psi ,`$ (43)
$`O_5`$ $`=`$ $`g\overline{\psi }\stackrel{~}{G}_\lambda n^\lambda \widehat{n}\gamma _5\psi `$ (44)
(in this section we will use a notation: $`O_{}O_\alpha ,\alpha =1,2`$). For this purpose we first consider the nondiagonal quark matrix elements of the operators (43), (44). A method of calculation of such matrix elements in the model of the instanton vacuum was suggested in Refs. . It was developed in the framework of the Effective Chiral Theory of the nucleon  (see also recent reviews  and the references therein). The essence of the method is in the following.

In the model of the instanton vacuum  the fermion partition function in the presence of $`N_+`$ instantons and $`N_{}`$ anti-instantons can be represented as a functional integral over the fermion fields with a fermion effective action $`S_{\text{eff}}[\overline{\psi },\psi ]`$:

$$Z_{N_\pm }^{\text{fermions}}=Z_{N_\pm }C𝒟\overline{\psi }𝒟\psi \mathrm{exp}\left\{iS_{\text{eff}}[\overline{\psi },\psi ]\right\},$$
(45)
where in the case of one quark flavor

$$S_{\text{eff}}[\overline{\psi },\psi ]=\frac{d^4k}{(2\pi )^4}\overline{\psi }(k)\left[\widehat{k}M(ik)\right]\psi (k)$$
(46)
and $`Z_{N_\pm }`$ is the variational partition function of the gluodynamics. Eq. (46) points out that a quark propagating in the instanton medium acquires a dynamical momentum-dependent mass $`M(ik)`$. It is this mass that leads to the nonvanishing value of the vacuum condensate $`\overline{\psi }\psi `$ and, consequently, to the spontaneous chiral symmetry breaking in QCD . The quark propagator in the instanton vacuum model can be written in the momentum representation as

$$\widehat{S}(k)=i\frac{\widehat{k}+M(ik)}{k^2M^2(ik)+i\epsilon }.$$
(47)
The dynamical mass $`M(ik)`$ has the form:

$$M(ik)=M|F(ik)|^2.$$
(48)
Here the mass $`M\rho /R^2345`$ MeV, where $`\rho `$ is the average size of the instantons in the instanton medium and $`R`$ is the average distance between instantons. The instanton medium is characterized by its “packing fraction” $`\rho /R`$ with a phenomenological value $`\rho /R1/3`$. Thus, the product $`M\rho \rho ^2/R^2`$ is parametrically small. The function $`F(k)`$ in (48) is a form factor proportional to the Fourier transformation of the wave function of the fermion zero-mode:

$`F(k)`$ $`=`$ $`t{\displaystyle \frac{d}{dt}}\left[I_0(t)K_0(t)I_1(t)K_1(t)\right],t={\displaystyle \frac{1}{2}}k\rho ,`$ (49)
$`F(k)`$ $`=`$ $`1,\rho 0,`$
where $`I_n(t)`$ are the modified Bessel functions and $`K_n(t)`$ are the modified Hankel functions.

The effective action (46) can be used in the calculation of the quark matrix elements of the quark-gluon operators (43)-(44), the operators being replaced by the effective ones:

$$p^{},s^{}\left|O\right|p,s_{\text{inst.vac.}}=p^{},s^{}\left|\mathrm{`}\mathrm{`}O\mathrm{"}\right|p,s_{\text{eff}}.$$
(50)
Here $`\mathrm{`}\mathrm{`}O\mathrm{"}`$ is the fermion operator, which represents the quark-gluon operator $`O`$ in the effective fermion theory. The explicit form of $`\mathrm{`}\mathrm{`}O\mathrm{"}`$ was obtained in . The effective action and effective operators being known, it is possible to calculate the quark matrix elements of interest.

The direct calculation for the case of the operators (43), (44) results in the following expressions:

$`p^{},s^{}\left|\mathrm{`}\mathrm{`}O_{}\mathrm{"}\right|p,s_{\text{eff}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}M\rho ^2{\displaystyle \frac{d^4k}{(2\pi )^4}G(ik)}`$ (51)
$`\overline{u}(p^{},s^{})[F^{}(i(kp))F(ip)\widehat{n}\gamma _5\widehat{S}(kp)\widehat{\mathrm{\Sigma }}_{}(k,n)`$
$`F^{}(ip^{})F(i(kp^{}))\widehat{\mathrm{\Sigma }}_{}(k,n)\widehat{S}(kp^{})\widehat{n}\gamma _5]u(p,s)`$
$`+(\text{corrections}),`$
$`p^{},s^{}\left|\mathrm{`}\mathrm{`}O_5\mathrm{"}\right|p,s_{\text{eff}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}M\rho ^2{\displaystyle \frac{d^4k}{(2\pi )^4}G(ik)}`$ (52)
$`\overline{u}(p^{},s^{})[F^{}(i(kp))F(ip)\widehat{n}\widehat{S}(kp)\widehat{\mathrm{\Sigma }}_{}(k,n)`$
$`F^{}(ip^{})F(i(kp^{}))\widehat{\mathrm{\Sigma }}_{}(k,n)\widehat{S}(kp^{})\widehat{n}]u(p,s)`$
$`+(\text{corrections}).`$
Here $`G(k)`$ is the form factor proportional to the Fourier transformation of the one-instanton dual field strength

$`G(k)`$ $`=`$ $`32\pi ^2\left[\left({\displaystyle \frac{1}{2}}+{\displaystyle \frac{4}{t^2}}\right)K_0(t)+\left({\displaystyle \frac{2}{t}}+{\displaystyle \frac{8}{t^3}}\right)K_1(t){\displaystyle \frac{8}{t^4}}\right],t=k\rho ,`$ (53)
$`G(k)`$ $`=`$ $`4\pi ^2,\rho 0,`$ (54)
and the following definitions have been used:

$$\widehat{\mathrm{\Sigma }}_\alpha (k,n)=\frac{k_\alpha k^\mu }{k^2}n^\nu \sigma _{\mu \nu }+\frac{k^\nu k^\lambda }{k^2}n_\lambda \sigma _{\alpha \nu }\frac{1}{2}n^\nu \sigma _{\alpha \nu },$$

$$\sigma _{\mu \nu }=\frac{i}{2}[\gamma _\mu ,\gamma _\nu ],n^\lambda =\frac{1}{\overline{p}_+\sqrt{2}}(1;0,0,1),$$

$`u(p,s)`$ \- is the quark Dirac bispinor.

Using equations of motion, the expressions (51), (52) can be reduced to the form:

$`p^{},s^{}\left|\mathrm{`}\mathrm{`}O_{}\mathrm{"}\right|p,s_{\text{eff}}`$ $`=`$ $`MI(p)\overline{u}(p^{},s^{})\widehat{n}\gamma _{}\gamma _5u(p,s),`$ (55)
$`p^{},s^{}\left|\mathrm{`}\mathrm{`}O_5\mathrm{"}\right|p,s_{\text{eff}}`$ $`=`$ $`\xi MI(p)\overline{u}(p^{},s^{})\widehat{n}\gamma _{}u(p,s),`$ (56)
where $`\xi =\frac{1}{2}(p^{}p)n`$, and the integral $`I(p)`$ is defined in leading oder of $`\rho ^2/R^2`$ as

$$I(p)=\rho ^2\frac{1}{3}\frac{d^4k}{(2\pi )^4}\frac{G(k)F^3(kp)}{(kp)^2+M^2F^2(kp)}\left[1+4\frac{(kp)^2}{k^2p^2}\right],$$
(57)
where the integration is performed in Euclidean space with $`p^2=M^2`$. In the case of diagonal matrix elements we reproduce (up to a sign) the results of Ref. :

$`p,s\left|\mathrm{`}\mathrm{`}O_{}\mathrm{"}\right|p,s_{\text{eff}}`$ $`=`$ $`4Md_{\text{quark}}^{(2)}s_{},`$ (58)
$`p,s\left|\mathrm{`}\mathrm{`}O_5\mathrm{"}\right|p,s_{\text{eff}}`$ $`=`$ $`0,`$ (59)
where

$$d_{\text{quark}}^{(2)}=\frac{1}{2}I(p)$$
(60)
is a quark constant analogous to the hadron constant (36) in Eq. (37). Thus, the nondiagonal quark matrix elements (55), (56) of the twist-3 quark-gluon operators can be expressed through the quark constant:

$`p^{},s^{}\left|\mathrm{`}\mathrm{`}O_{}\mathrm{"}\right|p,s_{\text{eff}}`$ $`=`$ $`2Md_{\text{quark}}^{(2)}\overline{u}(p^{},s^{})\widehat{n}\gamma _{}\gamma _5u(p,s),`$ (61)
$`p^{},s^{}\left|\mathrm{`}\mathrm{`}O_5\mathrm{"}\right|p,s_{\text{eff}}`$ $`=`$ $`2\xi Md_{\text{quark}}^{(2)}\overline{u}(p^{},s^{})\widehat{n}\gamma _{}u(p,s).`$ (62)
In the limit of a small average size of instantons $`\rho `$ we get the following estimate of the integral $`I(p)`$ (57):

$$d_{\text{quark}}^{(2)}I(p)\rho ^2M^2\mathrm{ln}\rho ^2M^2,\rho R.$$
(63)
The quantity $`M\rho \pi \rho ^2/R^2`$ is a small parameter in the effective fermion theory and, thus, we notice the quark matrix elements of the quark-gluon operators to be parametrically suppressed as a forth power of the “packing fraction” $`\rho /R`$ of the instanton vacuum. The numerical calculation  of $`I(p)`$ gives $`d_{\text{quark}}^{(2)}=0.011`$ at the phenomenological value $`M\rho =0.58`$.

The non-zero value of $`d_{\text{quark}}^{(2)}`$ appears only at the order $`\pi ^2\rho ^4/R^4`$. At the present stage of development of our theory we cannot make precise predictions at this order, as this would require to go beyond the one-instanton approximation. Therefore we consider the value of $`d_{\text{quark}}^{(2)}=0.011`$ only as an order of magnitude estimate. Using this order of magnitude estimate we can predict the order of magnitude of $`d^{(2)}`$ for the nucleon :

$$d^{(2)}10^3,\mathrm{with}d_p^{(2)}d_n^{(2)}1/N_c,d_p^{(2)}+d_n^{(2)}1/N_c^2.$$
(64)
Without the suppression due to the instanton packing fraction a “natural” value of the nucleon $`d^{(2)}`$ would be of order $`_0^1𝑑xx^2g_1(x)\mathrm{few}\times 10^2`$. A recent experiment at SLAC  quotes the values $`d_p^{(2)}=(3.2\pm 1.7)10^3`$ and $`d_n^{(2)}=(7.9\pm 4.8)10^3`$ at $`Q^25`$ GeV<sup>2</sup>. These data hint on possible “anomalous suppression” of $`d^{(2)}`$ as predicted in the model of the instanton vacuum .

Using the results (61) and (62) for the quark matrix elements we can make a simple estimate of the corresponding nucleon matrix element at $`\mathrm{\Delta }^2=0`$. In such a kind of estimate we assume that the nucleon matrix element at $`\mathrm{\Delta }^2=0`$ has the same structure as in Eqs. (61,62). Such an assumption corresponds to the use of a simple quark model to “translate” the quark matrix elements to the nucleon ones. More sophisticated estimates can be done with help of the chiral quark-soliton model, see e.g. Refs. . Note, however, that for the very fact of the suppression of the twist-3 nucleon matrix elements the usage of the simple quark model is sufficient. In the simple quark model of the nucleon we arrive at:

$`P^{},S^{}\left|\mathrm{`}\mathrm{`}O_{}\mathrm{"}\right|P,S_{\text{eff}}`$ $`=`$ $`2M_Nd^{(2)}\overline{N}(P^{},S^{})\widehat{n}\gamma _{}\gamma _5N(P,S),`$ (65)
$`P^{},S^{}\left|\mathrm{`}\mathrm{`}O_5\mathrm{"}\right|P,S_{\text{eff}}`$ $`=`$ $`2\xi M_Nd^{(2)}\overline{N}(P^{},S^{})\widehat{n}\gamma _{}N(P,S),`$ (66)
where $`N(P,S)`$ is the nucleon Dirac bispinor. We come to the conclusion that the hadron matrix elements of the “genuine twist-3” (quark-gluon) opeators are parametrically suppressed in the effective chiral theory by the “packing fraction” of the instanton vacuum.

## 6 GPDs decomposition

The quark twist-2 (WW) contributions to the nucleon GPD $`_{}`$ and $`\stackrel{~}{}_{}`$ have been obtained in Refs. . In this section we would like to compare the WW contributions with the genuine twist-3 ones, calculated in the previous section (to be precise, we will compare their $`x^2`$-moments). Following Ref.  we present the “transverse” part of the GPDs in the form<sup>1</sup><sup>1</sup>1We slightly changed the notations of Ref.  and added supplementary functions which were missed in Ref. , see also .:

$`_\mu (x,\xi ,\mathrm{\Delta })`$ $`=`$ $`\overline{N}(P^{},S^{})\{(H+E)\gamma _\mu ^{}`$ (67)
$`+`$ $`G_1{\displaystyle \frac{\mathrm{\Delta }_\mu ^{}}{2M}}+G_2\gamma _\mu ^{}+G_3\mathrm{\Delta }_\mu ^{}\widehat{n}+G_4i\epsilon _{\mu \nu }^{}\mathrm{\Delta }_{}^\nu \widehat{n}\gamma _5\}N(P,S),`$
$`\stackrel{~}{}_\mu (x,\xi ,\mathrm{\Delta })`$ $`=`$ $`\overline{N}(P^{},S^{})\{\stackrel{~}{H}\gamma _\mu ^{}\gamma _5+\stackrel{~}{E}{\displaystyle \frac{\mathrm{\Delta }_\mu ^{}}{2M}}\gamma _5`$ (68)
$`+`$ $`\stackrel{~}{G}_1{\displaystyle \frac{\mathrm{\Delta }_\mu ^{}}{2M}}\gamma _5+\stackrel{~}{G}_2\gamma _\mu ^{}\gamma _5+\stackrel{~}{G}_3\mathrm{\Delta }_\mu ^{}\widehat{n}\gamma _5+\stackrel{~}{G}_4i\epsilon _{\mu \nu }^{}\mathrm{\Delta }_{}^\nu \widehat{n}\}N(P,S).`$
Note that in the above equations all functions depend on $`x,\xi `$ and $`\mathrm{\Delta }^2`$. The first three moments of the quark contributions have the following structure.

The $`x^0`$-moments are:

$`{\displaystyle _1^1}𝑑x_\mu (x,\xi ,\mathrm{\Delta })`$ $`=`$ $`\overline{N}(P^{},S^{})\left\{G_M(\mathrm{\Delta }^2)\left[\gamma _\mu ^{}\right]\right\}N(P,S),`$ (69)
$`{\displaystyle _1^1}𝑑x\stackrel{~}{}_\mu (x,\xi ,\mathrm{\Delta })`$ $`=`$ $`\overline{N}(P^{},S^{})\left\{G_P(\mathrm{\Delta }^2)\left[{\displaystyle \frac{\mathrm{\Delta }_\mu ^{}}{2M}}\gamma _5\right]+G_A(\mathrm{\Delta }^2)\left[\gamma _\mu ^{}\gamma _5\right]\right\}N(P,S).`$
Here the well-known sum rules  have been used:

$`{\displaystyle _1^1}𝑑xH(x,\xi ,\mathrm{\Delta }^2)`$ $`=`$ $`F_1(\mathrm{\Delta }^2),`$
$`{\displaystyle _1^1}𝑑xE(x,\xi ,\mathrm{\Delta }^2)`$ $`=`$ $`F_2(\mathrm{\Delta }^2),`$
$`{\displaystyle _1^1}𝑑x\stackrel{~}{H}(x,\xi ,\mathrm{\Delta }^2)`$ $`=`$ $`G_A(\mathrm{\Delta }^2),`$
$`{\displaystyle _1^1}𝑑x\stackrel{~}{E}(x,\xi ,\mathrm{\Delta }^2)`$ $`=`$ $`G_P(\mathrm{\Delta }^2),`$ (71)
where $`F_1(\mathrm{\Delta }^2)`$ and $`F_2(\mathrm{\Delta }^2)`$ are the Dirac and Pauli form factors, $`G_P(\mathrm{\Delta }^2)`$ and $`G_A(\mathrm{\Delta }^2)`$ are the pseudo-scalar and axial-vector form factors; and, by definition, the magnetic form factor is

$$G_M(\mathrm{\Delta }^2)=F_1(\mathrm{\Delta }^2)+F_2(\mathrm{\Delta }^2).$$
(72)
The results (69) and (LABEL:A0mom) imply that

$$_1^1𝑑xG_i(x,\xi ,\mathrm{\Delta })=0,_1^1𝑑x\stackrel{~}{G}_i(x,\xi ,\mathrm{\Delta })=0.$$

These sum rules can be considered as a non-forward generalization of Burkhardt-Cottingham  sum rules (32).

For the $`x`$-moments of the vector GPD we get:

$`{\displaystyle _1^1}𝑑xxG_1(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{}{\xi }}{\displaystyle _1^1}𝑑xxE(x,\xi ,\mathrm{\Delta }),`$ (73)
$`{\displaystyle _1^1}𝑑xxG_2(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[G_A(\mathrm{\Delta }^2){\displaystyle _1^1}𝑑xx(H+E)(x,\xi ,\mathrm{\Delta })\right],`$ (74)
$`{\displaystyle _1^1}𝑑xxG_3(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`0,`$ (75)
$`{\displaystyle _1^1}𝑑xxG_4(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`0.`$ (76)
Note that the $`x`$-moment of the GPD $`G_2`$ does not vanish in the forward limit and gives the quark orbital momentum of the nucleon :

$$\underset{\mathrm{\Delta }0}{lim}_1^1𝑑xxG_2(x,\xi ,\mathrm{\Delta })=J_q+\frac{1}{2}\mathrm{\Delta }q=L_q.$$

For the $`x`$-moments of the axial-vector GPD we have:

$`{\displaystyle _1^1}𝑑xx\stackrel{~}{G}_1(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[F_2(\mathrm{\Delta }^2)+\left(\xi {\displaystyle \frac{}{\xi }}1\right){\displaystyle _1^1}𝑑xx\stackrel{~}{E}(x,\xi ,\mathrm{\Delta })\right],`$ (77)
$`{\displaystyle _1^1}𝑑xx\stackrel{~}{G}_2(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[\xi ^2G_E(\mathrm{\Delta }^2){\displaystyle \frac{\mathrm{\Delta }^2}{4M^2}}F_2(\mathrm{\Delta }^2){\displaystyle _1^1}𝑑xx\stackrel{~}{H}(x,\xi ,\mathrm{\Delta })\right],`$ (78)
$`{\displaystyle _1^1}𝑑xx\stackrel{~}{G}_3(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{1}{4}}\xi G_E(\mathrm{\Delta }^2),`$ (79)
$`{\displaystyle _1^1}𝑑xx\stackrel{~}{G}_4(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{1}{4}}G_E(\mathrm{\Delta }^2),`$ (80)
where

$$G_E(\mathrm{\Delta }^2)=F_1(\mathrm{\Delta }^2)+\frac{\mathrm{\Delta }^2}{4M^2}F_2(\mathrm{\Delta }^2)$$
(81)
is the electric form factor. These moments have no gluon contribution due to (40) and they have already been discussed in  and . Similarly to the previous case one can consider the relations (73)-(76) and (77)-(80) as a non-forward generalization of Efremov-Leader-Teryaev  sum rule (33).

We are passing now to a consideration of the $`x^2`$-moments. For these moments the “genuine twist-3” (quark-gluon) contribution (42) does not vanish and one can decompose the GPDs $`G`$ and $`\stackrel{~}{G}`$ from (67), (68) into the WW and genuine twist-3 parts:

$$G_i^{(2)}(\xi ,\mathrm{\Delta })=G_i^{(2)WW}(\xi ,\mathrm{\Delta })+G_i^{(2)tw3}(\xi ,\mathrm{\Delta }),i=1,2,3,4.$$
(82)
For the $`x^2`$-moments of the vector GPD “kinematical” part we get:

$`{\displaystyle _1^1}𝑑xx^2G_1^{WW}(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`\left({\displaystyle \frac{1}{3}}\right)\left[\xi F_2(\mathrm{\Delta }^2)+\xi \left(\xi {\displaystyle \frac{}{\xi }}1\right)\stackrel{~}{E}^{(1)}(\xi ,\mathrm{\Delta }){\displaystyle \frac{}{\xi }}E^{(2)}(\xi ,\mathrm{\Delta })\right],`$
$`{\displaystyle _1^1}𝑑xx^2G_2^{WW}(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{1}{3}}[\xi ^2G_M(\mathrm{\Delta }^2)+{\displaystyle \frac{\mathrm{\Delta }_{}^2}{4M^2}}(\xi {\displaystyle \frac{}{\xi }}1)\stackrel{~}{E}^{(1)}(\xi ,\mathrm{\Delta })`$ (84)
$`+\stackrel{~}{H}^{(1)}(\xi ,\mathrm{\Delta })2(H+E)^{(2)}(\xi ,\mathrm{\Delta })],`$
$`{\displaystyle _1^1}𝑑xx^2G_3^{WW}(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{1}{6}}[\xi G_M(\mathrm{\Delta }^2)+(1{\displaystyle \frac{\mathrm{\Delta }^2}{4M^2}})\xi (\xi {\displaystyle \frac{}{\xi }}1)\stackrel{~}{E}^{(1)}(\xi ,\mathrm{\Delta })`$ (85)
$`{\displaystyle \frac{}{\xi }}(H+E)^{(2)}(\xi ,\mathrm{\Delta })],`$
$`{\displaystyle _1^1}𝑑xx^2G_4^{WW}(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`\left({\displaystyle \frac{1}{6}}\right)\left(\xi {\displaystyle \frac{}{\xi }}1\right)\left[\stackrel{~}{H}^{(1)}(\xi ,\mathrm{\Delta })+{\displaystyle \frac{\mathrm{\Delta }^2}{4M^2}}\stackrel{~}{E}^{(1)}(\xi ,\mathrm{\Delta })\right],`$ (86)
where we used shorthand notations for the Mellin moments of the twist-2 vector GPDs:

$`H^{(m)}(\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle _1^1}𝑑xx^mH(x,\xi ,\mathrm{\Delta }),`$
$`E^{(m)}(\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle _1^1}𝑑xx^mE(x,\xi ,\mathrm{\Delta });`$
and analogous notations for the twist-2 pseudovector GPDs $`\stackrel{~}{H}`$ and $`\stackrel{~}{E}`$.

For the $`x^2`$-moments of the axial-vector GPD we get:

$`{\displaystyle _1^1}𝑑xx^2\stackrel{~}{G}_1^{WW}(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`\left({\displaystyle \frac{1}{3}}\right)\left(\xi {\displaystyle \frac{}{\xi }}2\right)\left[E^{(1)}(\xi ,\mathrm{\Delta })\stackrel{~}{E}^{(2)}(\xi ,\mathrm{\Delta })\right],`$ (87)
$`{\displaystyle _1^1}𝑑xx^2\stackrel{~}{G}_2^{WW}(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{1}{3}}[\xi ^2(G_A(\mathrm{\Delta }^2)+(H+E)^{(1)}(\xi ,\mathrm{\Delta }))`$ (88)
$`+{\displaystyle \frac{\mathrm{\Delta }_{}^2}{4M^2}}(\xi {\displaystyle \frac{}{\xi }}2)E^{(1)}(\xi ,\mathrm{\Delta })2\stackrel{~}{H}^{(2)}(\xi ,\mathrm{\Delta })],`$
$`{\displaystyle _1^1}𝑑xx^2\stackrel{~}{G}_3^{WW}(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`{\displaystyle \frac{1}{6}}[\xi (G_A(\mathrm{\Delta }^2)+(H+E)^{(1)}(\xi ,\mathrm{\Delta }))`$ (89)
$`+`$ $`(1{\displaystyle \frac{\mathrm{\Delta }^2}{4M^2}})\xi (\xi {\displaystyle \frac{}{\xi }}2)E^{(1)}(\xi ,\mathrm{\Delta }){\displaystyle \frac{}{\xi }}\stackrel{~}{H}^{(2)}(\xi ,\mathrm{\Delta })],`$
$`{\displaystyle _1^1}𝑑xx^2\stackrel{~}{G}_4^{WW}(x,\xi ,\mathrm{\Delta })`$ $`=`$ $`\left({\displaystyle \frac{1}{6}}\right)\left(\xi {\displaystyle \frac{}{\xi }}2\right)\left[H^{(1)}(\xi ,\mathrm{\Delta })+{\displaystyle \frac{\mathrm{\Delta }^2}{4M^2}}E^{(1)}(\xi ,\mathrm{\Delta })\right].`$ (90)
Let us consider the quark-gluon contributions to the moments of $`_{}`$ and $`\stackrel{~}{}_{}`$. We have already shown that their $`x^0`$\- and $`x^1`$-moments vanish (40), with the $`x^2`$-moment being given by (42) in the vector case. Substituting (65) into (42) and decomposing the result according to (67) we get the following expressions for the $`x^2`$-moments of the genuine twist-3 vector GPD at $`\mathrm{\Delta }^2=0`$:

$`{\displaystyle _1^1}𝑑xx^2G_1^{tw3}(x,\xi )`$ $`=`$ $`0,`$ (91)
$`{\displaystyle _1^1}𝑑xx^2G_2^{tw3}(x,\xi )`$ $`=`$ $`{\displaystyle \frac{2}{3}}(1\xi ^2)d^{(2)},`$ (92)
$`{\displaystyle _1^1}𝑑xx^2G_3^{tw3}(x,\xi )`$ $`=`$ $`{\displaystyle \frac{1}{3}}\xi d^{(2)},`$ (93)
$`{\displaystyle _1^1}𝑑xx^2G_4^{tw3}(x,\xi )`$ $`=`$ $`{\displaystyle \frac{1}{3}}d^{(2)}.`$ (94)
Here $`d^{(2)}`$ is given by Eq. (36). For the axial-vector GPD we get at $`\mathrm{\Delta }^2=0`$:

$`{\displaystyle _1^1}𝑑xx^2\stackrel{~}{G}_1^{tw3}(x,\xi )`$ $`=`$ $`0,`$ (95)
$`{\displaystyle _1^1}𝑑xx^2\stackrel{~}{G}_2^{tw3}(x,\xi )`$ $`=`$ $`{\displaystyle \frac{2}{3}}(1\xi ^2)d^{(2)},`$ (96)
$`{\displaystyle _1^1}𝑑xx^2\stackrel{~}{G}_3^{tw3}(x,\xi )`$ $`=`$ $`{\displaystyle \frac{1}{3}}\xi d^{(2)},`$ (97)
$`{\displaystyle _1^1}𝑑xx^2\stackrel{~}{G}_4^{tw3}(x,\xi )`$ $`=`$ $`{\displaystyle \frac{1}{3}}d^{(2)}.`$ (98)
The expressions (91)-(94) and (95)-(98) for the the genuine twist-3 contributions should be compared to the expressions (6)-(86) and (87)-(90) of the WW part for the quark contributions correspondingly. As we saw above, in the instanton vacuum the constant $`d^{(2)}`$ is strongly suppressed by the packing fraction of instantons . Here we expressed the $`x^2`$-moments of the genuine twist-3 GPDs in terms of the constant $`d^{(2)}`$, see expressions (91)-(94) and (95)-(98). We conclude that the twist-3 vector and axial-vector GPDs are parametrically suppressed relative to the WW ones due to the packing fraction of the instanton vacuum:

$$\frac{G_i^{(2)tw3}}{G_i^{(2)WW}}\pi ^2\left(\frac{\rho ^2}{R^2}\right)^2\mathrm{ln}\left(\frac{\rho ^2}{R^2}\right)1.$$
(99)
This indicates that in the QCD vacuum the “genuine twist-3” (quark-gluon) contributions turn out to be small relative to the purely quark ones.

## Acknowledgments

We acknowledge useful discussions with N. Kivel, A. Mukherjee and M. Vanderhaeghen. We thank Vadim Guzey for careful reading of the manuscript. This work is supported by the Sofja Kovalevskaja Programme of the Alexander von Humboldt Foundation, the Federal Ministry of Education and Research and the Programme for Investment in the Future of German Government.