text:
# Comment about quasi-isotropic solution of Einstein equations near cosmological singularity

## Abstract

We generalize for the case of arbitrary hydrodynamical matter the quasi-isotropic solution of Einstein equations near cosmological singularity, found by Lifshitz and Khalatnikov in 1960 for the case of the radiation-dominated Universe. It is shown that this solution always exists, but dependence of terms in the quasi-isotropic expansion acquires a more complicated form.

1L.D. Landau Institute for Theoretical Physics of Russian Academy of Sciences, Kosygin str. 2, 117334, Moscow, Russia
2Landau Network - Centro Volta, Villa Olmo, via Cantoni 1, 22100 Como, Italy
3Tel Aviv University, Tel Aviv University, Raymond and Sackler Faculty of Exact Sciences, School of Physics and Astronomy, Ramat Aviv, 69978, Israel

PACS: 04.20-q; 04.20.Dw

In the paper  by Lifshitz and Khalatnikov, the quasi-isotropic solution of the Einstein equations near a cosmological singularity was found provided the Universe was filled by radiation with the equation of state $`p=\frac{\epsilon }{3}`$. The metric of this solution was written down in the synchronous system of reference

$$ds^2=dt^2\gamma _{\alpha \beta }dx^\alpha dx^\beta ,$$
(1)
where spatial metric $`\gamma _{\alpha \beta }`$ near the singularity has the form

$$\gamma _{\alpha \beta }=ta_{\alpha \beta }+t^2b_{\alpha \beta }+\mathrm{}$$
(2)
where $`a_{\alpha \beta }`$ and $`b_{\alpha \beta }`$ are functions of spatial coordinates. The functions $`a_{\alpha \beta }`$ are chosen arbitrary, and then the functions $`b_{\alpha \beta }`$ and also the energy and velocity distributions for matter can be expressed through these functions (for details see  and also ).

In correspondence with the standard cosmological model of the hot Universe, it was supposed that the natural equation of state for the matter near the cosmological singularity is that of radiation: $`p=\epsilon /3`$. However, nowadays the situation has changed in connection with the development of inflationary cosmological models, which as an important ingredient contain inflaton scalar field or/and other exotic types of matter . One can add also that the appearance of brane and M theory cosmological models  and the discovery of the cosmic acceleration  suggests that the matter playing essential role on different stages of cosmological evolution can obey very different equations of state . Thus, generalization of the old quasi-isotropic solution of the Einstein equations near the cosmological singularity can be useful in this new context.

In this note we make such a generalization for the equation of state:

$$p=k\epsilon ,$$
(3)
where $`p`$ denotes pressure and $`\epsilon `$ denotes energy density 111Actually, the authors have known this generalization for a long time but have never published it in detail in regular journals.. The Friedmann isotropic solution near the singularity for such a matter behaves as

$$aa_0t^m,$$
(4)
where

$$m=\frac{4}{3(1+k)}.$$
(5)
We look for an expression for a spatial metric in the following form:

$$\gamma _{\alpha \beta }=t^ma_{\alpha \beta }+t^nb_{\alpha \beta },$$
(6)
where the power index $`m`$ is given by Eq. (5). We leave the power index $`n`$ free for some time, requiring only that

$$n>m.$$
(7)
The inverse metric reads

$$\gamma ^{\alpha \beta }=\frac{a^{\alpha \beta }}{t^m}\frac{b^{\alpha \beta }}{t^{2mn}},$$
(8)
where $`a^{\alpha \beta }`$ is defined by the relation

$$a^{\alpha \beta }a_{\beta \gamma }=\delta _\gamma ^\alpha $$
(9)
while the indices of all the other matrices are lowered and raised by $`a_{\alpha \beta }`$ and $`a^{\alpha \beta }`$, for example,

$$b_\beta ^\alpha =a^{\alpha \gamma }b_{\gamma \beta }.$$
(10)
Let us write down also expressions for the extrinsic curvature, its contractions and its derivatives:

$$\kappa _{\alpha \beta }\frac{\gamma _{\alpha \beta }}{t}=mt^{m1}a_{\alpha \beta }+nt^{n1}b_{\alpha \beta },$$
(11)
$$\kappa _\alpha ^\beta =\frac{m\delta _\alpha ^\beta }{t}+\frac{(nm)b_\alpha ^\beta }{t^{mn+1}},$$
(12)
$$\kappa _\alpha ^\alpha =\frac{3m}{t}+\frac{(nm)b}{t^{mn+1}},$$
(13)
$$\frac{\kappa _\alpha ^\beta }{t}=\frac{m\delta _\alpha ^\beta }{t^2}\frac{(mn+1)(nm)b_\alpha ^\beta }{t^{mn+2}},$$
(14)
$$\frac{\kappa _\alpha ^\alpha }{t}=\frac{3m}{t^2}\frac{(mn+1)(nm)b}{t^{mn+2}},$$
(15)
$$\kappa _\alpha ^\beta \kappa _\beta ^\alpha =\frac{3m^2}{t^2}+\frac{2m(nm)b}{t^{mn+2}}.$$
(16)
We need also an explicit expression for the determinant of the spatial metric:

$$\gamma det\gamma _{\alpha \beta }=t^{3m}(1+t^{nm}b)deta,$$
(17)
$$\dot{\gamma }\frac{\gamma }{t}=(3mt^{3m1}+b(2m+n)t^{2m+n1})deta,$$
(18)
$$\frac{\dot{\gamma }}{\gamma }=\frac{3m}{t}\left(1+\frac{b(nm)t^{nm}}{3m}\right).$$
(19)
Now, using well-known expressions for the components of the Ricci tensor :

$$R_0^0=\frac{1}{2}\frac{\kappa _\alpha ^\alpha }{t}\frac{1}{4}\kappa _\alpha ^\beta \kappa _\beta ^\alpha ,$$
(20)
$$R_\alpha ^0=\frac{1}{2}(\kappa _{\alpha ;\beta }^\beta \kappa _{\beta ;\alpha }^\beta ),$$
(21)
$$R_\alpha ^\beta =P_\alpha ^\beta \frac{1}{2}\frac{\kappa _\alpha ^\beta }{t}\frac{\dot{\gamma }}{4\gamma }\kappa _\alpha ^\beta ,$$
(22)
where $`P_\alpha ^\beta `$ is a three-dimensional part of the Ricci tensor, and substituting into Eqs. (20)-(22) the expressions (11)-(19), one get

$$R_0^0=\frac{3m(2m)}{4t^2}\frac{(n1)(nm)b}{2t^{mn+2}},$$
(23)
$$R_\alpha ^0=\frac{nm}{2t^{mn+1}}(b_{\alpha ;\beta }^\beta b_{;\alpha }),$$
(24)
$`R_\alpha ^\beta ={\displaystyle \frac{\stackrel{~}{P}_\alpha ^\beta }{t^m}}+{\displaystyle \frac{m(23m)\delta _\alpha ^\beta }{4t^2}}`$
$`+{\displaystyle \frac{(nm)(22nm)b_\alpha ^\beta }{4t^{mn+2}}}{\displaystyle \frac{m(nm)b\delta _\alpha ^\beta }{4t^{mn+2}}}.`$ (25)
Notice, that in Eq. (25) $`\stackrel{~}{P}_\alpha ^\beta `$ denotes a three-dimensional Ricci tensor constructed by using the metrics $`a_{\alpha \beta }`$. The terms in the curvature tensor $`P_\alpha ^\beta `$, which are proportional to $`\beta _{\alpha \beta }`$ have the time dependence $`\frac{1}{t^{2mn}}`$ and are less divergent than the first term in the right-hand side of Eq. (25) provided the condition (7) is satisfied.

Now, let us write down the expressions for the components of the energy-momentum tensor of the perfect fluid

$$T_{ik}=(\epsilon +p)u_iu_kpg_{ik},$$
(26)
satisfying the equation of state (3). Up to higher-order corrections, they have the following form :

$$T_0^0=\epsilon $$
(27)
$$T_\alpha ^0=\epsilon (k+1)u_\alpha ,$$
(28)
$$T_\alpha ^\beta =k\epsilon \delta _\alpha ^\beta ,$$
(29)
$$T=T_i^i=\epsilon (13k).$$
(30)
Using the Einstein equations

$$R_i^j=8\pi G(T_i^j\frac{1}{2}\delta _i^jT),$$
(31)
one has from $`00`$-component of these equations:

$$8\pi G\epsilon =\frac{1}{3k+1}\left(\frac{3m(2m)}{2t^2}\frac{(n1)(nm)b}{t^{mn+2}}\right),$$
(32)
and from $`0\alpha `$-component of these equations one has

$$u_\alpha =\frac{(nm)(3k+1)(b_{\alpha ;\beta }^\beta b_{;\alpha })t^{n+1m}}{3m(2m)(k+1)}.$$
(33)
Now, writing down the spatial components of the Einstein equations, using the expressions (29)-(30) and the expression (32) for the energy density $`\epsilon `$, one get:

$`{\displaystyle \frac{\stackrel{~}{P}_\alpha ^\beta }{t^m}}+{\displaystyle \frac{m(23m)\delta _\alpha ^\beta }{4t^2}}+{\displaystyle \frac{(nm)(22nm)b_\alpha ^\beta }{4t^{mn+2}}}`$
$`{\displaystyle \frac{m(nm)b\delta _\alpha ^\beta }{4t^{mn+2}}}={\displaystyle \frac{(k1)\delta _\alpha ^\beta }{3k+1}}\left({\displaystyle \frac{3m(2m)}{2t^2}}{\displaystyle \frac{(n1)(nm)b}{t^{mn+2}}}\right).`$ (34)
Using the relation (5) it is easy to check that the terms proportional to $`\frac{1}{t^2}`$ in the left- and right-hand sides of Eq. (34) cancel each other. On the other hand, the only way to cancel the term $`\frac{\stackrel{~}{P}_\alpha ^\beta }{t^m}`$ is to require that the terms proportional to $`\frac{1}{t^{mn+2}}`$ behave as the term $`\frac{1}{t^m}`$, i.e.

$$n=2.$$
(35)
In this case, the condition of cancellation of terms proportional to $`\frac{1}{t^m}`$ gives the following expression for the tensor $`b_\alpha ^\beta `$:

$$b_\alpha ^\beta =\frac{4\stackrel{~}{P}_\alpha ^\beta }{m^24}+\frac{\stackrel{~}{P}\delta _\alpha ^\beta (3m^2+12m4)}{3m(m3)(m^24)}.$$
(36)
Using the relation (5) one can rewrite the expression (36) in the following form:

$$b_\alpha ^\beta =\frac{9(k+1)^2}{(3k+5)(3k+1)}\left(\stackrel{~}{P}_\alpha ^\beta +\frac{(3k^26k5)\delta _\alpha ^\beta \stackrel{~}{P}}{9k+5}\right).$$
(37)
It is easy to see that the Eq. (37) expressing the second-order correction for the spatial metric (2) is well defined for all the types of hydrodynamical matter with $`0k1`$, including the stiff matter, i.e. a fluid with the equation of state $`p=\epsilon `$.

Now, using the relation

$$\stackrel{~}{P}_\alpha ^{\beta ;\alpha }=\frac{1}{2}\stackrel{~}{P}_{;\beta },$$
(38)
and the formulae (37) and (33) we arrive to the following expression for the three-dimensional velocity $`u_\alpha `$:

$$u_\alpha =\frac{27k(k+1)^3\stackrel{~}{P}_{;\alpha }}{8(3k+5)(9k+5)}t^{3\frac{4}{3(k+1)}}.$$
(39)
Note that the velocity flow is potential. Actually, it can be shown that this important property remains in all higher orders of perturbative expansion for the quasi-isotropic solution. Similarly, the expression for the energy density of matter (32) can be rewritten as

$$8\pi G\epsilon =\frac{4}{3(k+1)^2t^2}+\frac{9(k+1)\stackrel{~}{P}}{2(9k+5)t^{\frac{4}{3(k+1)}}}.$$
(40)
It is straightforward to check that further terms of perturbative expansion (2) for the metric $`\gamma _{\alpha \beta }`$ have the following form:

$$\gamma _{\alpha \beta }=t^ma_{\alpha \beta }+t^2b_{\alpha \beta }+t^{2+(2m)}c_{\alpha \beta }+\mathrm{}.$$
(41)
Thus, we have seen that this expansion has a curious feature. The order of its first term $`t^m`$ is defined by the equation of state of the matter (3), the second order term always has the behavior $`t^2`$, while logarithmic distance between orders is equal to $`2m`$.

It is easy to understand that the quasi-isotropic expansion does work if and only if the first term in the right-hand side of Eq. (6) is smaller than the second one. Remembering that $`n=2`$ for any value of $`m`$ and using the equation (5), we get the following restriction on the parameter $`k`$ from the equation of state (3):

$$k>\frac{1}{3}.$$
(42)
Thus, for the values $`k1/3`$ the quasi-isotropic expansion at small times may not be constructed. It is interesting to notice, however, that for $`k=const<\frac{1}{3}`$, a quasi-isotropic-like solution arises as a late-time ($`t\mathrm{}`$) attractor for generic inhomogeneous evolution of space-time  (this regime is the power-law inflation  actually). Of course, perturbative expansion is made in inverse powers of $`t`$ in that case. Also, we would like to note that different aspects of relation between the quasi-isotropic expansion and other approximation schemes were considered in detail in the papers .

This work was supported by RFBR via grants No 02-02-16817 and 00-15-96699. A.K. is grateful to the CARIPLO Scientific Foundation for the financial support.


format: text

text:
# Surface Quasigeostrophic Turbulence : The Study of an Active Scalar

## I Introduction

In the Quasigeostrophic (QG) framework , a simplification of the Navier Stokes equations for describing the motion of a stratified and rapidly rotating fluid in a 3D domain, there are two classes of problems that immediately come to attention. The first (Charney type) are the ones where attention is focussed on the interior of the domain; the temperature is uniform along the boundaries and they play no dynamical role in the evolution of the system. The other class (Eady type) of problems lead to the Surface Quasigeostrophic (SQG) equations. The potential vorticity in the 3D interior is forced to be zero and the dynamical problem is controlled by the evolution of the potential temperature at the 2D boundaries. Working with a single lower boundary (assuming all fields to be well behaved as $`z\mathrm{}`$), the equations making up the SQG system can be expressed as ,,

$`_t\theta +u^i_i\theta =D_t\theta =0,z=0`$
$`\theta =_z\psi =(\mathrm{})^{\frac{1}{2}}\psi `$
$`\text{where}^2\psi =0,z>0\text{and}\stackrel{}{u}=^{}\psi `$ (1)
Here $`\theta `$ is the potential temperature (it is a dynamically active scalar due to the coupling of $`\theta `$ to $`\psi `$), $`\psi `$ is the geostrophic streamfunction, $`\stackrel{}{u}`$ the geostrophic velocity, $`^{}(_y,_x)`$, $`\mathrm{}`$ is the horizontal Laplacian, $`^2`$ is the full 3D Laplacian, the operator $`(\mathrm{})^{\frac{1}{2}}`$ is defined  in Fourier space via $`\widehat{(\mathrm{})^{\frac{1}{2}}\psi (k)}=|k|\widehat{\psi }(k)`$ and $`i=x,y`$. Recall that the 2D Euler equations (for an incompressible fluid) in vorticity form are,

$`_t\xi +u^i_i\xi =D_t\xi =0`$
$`\xi =\mathrm{}\psi \text{with}\stackrel{}{u}=^{}\psi `$ (2)
where $`\xi `$ is the vertical component vorticity. The similarity in the evolution equations for $`\theta `$ and $`\xi `$ has been explored in detail by Constantin, Majda and Tabak ,. It can be seen that the structure of conserved quantities in both equations is exactly the same. To be precise, just as $`f(\xi )`$,$`\psi \xi `$ are conserved by the 2D Euler equations similarly $`f(\theta )`$,$`\psi \theta `$ are conserved in the SQG system. The basic difference in the above two systems is the degree of locality of the active scalar. For the 2D Euler equations the free space Green’s function behaves as $`\mathrm{ln}(r)`$ implying a $`\frac{1}{r}`$ behaviour for the velocity field due to point vortex at the origin. In contrast, in the SQG equations the free space Green’s function has the form $`\frac{1}{r}`$ implying a much more rapidly decaying $`\frac{1}{r^2}`$ velocity field due to a point "PT vortex" at the origin . Or, in Fourier space one has $`\widehat{\xi }=|k|^2\widehat{\psi }`$ and $`\widehat{\theta }=|k|\widehat{\psi }`$ for the 2D Euler and SQG equations respectively . Hence the nature of interactions is much more local in the SQG case as compared to the 2D Euler equations.

Studying the properties of active scalars with different degrees of locality would be an interesting question in its own right  but the specific interest in the SQG equations comes from an analogy with the 3D Euler equations. This can be seen by a comparison with the 3D Euler equations, which in vorticity form read,

$$[\frac{D\omega }{Dt}]^j=_iv^j\omega ^i$$
(3)
where $`\stackrel{}{v}`$ is a divergence free velocity field, $`\stackrel{}{\omega }(=\times \stackrel{}{v})`$ is the vorticity and $`i,j=x,y,z`$. Introducing $`\stackrel{}{V}=\stackrel{}{^{}\theta }`$, a "vorticity" like quantity for the SQG system which satisfies (differentiating Eq. (1) and using incompressibilty),

$$[\frac{DV}{Dt}]^j=_iu^jV^i$$
(4)
Identifying $`\stackrel{}{V}`$ in Eq. (4) with $`\stackrel{}{\omega }`$ in Eq. (3) it can be seen  that the level sets of $`\theta `$ are geometrically analogous to the vortex lines for the 3D Euler equations. Similar to the question of a finite time singularity in the 3D Euler equations (which is thought to be physically linked to the stretching of vortex tubes), in the SQG system one can think of a scenario where the intense stretching (and bunching together) of level set lines during the evolution of a front leads to the development of shocks in finite time. The issue of treating the SQG system as a testing ground for finite time singularities has generated interest ,,,, in the mathematical community and the reader is referred to the aforementioned papers for details regarding this issue.

In view of the similarities between the SQG and 2D Euler equations and the level set stretching analogy with the 3D Euler system, it is natural to inquire into the statistical/geometrical properties of the SQG active scalar within an appropriately defined "inertial range". The broad aim is to compare these properties with the large body of work available for the 2D and 3D Euler equations. In the second section we examine the PT field via global (power spectrum, structure functions, ($`\beta ,g(\beta )`$) spectrum) and local (wavelets) methods. One of the few rigorous estimates that exist in fluid turbulence is that for the level set dimensions. The extraction of these dimensions and their agreement with analytical bounds is demonstrated. The third section is devoted to the examination of the dissipation field, the generalized dimensions of a measure based upon the dissipation field are calculated and commented upon. In the fourth section we focus our attention on the gradient fields, a simple derivation of the relation between the variety of scaling exponents is presented and the underlying assumptions are clearly stated. The failure of the cancellation exponent is demonstrated and a simple example is presented so as to put some of the ideas in perspective.

## II The Potential Temperature Field

### II.1 The Power Spectrum and the Structure Functions

A pseudo-spectral technique was employed to solve Eq. (1) numerically on a $`2048\times 2048`$ grid. Linear terms are handled exactly using an integrating-factor method, and nonlinear terms are handled by a third-order Adams-Bashforth scheme (fully de-aliased by the 2/3 rule method). The calculations were carried out for freely decaying turbulence. The initial conditions consisted of a large-scale random field, specifically a random-phase superposition of sinusoids with total wavenumber approximately equal to 6, in units where the gravest mode has unit wavenumber. Potential temperature variance is dissipated at small scales by $`\nu ^2\theta `$ diffusion. Based on the typical velocity $`U`$ and scale $`L`$ of the initial condition, one may define a Peclet number $`UL/\nu `$. The calculations analyzed here were carried out for a Peclet number of 2500. After a short time, the spectrum develops a distinct inertial range. As time progresses, energy and $`\theta `$ variance are dissipated at small scales, the amplitude decreases, and the effective Peclet number also decreases. After sufficient time, the flow becomes diffusion-dominated and the inertial range is lost. Analysis of other cases, not presented here, indicates that the results are not sensitive to the time slice or the Peclet number, so long as the Peclet number is sufficiently large and the time slice is taken at a time when there is an extensive inertial range.

The mean 1D power spectra from different stages of evolution can be seen in Fig. 1. As these are decaying simulations the structure in the PT field is slowly dying out. The resulting increase in smoothness of the PT field can be seen via the roll off of the spectrum during the later stages. In spite of this a fairly clean power law is visible for a sizeable "inertial range" (other runs with large scale initial conditions posessing various amounts of energy show similar behavior). We choose to concentrate on the particular stage which has the largest inertial range. The 2D power spectrum for this stage can be seen in Fig. 2 and a snapshot of the PT field itself can be seen in Fig. 3. Interestingly the 2D power spectrum seems to roll off at larger wavenumbers as compared to the 1D spectrum. In this stage the spectral slope (from the 1D spectra) between the scales 256 to 8 (the scales are in terms of grid size) is $`2.15`$ (the other runs also showed slopes steeper than -2). The slope from the 2D spectrum is $`2.11`$ (due to the early roll off of the 2D spectrum, this slope is extracted between the scales 128 to 8). Previous decaying simulations  obtained values near $`2`$ and seem to be consistent with our observations. A slope as steep as this suggests that the field being examined is smoother than expectations from a similarity hypothesis (which yields a $`\frac{5}{3}`$ slope ). A closer look indicates (Fig. 3) that the field is composed of a small number of "coherent structures" superposed upon a background which has a filamentary structure consisting of very fine scales. This immediately brings to mind the studies on vorticity in decaying 2D turbulence ,, wherein a similar coherent structure/background picture was found to exist. Further analysis indicated that the vorticity field possessed normal scaling whereas a measure based upon the gradient of the vorticity (precisely the enstrophy dissipation) was multifractal. To proceed in this direction we introduce the generalized structure functions of order $`q\mathrm{}^+`$,

$$S_q(r)=<|\theta (x+r)\theta (x)|^q>$$
(5)
Here $`<.>`$ represents an ensemble average. The directional dependence is suppressed due to the assumed isotropy of the PT field. Scaling behaviour in the field implies that one can expect the generalized structure functions to behave as ,

$$S_q(r)=C_1(q)|\mathrm{}\theta (L_1)|^q(\frac{r}{L_1})^{\zeta _q};r_1rL_1$$
(6)
where $`\zeta _q`$ are the generalized scaling exponents, $`C_1(q)`$ is of order unity for all $`q`$, $`|\mathrm{}\theta (L_1)|`$ is the absolute value of the difference in $`\theta `$ over a scale $`L_1`$. $`r_1`$ and $`L_1`$ are the inner and outer scale (8 and 256 respectively) over which the power law in the spectrum was observed.

If the field being examined is smooth at a scale $`r`$ then the gradient at this scale would be finite and as a consequence $`\zeta _q=q`$ (due to the domination of the linear term in the Taylor expansion about the point of interest) ie. the scaling would be trivial. Conventionally normal scaling is a term reserved for linear $`\zeta _q`$ and any nonlinearity in $`\zeta _q`$ is referred to as anomalous scaling. In 2D turbulence the velocity field is known to be smooth for all time if the initial conditions are smooth  and hence  $`\zeta _q`$(velocity)$`=q`$. Also, as mentioned, from the analysis of the vorticity field  the scaling exponents for the vorticity structure functions were found to depend on $`q`$ in a linear fashion. Plots of $`\text{log}(S_q(r))`$ Vs $`\text{log}(r)`$ for the PT can be seen in the upper panel of Fig. 4. In all cases the scaling is valid upto $`r128`$, using these plots we extracted $`\zeta _q`$ which are presented in the lower panel of Fig. 4. It is seen that the scaling is anomalous and in fact a best fit to the scaling exponents is of the form $`\zeta _q=Aq^B,q>0`$ with $`B=0.82`$.

For the special case of $`q=2`$, one can in principle relate the scaling exponent $`\zeta _2`$ to the slope of the power spectrum ($`n`$) via  $`n=(1+\zeta _2)`$. This relation is only valid for $`3<n<1`$ (note that this does not prevent the spectral slope from being steeper than $`3`$; it just implies that $`\zeta _2`$ saturates at $`2`$ for smooth fields and the particular relation between $`\zeta _2`$ and $`n`$ breaks down). In our case $`\zeta _2=1.05`$ so the predicted spectral slope is $`n=2.05`$ which is near the observed mean value of $`2.15`$ (or $`2.11`$ from the 2D spectrum). Even though the scaling exponents give an idea of the roughness in the field (anomalous scaling implying differing degrees of roughness) there is a certain unsatisfactory aspect about the structure functions, namely, there is no estimate of "how much" of the field is rough. The following subsection aims to address this very issue.

### II.2 The $`(\beta ,g(\beta ))`$ spectrum

In scaling literature the roughness of a field is specified by means of an exponent $`\beta `$ ($`>0`$) defined as ,

$$|\theta (x+r)\theta (x)||r|^{\beta (x)}$$
(7)
Here $`\beta `$ is a function of position and it refers to the fact that the derivative of $`\theta `$ will be unbounded as $`r0`$ if $`\beta <1`$. As mentioned previously there is a lower scale associated with the problem so technically nothing is blowing up and in effect $`\beta <1`$ represents the regions where the derivative will be large as compared to the rest of the field. Note that $`\theta `$ itself cannot be singular due its conserved nature. The focus is on whether an initially smooth $`\theta `$ field becomes rough so as to cause the gradient fields to experience a singularity. It is clear that Eq. (7) is by itself of not much use in characterizing $`\theta `$ (as $`\beta `$ depends upon position), in fact a global view of the specific degrees of rougness of $`\theta `$ can be attained via the $`(\beta ,g(\beta ))`$ spectrum which we introduce next.

An iso-$`\beta `$ set is defined as the set of all $`x`$’s where $`\beta (x)=\beta `$ and $`g(\beta )`$ is the dimension (to be precise $`g(\beta )`$ should be viewed in an probabilistic fashion ) of an iso-$`\beta `$ set. The dimensions of iso-$`\beta `$ sets are derived by Frisch . Briefly, at a scale $`r`$ the probability of encountering a particular value of $`\beta `$ is proportional to $`r^{dg(\beta )}`$ (where $`d=2`$ for the 2D PT field and $`d=1`$ for 1D cuts of the PT field). By using a steepest descent argument in the integral for the expectation value of $`|\theta (x+r)\theta (x)|^q`$ one obtains ,,

$$\zeta _q=\underset{\beta }{\mathrm{min}}[q\beta +dg(\beta )]\text{or}g(\beta )=\underset{q}{\mathrm{max}}[q\beta +d\zeta _q]$$
(8)
Hence given $`\zeta _q`$ for a fixed $`q=q_{}`$ using the first part of Eq. (8),

$$q_{}=\frac{dg(\beta )}{d\beta }$$
(9)
Denoting the value of $`\beta `$ for which Eq. (9) is satisfied by $`\beta _{}`$ we have,

$$g(\beta _{})=d+q_{}\beta _{}\zeta _q_{}$$
(10)
Notice that, as the structure functions involve moments with positive $`q`$, they only pick out $`\beta `$’s such that $`\beta <1`$ and $`g(\beta )<d`$. The $`(\beta ,g(\beta ))`$ spectrum seen in Fig. 5 hints at a hierarchy in roughness of the PT field. From the calculations we see that $`\beta [0.26,0.6]`$. In Fig. 4 along with $`\zeta _q`$ we have plotted the lines corresponding to $`\zeta _q=q\beta _{min}`$ and $`\zeta _q=q\beta _{max}`$. As is expected these lines straddle the actual scaling exponents. For smaller values of $`q`$ the scaling exponents are close to $`q\beta _{max}`$ as $`\beta _{max}`$ has the largest associated dimension whereas for higher values of $`q`$ the scaling exponents reflect the roughest regions and hence tend towards $`q\beta _{min}`$. Note that from Eq. (9) and Eq. (10) the approximation $`\zeta _q=Aq^B,B<1,q>0`$ implies that $`\beta `$ becomes large and $`g(\beta )d`$ as $`q0`$. In essence the picture that emerges is that even though $`\theta `$ appears to become rough (with differing degress of roughness), in fact, it is the smooth regions that occupy most of the space.

### II.3 The roughness of $`\theta `$ : a qualitative local view

All of the previous tools, the power spectrum, the $`(\beta ,g(\beta ))`$ spectrum and the structure functions extracted information in a global sense. To get a picture of the actual positions where the signal may have unbounded derivatives, and to get a qualitative feel of the spareseness of these regions, one has to determine the local behaviour of the signal in question. Recently the use of wavelets has allowed the identification of local Holder exponents in a variety of signals. The Holder exponents are extracted by a technique known as the wavelet transform modulus maxima (WTMM) method ,, ,. The modulus maxima refers to the spatial distribution of the local maxima (of the modulus) of the wavelet transform. In a crude sense the previous methods used ensemble averages of the moments of differences in $`\theta (x)`$ as "mathematical microscopes" whereas in the wavelet method it is the scale parameter of the wavelet transform that performs this task.

By using wavelets whose higher moments vanish one can detect sigularities in the higher order derivatives of the signal being analyzed . Our previous results indicate that the PT field becomes rough, i.e., the first derivatives of the PT field should be unbounded. Hence, the particular wavelet we use is the first derivative of a Gaussian whose first moment vanishes ,i.e., it picks up points where the signal becomes rough. The wavelet transform of a 1D cut of the PT field can be seen in Fig. 6. The cone like features imply the presence of a rough spot . The modulus maxima lines are extracted from this transformed field and can be seen in Fig. 7. The value of the local Holder exponent can be extracted via a log plot of the magnitude of a particular maxima line . In essence, the presence of the cones in the wavelet transform indicate roughness in the PT field and the WTMM lines locate the positions of the rough spots. However, a closer look (the lower panel of Fig. 7) suggests, qualitatively, that the rough regions are sparsely distributed (for eg. comparing with Fig. 8 in Arneodo et. al ). This goes along with the observation in the last subsection that the rough regions were non space filling. In contrast, it has been found ,, that the local Holder exponents ($`h(x)`$) associated with a 1D cut of the velocity field in 3D turbulence satisfy $`0.3h(x)0.7`$ for almost all $`x`$ (the peak of the histogram being at $`\frac{1}{3}`$) implying that the velocity field in 3D turbulence is very rough or has unbounded first derivatives at almost all points.

### II.4 The Dimension of Level Sets

Apart from being physically interesting, the level set (iso-$`\theta `$ contour) dimensions provide another link to the roughness of the scalar field. As the initial condition was a smooth 2D field, initially, any given level set ($`E_{\theta _0}`$ for $`\theta =\theta _0`$) is a non space filling curve, i.e., the level set dimension ($`D(E_{\theta _0})`$) is one. In the case of 3D turbulence there exist analytical estimates of the scalar level set dimensions . These estimates are seen to be numerically satisfied by a host of fields (both passive and active) in isotropic 3D and Magnetohydrodynamic (MHD) turbulence . It must be emphasized that these estimates come directly from the equations of evolution and are much more powerful than the phenomenological ideas we have been working with so far. The actual calculation , is of the area of an isosurface contained in a ball of specified size, given a Holder condition on the velocity field,

$$|u(x+y,t)u(x,t)|U_L(\frac{y}{L})^{\zeta _u}$$
(11)
this area estimate leads to the bound,

$$D(E_\theta )2.5+\frac{\zeta _u}{2}$$
(12)
The bound in Eq. (12) is expected to be saturated  above a certain cutoff scale. A valid extrapolation  for the level sets of PT in the SQG system reads, $`D(E_\theta )1.5+(\zeta _1)/2`$, where $`\zeta _1`$ is given by Eq. (6) due to the equality of the scaling exponents for the velocity field and the PT in the SQG system.

In passing we mention that the analytical estimates are for the Hausdorff dimension whereas practically we compute the box counting dimension. We performed calculations for a variety of nominal (i.e., near the mean) level sets and Fig. 8 shows the log-log plots used in the calculation of the box counting dimension. As is the case for 3D , there appears to be a crossover in $`D(E_\theta )`$, we find that the even though the dimension undergoes a change, the fractal nature seems to persist at smaller scales. For small $`r`$, $`D(E_\theta )1.3`$ whereas for larger values of $`r`$ we find $`D(E_\theta )1.74`$ (see the caption of Fig. 8 for details on the actual values of $`r`$). \>From our previous calculations $`\zeta _1=0.57`$, hence the analytical prediction is $`D(E_\theta )1.785`$. Furthermore, as the bound is expected to be saturated above the cutoff we see that the computed value of $`D(E_\theta )`$ for large $`r`$ is quite close to the analytical prediction. In all, apart from satisfying the analytically prescribed bounds (and indirectly indicating roughness in the PT field), the level set dimensions indicate that initially non space filling level sets acquire a fractal nature in finite time.

## III Gradient Field Characterization

The effect of the inferred roughness in the PT field will, as mentioned previously, be reflected in the singular nature of the gradient fields. In this section we proceed to examine a variety of fields which are functions of the PT gradient. The aim is to see if we can actually detect the expected singularities, and if so, to characterize them.

### III.1 The Dissipation Field

A physically interesting function of the gradient field is the PT dissipation, as it is connected to the variance of the PT field. The equation for the dissipation of $`\theta `$ can be obtained by multiplying Eq. (1) by $`\theta `$ and averaging over the whole domain,

$$\frac{<\theta ^2>}{t}=2\nu <(\theta )^2>$$
(13)
Here $`\nu (\theta )^2`$ is the dissipation field. Now consider the quantity $`\mu (x,r)`$,

$$\mu (x,r)=\frac{1}{r^d}_{B(x,r)}\nu (\theta )^2d^dx$$
(14)
Physically this is the average dissipation in a ball of size $`r`$ centered at $`x`$. Due to the smoothing via integration it is expected that $`\mu (x,r)`$ will be fairly well behaved through most of the domain with intermittent bursts of high values concentrated in the regions where the PT field is rough. The multifractal formalism , , provides a convenient way to characterize such "erratic" or singular measures. The technique  consists of constructing a measure ($`\mu `$ with suitable normalization) and using its moments to focus on the singularities of the measure. The domain in which the field is defined is partitioned into disjoint boxes of size $`r`$ and it is postulated (see for eg. ) that moments of $`\mu `$ will scale as,

$$\mu (x,r)^qr^{\tau _qqd}$$
(15)
where the sum goes over all the boxes into which the domain was partitioned. Consider the set of points where the measure scales $`r^\alpha `$ and denote the dimension of this set by $`f(\alpha )`$ (again a probabilistic view is more precise). By similar considerations as for the $`g(\beta )`$ spectrum it can be seen that ,

$$\tau _q=\underset{\alpha }{\mathrm{min}}(q\alpha f(\alpha ))\text{and}f(\alpha )=\underset{q}{\mathrm{max}}(q\alpha \tau _q)$$
(16)
The function $`\tau _q`$ is further related to the generalized dimensions  via $`D_q=\tau _q/(q1)`$. Practically  a log-log plot of the ensemble average of $`\mu (x,r)^q`$ for different values of $`r`$ gives $`(D_qd)(q1)`$. By knowing $`(q,D_q)`$ one can use Eq. (16) to obtain the $`(\alpha ,f(\alpha ))`$ spectrum.

Again, as we expect the scaling of any physical quantity to be restricted to a range of length scales (say $`r_a`$ to $`r_b`$), it is preferable to work in terms of ratios of $`\frac{r}{r_b}`$ where $`r_b`$ is the outer scale, $`r_a`$ is the inner scale and $`r_brr_a`$. In these terms the generalized dimensions can be expressed as ,

$$<\mu (r)^q>=C_2(q)(\mu (r_b)^q)(\frac{r}{r_b})^{(D_qd)(q1)}$$
(17)
where $`\mu (r_b)`$ is the measure on the outer scale and $`C_2(q)`$ is order unity for all $`q`$ (a similar relation would hold if we replaced the $`(\mu (r_b)^q)`$ in the right hand side by $`<\mu (r)>^q`$ but with different $`C_2(q)`$’s). For $`r<r_a`$ the dissipation field is assumed to have become smooth via the action of viscosity. The $`f(\alpha )`$ spectrum and the generalized dimensions for the dissipation field can be seen in Fig. 9. The nontrivial behaviour of the generalized dimensions demonstrates that the dissipation field is singular (with different singularity strengths) within these scales. Also, from the calculations we find that $`f(1)<2`$ which acts as a check that the dissipation in a finite volume is bounded (as $`f(1)`$ is the dimension of the support of the singular regions) .

### III.2 General Considerations

The gradient squared nature of the dissipation field implied that it was positive definite. In principle one could conceive of fields which possess singularities but aren’t sign definite. In order to characterize this possible sign indefiniteness Ott and co-workers studied , sign singular measures and a related family of exponents called the cancellation exponents ($`\kappa _q`$). It is immediately clear that singular nature is a prerequisite for the phenomenon of cancellation, the reason being that a nonsingular field is bounded and we can always add a constant so as to make the field positive definite and hence eliminate cancellation. Consider a sign indefinite 1D field $`\theta (x)`$ (which will later be interpreted as a 1D cut of the PT field) and construct,

$$\mu ^{}(x,r)=(\frac{1}{r}_x^{x+r}|\frac{d\theta }{dx^{}}|𝑑x^{});$$
(18)
$$\eta (x,r)=|\frac{1}{r}_x^{x+r}(\frac{d\theta }{dx^{}})𝑑x^{}|$$
(19)
As $`\mu ^{}`$ is defined for the magnitude of the gradient field we can postulate (similar to Eq. (17) for $`\mu `$),

$$<\mu ^{}(r)^q>=C_3(q)(\mu ^{}(r_d)^q)(\frac{r}{r_d})^{(D_q^{}d)(q1)}$$
(20)
Here $`d=1`$ and $`D_q^{}`$ are the generalized dimensions associated with $`\mu ^{}`$. As usual the scaling is restricted to a range of scales ($`r_c<r<r_d`$) and $`r_c`$ is the scale below which $`\mu ^{}`$ appears smooth. Formally, the entity $`\eta `$ after suitable normalization is a sign singular measure (see Ott et. al.  for a rigorous definition). It was conjectured  that the sign singular entity $`\eta `$ might also possess scaling properties in analogy with $`\mu ^{}`$, implying,

$$<\eta (r)^q>=C_4(q)(r)^{\kappa _q^{}};r>r_{cc}$$
(21)
Using Eq. (20) we can express Eq. (21) as,

$$<\eta (r)^q><\mu ^{}(r)^q>(r)^{\kappa _q};r_d>r>r_{cc}r_c$$
(22)
Where $`\kappa _q=\kappa _q^{}+(D_q^{}d)(q1)`$ and $`r_{cc}`$ is the lower oscillatory scale below which the derivative does not oscillate. We prefer to call $`\kappa _q`$ from Eq. (22) the cancellation exponents as they directly reflect the difference in the scaling properties of $`\frac{d\theta }{dx}`$ and $`|\frac{d\theta }{dx}|`$.

Apriori there is no justification in assuming that the scale at which cancellation ceases ($`r_{cc}`$) is the same as the scale at which $`\mu ^{}`$ becomes smooth ($`r_c`$). Consider the example where the derivative is a discrete signal composed of a train of delta functions (zero elsewhere) where the minimum separation between the delta functions is $`l`$. Furthermore let us assign the sign of the delta functions in a random fashion. In this case the $`r_{cc}=l`$, in fact, $`\kappa _1=0.5`$ due to the random distribution of the signs. But as the delta functions are supported at points we have $`r_c0`$. As an aside we point out that if there is a maximum scale of separation between the delta functions (say $`L`$) then at scales greater than $`L`$ we will see $`D_q^{}=d`$ for $`\mu ^{}`$. The reason for pointing this out is to give a feel for fields that exhibit scaling, on one hand smooth fields have $`D_q^{}=d`$ whereas on the other extreme random fields with small correlations also have $`D_q^{}=d`$ (at scales larger than their correlation lengths). Fields with nontrivial scaling over a significant range are in effect random but with large correlation lengths. The reader is referred to Marshak et. al.  for a detailed examination of multiplicative processes and the resulting characterization by structure functions and generalized dimensions.

In order to get a unified picture of scaling in both the gradient and the field itself there have been attempts (for eg. ,) to link $`\zeta _q`$, $`\kappa _q`$ and $`D_q^{}`$ to each other. The view that seems to have emerged is that there exist simple relations linking the various exponents and that these relations are valid under very general conditions. We present an alternate derivation of some of these relations which makes the implicit assumptions explicit. Proceeding from Eq. (22), assuming that $`\frac{d\theta }{dx}`$ has integrable singularities we obtain,

$$<|\theta (x+r)\theta (x)|^q>=<\mu ^{}(r)^q>r^{q\kappa _q}$$
(23)
Substituting from Eq. (20) and on comparing with Eq. (6) we get,

$$\zeta _q=(D_q^{}d)(q1)+q\kappa _q$$
(24)
Writing Eq. (23) for $`q=1`$ we have,

$$<|\theta (x+r)\theta (x)|>=<\mu ^{}(r)>r^{1\kappa _1}$$
(25)
which yields $`\kappa _1+\zeta _1=1`$. Now if we make a strong assumption regarding the uniform nature of the cancellation, it is possible to claim that Eq. (25) holds not only in average but on every interval, ie,

$$|\theta (x+r)\theta (x)|=\mu ^{}(r)r^{1\kappa _1}$$
(26)
Raising this to the $`q^{th}`$ power and performing an ensemble average yields,

$$<|\theta (x+r)\theta (x)|^q>=<\mu ^{}(r)^q>r^{q(1\kappa _1)}$$
(27)
which implies,

$$\zeta _q=(D_q^{}d)(q1)+q(1\kappa _1)$$
(28)
The implications of Eq. (28) are quite severe in that it shows $`\zeta _q`$ to be dependent on $`D_q^{}`$ and the knowledge of only the first cancellation exponent allows the derivation of one from the other. In general the scaling exponents of $`\theta `$ and the generalized dimensions of $`\mu ^{}`$ provide exclusive information. It is only in the presence of integrable singularities that one can link the two via Eq. (24), furthermore the stronger relation (Eq. (28)) is valid under the added assumption of uniform cancellation.

### III.3 The Gradient of the PT and its absolute value

We proceed to check if scaling is observed (as postulated) for the PT field gradient and its absolute value and whether one can extract the aforementioned exponents. In the upper panel of Fig. 10 we show the log-log plots of $`<\mu ^{}(r)^q>`$ Vs. $`r`$ for different $`q`$. The scaling relations certainly appear to hold true (which was expected as they held for the dissipation field). The generalized dimensions for $`\mu ^{}`$ can be seen the lower panel of Fig. 10. On the other hand the scaling for $`<\eta (r)>`$ seen in Fig. 11 fails to exhibit a power law in $`r`$. Hence there is no meaningful way of extracting the cancellation exponents as they have been defined. Unfortunately this implies that the relations derived in the previous section (Eqs. (24) and (28)) cannot be used in this situation. In spite of this, we can see that as $`r`$ decreases the tangent to $`\text{log}(<\eta (r)>)`$ has a smaller slope which is consitent with the existence of a oscillatory cutoff at small scales.

## IV Conclusion and Discussion

In summary, we have found that the PT field in the SQG equations appears to become rough within a specified range of scales. Moreover, not only is there a heirarchy in the degree of roughness, the roughness is distributed sparsely in a qualitative sense. These conclusions are based on a combination of factors, namely, the algebraic power spectrum, anomalous scaling in the structure functions, a nontrivial $`(\beta ,g(\beta ))`$ spectrum, the nature of the WTMM map and the wrinkling of the PT level sets. The roughness in the PT field is expected to have an adverse effect on functions of the gradient field. This expectation is bourne out in the multifractal nature of the dissipation field. Also, the singular nature of the gradient field in combination with its sign indefinetness led us to examine a sign singular measure based upon the gradient field. The failure to observe scaling in the sign singular measure serves, in our opinion, as a reminder that most scaling arguments are postulated at a phenomenological level and the underlying basis of why scaling is observed in the first place is a fairly subtle and unsettled issue. Similarly the simple derivation of the relation between the variety of scaling exponents makes explicit some of the assumptions that are required for the validity of similar relations proposed in earlier studies.

Regarding the more general question we posed in the very beginning of this paper, i.e., where does the SQG active scalar stand with respect to both 2D and 3D turbulence, we have the following comments. With respect to the vorticity in the 2D Euler equations the corresponding quantity in the SQG system is the PT. The coherent structure/fine background nature of the vorticity field ,, carries over qualitatively to the PT field. The vorticity structure functions which showed normal scaling in the 2D Euler system  cross over to anomalous scaling in the SQG system. Our conjecture is that the stronger local interactions in the SQG system are responsible for this anomalous scaling. In both these systems the vorticity and PT respectively are conserved quantities and hence any singularities one might experience are actually in the gradient fields, as is seen in the multifractal nature of the enstrophy dissipation in the Euler system  and the PT dissipation in the SQG system.

In the 3D Euler case, the PT from the SQG equations is analogous to the velocity field and $`\stackrel{}{^{}\theta }`$ from SQG corresponds to the vorticity of the 3D Euler equations. The roughness of the PT field is similar to the postulated roughness of the velocity field in the inertial range, a subtle difference being that the roughness in the PT appears to be sparse whereas indications are that the roughness in the 3D velocity field is present almost everywhere ,,. We re-emphasize that in both these cases the roughness is restricted to a range of scales, i.e., no claim is made for an actual singularity in the corresponding gradients. Similarly the anomalous scaling of the PT follows that of the velocity field in 3D but again it is not as strong as in the 3D case. The general theory developed  for the deformation of scalar level sets in the 3D case is seen to carry over to the SQG equations. In essence the SQG equations follow the 3D Euler equations but in a somewhat weaker sense. This "weakness" is clearly manifested in the behaviour of the gradients. In the 3D Euler equations a sign singular measure constructed from the vorticity field shows good scaling properties and a cancellation exponent can be meaningfully extracted , whereas in the SQG equations a similarly constructed entity lacks scaling.

The results presented in this paper have been for the most part diagnostic, in that they characterize the nature of the roughness of the potential temperature field in SQG dynamics. Although we have exhibited anomalous scaling of the potential temperature fluctuations, we do not have a theory accounting for the observed form of $`\zeta _q`$. Arriving at such a theory will be a major challenge for future work. Our results point efforts in the direction of considering the dichotomy between smooth fields within large organized vortices, and a rather sparse set at the boundaries of and between vortices which exhibits a greater degree of roughness. The diagnostic results also suggest a means for distinguishing between SQG and Euler dynamics in Nature, in cases where only a tracer field can be observed, as in the gas giant planets (notably Saturn, Jupiter and Neptune, which exibit a rich variety of turbulent patterns). SQG dynamics should yield anomalous scaling corresponding to sparse roughness, whereas Euler dynamics should yield normal scaling.

###### Acknowledgements.

We thank Prof. F. Cattaneo for performing the actual numerical simulations of the SQG system. One of the authors (J.S.) would like to thank Prof. N. Nakamura for advice regarding the manuscript and Prof. P. Constantin for clarifying certain points. Wavelab 805 (http://www-stat.stanford.edu/ wavelab/) was used for certain parts of the wavelet analysis. This project is supported by the ASCI Flash Center at the University of Chicago under DOE contract B341495


text:
# Optical qubit generation by state truncation using an experimentally feasible scheme

## 1 Introduction

In recent years, many experiments falling into the quantum domain of optical fields have been performed in optics laboratories. The field has received much attention motivated not only by the excitement of studying the fundamentals of quantum mechanics, but also by the potential use of quantum optical light for the manipulation and transmission of information. Consequently, the quantum engineering of light, which consists of generation and measurement of non-classical states of the optical field, has seen a rapid development . Several schemes, which include the Fock filtering scheme of D’Ariano et al. , photon adding and displacing scheme of Dakna et al.  and the optical state truncation of Pegg et al.  have been proposed and studied.

The optical state truncation scheme, which is also referred to as the quantum scissors device (QSD), was originally proposed for preparing superposition of vacuum and one-photon state by truncating a coherent light, and it has been modified by Koniorczyk et al.  to generate superposition of vacuum, one-photon and two-photon states. The original QSD scheme is built by two beam splitters, one of which is fed by one-photon in one input port whereas the second port is left at vacuum. One of the output ports of this beam splitter is fed to the second beam splitter where it is mixed with the coherent light. The output modes of the second beam splitter are detected and the condition in which one photon is detected in one of the modes and none in the other mode corresponds to the conditional preparation of vacuum and one-photon states at the output of the first beam splitter. Paris  has modified this scheme by replacing each beam splitter by a Mach-Zehnder interferometer and proposed a fully interferometric correspondent of the QSD in which relative weights of the vacuum and single-photon states can be tuned by varying the internal phase-shift of the interferometers.

In a recent study , we have proposed an experimental scheme for the practical realization of the QSD scheme for state truncation taking into account the realistic description of single-photon-state generation and photon counting detectors. In this paper, that study is extended to include generation of arbitrary superposition of vacuum and one-photon states and to discuss the trade-off between the fidelity and relative weights of the vacuum and one-photon states in the superposition. The effect of detector efficiency is described and the rate of preparation of the desired state is discussed. This analysis goes beyond the analyses of  by considering the realistic descriptions of photodetection and single-photon generation.

## 2 Experimental scheme for qubit generation

The device proposed is schematically depicted in Fig.1. This scheme is based on the ideas developed in . It consists of a parametric down conversion crystal as the single-photon source, conventional photon counters for conditional measurement and 50:50 beam splitters for generation of entangled photon number states (BS1) and for the mixing of coherent state with the entangled state (BS2). The overall input to the QSD scheme is

$$\widehat{\rho }_{\mathrm{in1}}=\widehat{\rho }_{(a_1,c_1)}|0_{a_2a_2}0||\alpha _{b_3b_3}\alpha |$$
(1)
where $`|\alpha _{b_3}`$ is the coherent state to be truncated to prepare the desired qubit state, and $`\widehat{\rho }_{(a_1,c_1)}`$ is a mixed state density operator obtained at the output of the parametric down conversion after averaging over all possible phases and is given as

$$\widehat{\rho }_{(a_1,c_1)}=(1\gamma ^2)\left[|0000|+\gamma ^2|1111|+\gamma ^4|2222|+\mathrm{}\right]_{(a_1,c_1)}$$
(2)
where $`\gamma ^2`$, typically $`10^4`$ , corresponds to the rate of one photon pair generation per pulse of the pump field. The action of a beam splitter (BS1 and BS2) is represented by a unitary operator $`\widehat{U}`$. Given the input density operator $`\widehat{\rho _1}`$ of a lossless 50:50 beam splitter, the output can be calculated using $`\widehat{U}^{}\widehat{\rho }_1\widehat{U}`$. $`\widehat{U}`$ is conventionally given by

$`\widehat{U}`$ $`=`$ $`\mathrm{exp}[i(\psi _t\psi _r)\widehat{L}_3]\mathrm{exp}[i2\theta \widehat{L}_2]\mathrm{exp}[i(\psi _t+\psi _r)\widehat{L}_3],`$ (3)
where

$`\widehat{L}_2={\displaystyle \frac{1}{2i}}(\widehat{a}_1^{}\widehat{a}_2\widehat{a}_2^{}\widehat{a}_1),\widehat{L}_3={\displaystyle \frac{1}{2}}(\widehat{a}_1^{}\widehat{a}_1\widehat{a}_2^{}\widehat{a}_2),`$ (4)
with $`\theta =\pi /4`$. $`\psi _t=0`$ and $`\psi _r=\pi /2`$ correspond to the phases of the beam-splitter transmission and reflection parameters, respectively. After the action of the beam splitters (BS1 and BS2) on the input state $`\widehat{\rho }_{\mathrm{in1}}`$, the output state before the photodetection process becomes

$`\widehat{\rho }_{(c_1,c_2,c_3,b_1)}=\widehat{U}_2^{}\widehat{U}_1^{}(\widehat{\rho }_{(a_1,c_1)}|0_{a_2}{}_{a_2}{}^{}0||\alpha _{b_3}^{}{}_{b_3}{}^{}\alpha |)\widehat{U}_1\widehat{U}_2.`$ (5)
The normalized truncated output state density operator at mode $`\widehat{b}_1`$ of BS1 after the conditional measurement of coincidence detection (“click”) at detectors D1 and D2 and anti-coincidence (“no-click”) at D3 is obtained by

$`\widehat{\rho }_{\mathrm{out}}={\displaystyle \frac{\mathrm{Tr}_{(c_1,c_2,c_3)}(\mathrm{\Pi }_1^{c_1}\mathrm{\Pi }_1^{c_2}\mathrm{\Pi }_0^{c_3}\widehat{\rho }_{(c_1,c_2,c_3,b_1)})}{\mathrm{Tr}_{(b_1,c_1,c_2,c_3)}(\mathrm{\Pi }_1^{c_1}\mathrm{\Pi }_1^{c_2}\mathrm{\Pi }_0^{c_3}\widehat{\rho }_{(c_1,c_2,c_3,b_1)})}},`$ (6)
where $`\mathrm{\Pi }_1^{c_1}`$, $`\mathrm{\Pi }_1^{c_2}`$, and $`\mathrm{\Pi }_0^{c_3}`$ are elements of the positive-operator-valued measures (POVMs) for the detectors D1, D2, and D3, respectively, with $`0`$ and $`1`$ corresponding to the number of clicks recorded at the detectors. The measurement by conventional photon counters is described by the following two-value POVM

$`\mathrm{\Pi }_0`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(1\eta )^m|mm|,`$
$`\mathrm{\Pi }_{N1}`$ $`=`$ $`1\mathrm{\Pi }_0.`$ (7)
with $`\eta `$ being the quantum efficiency of the detectors. The effect of mean dark count $`\nu `$ of the detector in the POVM has been neglected because in a previous study, it was shown that when $`\nu \gamma ^2`$, which can be achieved in experiments using commercially available detectors, dark count rate does not have a significant effect on the fidelity and efficiency of the truncation process.

Due to non-ideal detection and non-ideal single-photon generation (output of parametric down conversion process contains vacuum and higher number of photons together with a low rate of single-photon pair), the conditional output state $`\widehat{\rho }_{\mathrm{out}}`$ is not a pure state. However, as was shown in , there are regimes of intensity of the input coherent state and the detector parameters where the output state approaches the superposition state of vacuum and one-photon with non-vanishing probability. To compare $`\widehat{\rho }_{\mathrm{out}}`$ with the desired qubit state of the form $`|\psi _{\mathrm{desired}}=N[c_0|0+c_1|1]`$ with $`N`$ being the normalization constant, we consider the fidelity $`F=\psi _{\mathrm{desired}}|\widehat{\rho }_{\mathrm{out}}|\psi _{\mathrm{desired}}`$.

## 3 Preparation of qubit states with the experimental QSD scheme

In the QSD scheme there are two kinds of free parameters, namely the intensity $`|\alpha |^2`$ of the coherent light and the beam splitter parameter (reflectance and transmittance), that can be tuned to properly set the relative weights of the vacuum and one-photon states in the superposition of $`N[c_0|0+c_1|1]`$. The tuning of the beam splitter parameters can be realized by using the interferometric scheme of Paris . This scheme makes the setup more complicated and introduces the problem of controlling the stability and balance of the interferometers. The adjustment of the $`|\alpha |^2`$ can be realized without introducing additional complexity to the setup except a controllable tuning of the intensity of the coherent light. However, this may also constitute a challenge for the preparation of some specific qubit states which needs very fine tuning of the intensity.

In this study, our aim is to optimize the intensity $`|\alpha |^2`$ of the input coherent light to generate an arbitrary qubit state with the highest fidelity with non-vanishing probability of the state generation. In that case, $`BS1`$ and $`BS2`$ are chosen as $`50:50`$ beam splitters because they give the highest probability and fidelity for state truncation using the ideal QSD scheme . The highest fidelity in generating arbitrary qubit states from a coherent state can be achieved if arg($`\alpha `$)=arg($`c_1`$), which is assumed to be the case in this study.

In the following, the results of numerical simulations for the proposed experimental QSD scheme are presented. The simulations were performed for different $`|\alpha |^2`$ and $`\eta `$ using the POVM given by (7). Parametric down conversion output was used in the simulations in the form of (2) and $`\gamma ^2`$ parameter was taken as $`4\times 10^4`$ for which the ratio of one-photon pair generation to that of the two-photon pair generation is $`O(10^4)`$ resulting in a low probability of having contributions from higher number of photons at the output mixed state.

### 3.1 Equally weighted superposition of vacuum and one-photon states

Effects of detection efficiency and the intensity of the coherent light on the preparation of the qubit state which consists of equally weighted vacuum and one-photon components of the form $`\left[|0+|1\right]/\sqrt{2}`$ are analyzed . The results of these simulations are depicted in figure 2. The relation between $`\eta `$ and $`|\alpha |^2`$ to achieve a given fidelity is clearly seen in this figure. For $`|\alpha |^20.36`$, fidelity is less than $`0.9`$ for any value of detector efficiency. At much smaller values of intensity, fidelity is not affected by the detector efficiency. At a constant $`|\alpha |^2`$, the number of correct detection events increases with increasing detector efficiency. It is seen that equally weighted superposition of vacuum and single-photon states can be prepared with high fidelity even for $`\eta `$ as low as $`0.1`$ provided that the intensity of the coherent light is chosen properly. For $`\eta =0.5`$, the optimum value for $`|\alpha |^2`$ is $`0.72`$ for which the fidelity of the prepared state to the desired one becomes $`0.89`$. Such a detection and state preparation can be achieved with a rate of 4533 s-1. For each detection efficiency, the optimized value of

$`|\alpha |^2`$ and the maximum attainable fidelity are different. Increasing $`\eta `$ from $`0.5`$ to $`0.7`$ will increase the value of $`|\alpha |^2`$ to $`1.06`$ and will change the fidelity only by $`0.11\%`$. On the other hand, the generation rate of such a state is almost doubled to a value of 8524 s-1.

### 3.2 Arbitrary superposition of vacuum and one-photon states

The intensity $`|\alpha |^2`$ of the coherent light has been optimized to give the maximum fidelity to prepare an arbitrary superposition of vacuum and one-photon state in the form of

$$|\psi _{\mathrm{desired}}=\frac{c_0|0+c_1|1}{\sqrt{|c_0|^2+|c_1|^2}}.$$
(8)
Figure 3 depicts the values of the optimized $`|\alpha |^2`$ and the maximum value of fidelity that can be achieved for any desired state of $`|c_1/c_0|`$ with that value of $`|\alpha |^2`$ for different detector efficiencies. It has been understood that an optimum $`|\alpha |^2`$, which will yield a maximum fidelity with non-zero preparation rate, can always be found for any arbitrary desired state. This optimum value of $`|\alpha |^2`$ depends on the detector efficiency $`\eta `$ and $`|c_1/c_0|`$. When $`|c_1/c_0|0.4`$, a fidelity value of $`0.99`$ can be achieved with the optimized $`|\alpha |^2`$ for $`\eta 0.5`$. The fidelity of the qubit state preparation depends significantly on the relative weights of the vacuum and single-photon states in the superposition. If the vacuum component is dominant in the desired state, then high fidelity values greater than $`0.9`$ can be achieved even for smaller detector efficiencies. But when the one-photon component becomes dominant, the fidelity of the state preparation starts decreasing. For $`|c_1/c_0|=0.5`$, $`F0.98`$ at the optimized value $`|\alpha |^20.21`$ and $`\eta =0.5`$, however when the desired state is $`|c_1/c_0|=2`$, the fidelity drops to 0.74 at the optimized value of $`|\alpha |^22.7`$. If $`\eta `$ is increased, fidelity increases to only $`0.77`$ at $`\eta =0.7`$, and to $`0.8`$ at $`\eta =0.9`$. All of these arbitrary qubit states can be generated with non-zero preparation rate. This preparation rate increases with increase in $`|c_1/c_0|`$ up to some maximum value which depends on $`\eta `$ and then starts decreasing when one-photon state becomes more dominant in the superposition.

It is observed that when detector efficiency is in the range of $`1.0\eta 0.9`$, for some values of $`|c_1/c_0|`$, two different values of optimum $`|\alpha |^2`$ at which fidelity takes the

same maximum value can be found (see figure 3(b)). Although fidelity reaches the maximum value for these two values of $`|\alpha |^2`$, the corresponding rates of the preparation of the state are different. As it can be seen in figure 3, for the preparation of the qubit state with $`|c_1/c_0|1.145`$ using detectors of $`\eta =0.9`$, the maximum fidelity of $`0.906`$ can be obtained at $`|\alpha |^21.464`$ and $`|\alpha |^28.644`$ for which the state preparation rates are 11869s-1 and 561s-1, respectively. In the same way, for $`|c_1/c_0|0.712`$, if we increase $`\eta `$ to unity, the state preparation rate will become 13330s-1 and 35s-1 at $`|\alpha |^20.548`$ and $`|\alpha |^212.532`$, respectively, with the same maximum attainable fidelity of $`0.973`$.

It must be noted that optimizing the intensity of the coherent light to obtain the maximum fidelity for a desired qubit state does not necessarily mean that the system is maximized for the highest preparation rate, too. The cost of maximizing the fidelity is the decrease in the preparation rate and vice versa. At $`\eta =0.7`$, to prepare a qubit state having $`|c_1/c_0|=0.4`$ with the highest attainable fidelity of $`0.992`$, the optimum $`|\alpha |^2=0.152`$ must be used which will result in preparation rate of 5790s-1. However, if the system is optimized to have the highest preparation rate, an optimum value of $`|\alpha |^2=1.538`$ needs to be used. This will result in a fidelity of $`0.859`$ which is much lower than that of the former case. At $`\eta 0.7`$, only for the case of balanced superposition state, the values of fidelity and state preparation rate obtained when the system is optimized for the highest fidelity are very close to those obtained when optimization is made for the highest preparation rate.

Figure 4 depicts the comparison of the ideal QSD scheme of Barnett et al.  and the proposed experimental scheme for $`|\alpha |^24`$ and $`\eta =0.5`$. We have not considered the larger values of $`|\alpha |^2`$, because at these higher values of light intensity, even though in some cases, fidelity has a higher value, the state preparation rate approaches very small values which are not experimentally feasible. The comparison has been done in three different regions of $`|c_1/c_0|`$ for various $`\eta `$: (a) For $`|c_1/c_0|<1.0`$, which means vacuum dominated qubit state, fidelity increases with increasing $`|\alpha |^2`$ until it reaches the maximum value for both schemes. After that maximum value, fidelity starts decreasing with respect to $`|\alpha |^2`$, however the experimental scheme overwhelms the results obtained in  for the ideal scheme and has much higher fidelity values. It is seen in figure 4 that for $`|c_1/c_0|=0.4`$, the fidelity values for both schemes are very close to each other until $`|\alpha |^20.28`$ beyond which fidelity of the ideal scheme decreases

more rapidly with increasing values of $`|\alpha |^2`$. (b) For the case of balanced superposition state having $`|c_1/c_0|=1.0`$, ideal QSD scheme has higher value of fidelity than the experimental scheme up to a critical value of $`|\alpha |^21.52`$ after which the experimental scheme starts to give slightly higher fidelity values; this difference is not as pronounced as the difference seen in case (a). (c) For $`|c_1/c_0|>1.0`$ corresponding to the one-photon-dominant qubit states, the ideal QSD scheme prepares states with higher fidelity than that of the proposed experimental scheme. When the total number of photons incident on the detectors is large, the output state approaches the vacuum state enhancing the advantage of the experimental scheme.

## 4 CONCLUSIONS

The experimentally feasible optical state truncation scheme proposed in  has been analyzed for the preparation of arbitrary qubit states of the form given by (8). State preparation is based on truncating a coherent light to contain only its vacuum and one-photon components. The relative weights of the vacuum and one-photon states in the qubit state can be adjusted by manipulating the intensity of the coherent light. It is understood that the qubit states, which are dominated by the vacuum component, can be prepared with high fidelity and efficiency with the proposed scheme considering the realistic descriptions of photodetection and single-photon generation. Moreover, one can always find an optimum value for the intensity of the coherent light which will maximize the fidelity. However, the price for the maximization of the fidelity is the decrease in the preparation rate of the state. This study reveals that the original QSD scheme and the proposed experimental scheme can be used not only for truncating an input coherent state but also for generating arbitrary qubit states with proper manipulation of the light intensity.

## Acknowledgments

We thank Prof. Stephen M. Barnett for pointing us this problem and for his invitation and hospitality during Quantum Electronics and Photonics 15 (QEP-15) conference. We also thank Takashi Yamamoto and Yu-xi Liu for their stimulating discussions. This work was supported by a Grant-in-Aid for Encouragement of Young Scientists (Grant No. 12740243) and a Grant-in-Aid for Scientific Research (B) (Grant No. 12440111) by Japan Society for the Promotion of Science.

## References

to appear in the Journal of Modern Optics (spring 2002).


text:
# Positivity bounds on generalized parton distributions in impact parameter representation

## I General form of the positivity bounds on GPDs

Generalized parton distributions (GPDs) MRGDH-94 ; Radyushkin-96 ; Ji-97 ; CFS-97 ; Radyushkin-97 ; Radyushkin-review ; GPV ; BMK-2001  also known as off-forward, skewed, nondiagonal etc. appear in the QCD description of various hard processes e.g. deeply virtual Compton scattering and hard exclusive meson production. Our knowledge about GPDs is poor and any additional theoretical information is of value. From this point of view the positivity bounds are rather important.

GPDs are defined in terms of matrix elements

$$G_{P_1i_1,P_2i_2}^{\alpha _1\alpha _2}(x_1,x_2)=P_2,i_2|\varphi _{\alpha _2}^{}\left(x_2\right)\varphi _{\alpha _1}\left(x_1\right)|P_1,i_1$$

$$P_2,i_2|P_1,i_10|\varphi _{\alpha _2}^{}\left(x_2\right)\varphi _{\alpha _1}\left(x_1\right)|0.$$
(1)
Here $`|P_k,i_k`$ is a hadron state (with momentum $`P_k`$ and spin/isospin indices $`i_k`$) and the field $`\varphi _{\alpha _k}\left(x_k\right)`$ describes the annihilation of a parton with momentum fraction $`x_k`$ and with spin/isospin labeled by $`\alpha _k`$. Indices $`\alpha _k`$ also contain information about the type of the parton (quark or gluon). The vacuum subtraction term in the RHS of Eq. (1) can be ignored in practical applications of GPDs with $`P_1P_2`$ but this term is important in the derivation of the positivity constraints on GPDs. Parton momentum fractions $`x_ik_i^+/P^+`$ are normalized with respect to some external fixed scale $`P^+`$ and not to the hadron momenta $`P_1`$ or $`P_2`$. Although notation (1) for GPDs differs from the standard one, we find the form (1) rather convenient for the derivation of positivity bounds and for the analysis of the interplay between the evolution and the positivity properties.

The positivity of the norm in the Hilbert space of states

$$\underset{i\alpha }{}\frac{dP^+d^2P^{}dx}{2P^+(2\pi )^3}\theta (x)f_{i\alpha }(x,P)\varphi _\alpha \left(x\right)|P,i^20$$
(2)
leads to the inequality

$$\underset{i_1\alpha _1i_2\alpha _2}{}\frac{dP_1^+d^2P_1^{}dx_1}{2P_1^+(2\pi )^3}\frac{dP_2^+d^2P_2^{}dx_2}{2P_2^+(2\pi )^3}\theta (x_1)\theta (x_2)$$

$$\times f_{i_2\alpha _2}^{}(x_2,P_2)f_{i_1\alpha _1}(x_1,P_1)G_{P_1i_1,P_2i_2}^{\alpha _1\alpha _2}(x_1,x_2)0.$$
(3)
The integration over $`x_{1,2}`$ is restricted to the positive region

$$x_1>0,x_2>0$$
(4)
where the vacuum term vanishes in the RHS of Eq. (1) so that its subtraction does not violate the positivity. In the case of the antiquark GPDs one should consider the region $`x_1,x_2<0`$.

Strictly speaking the above expressions have to be written more carefully: in the case of the gauge theories the $`P`$ exponent should be inserted between parton fields and certain conventions have to be chosen concerning the normalization of the momentum fractions $`x_1,x_2`$.

Various positivity inequalities for GPDs corresponding to specific choices of the function $`f_{i\alpha }(x,P)`$ have been already discussed in literature Martin-98 ; Radyushkin-99 ; PST-99 ; Ji-98 ; DFJK-00 ; Burkardt-01 ; Pobylitsa-01 ; Pobylitsa-02 ; Diehl-02 ; Burkardt-02-a ; Burkardt-02-b . The aim of this paper is to analyze inequality (3) with arbitrary functions $`f_{i\alpha }(x,P)`$.

## II Positivity bounds and evolution

It is well known that in the case of forward parton distributions the probabilistic interpretation of the one-loop DGLAP evolution GL-BLP ; AP-77 ; Dokshitzer-77  leads to the stability of the positivity properties under the one-loop evolution to higher normalization points. The generalization of this property for some particular positivity bounds on GPDs was considered in Ref. PST-99 . Let us show modifying the argument of Ref. PST-99  that if at some normalization point inequality (3) holds for all functions $`f_{i\alpha }(x,P)`$ then after the one-loop evolution

$$\mu \frac{}{\mu }G_{P_1i_1,P_2i_2}^{\alpha _1\alpha _2}(x_1,x_2;\mu )=g^2(\mu )𝑑y_1𝑑y_2$$

$$\times \underset{\beta _1\beta _2}{}K_{\beta _1\beta _2}^{\alpha _1\alpha _2}(x_1,x_2,y_1,y_2)G_{P_1i_1,P_2i_2}^{\beta _1\beta _2}(y_1,y_2;\mu )$$
(5)
to a higher normalization point inequality (3) is still valid for all $`f_{i\alpha }(x,P)`$.

Positivity bound (3) involves GPDs $`G_{P_1i_1,P_2i_2}^{\alpha _1\alpha _2}(x_1,x_2;\mu )`$ with positive $`x_1,x_2>0`$ (4). If $`x_1,x_2>0`$ then the one-loop evolution kernels $`K_{\beta _1\beta _2}^{\alpha _1\alpha _2}(x_1,x_2,y_1,y_2)`$ differ from zero only in the region

$$y_kx_k0.$$
(6)
This constraint has a simple physical meaning: a parton with momentum fraction $`y_k`$ can emit a parton with momentum fraction $`x_k`$ only under the condition (6).

The one-loop evolution kernels $`K_{\beta _1\beta _2}^{\alpha _1\alpha _2}(x_1,x_2,y_1,y_2)`$ for GPDs can be interpreted as perturbative parton-in-parton GPDs. For parton-in-parton GPDs one can repeat the derivation of the inequality (3) arriving at the following inequality holding for any function $`h_{\alpha \beta }(x,y)`$

$$\underset{i_k\alpha _k\beta _k}{}\underset{0}{\overset{\mathrm{}}{}}𝑑x_1\underset{x_1}{\overset{\mathrm{}}{}}𝑑y_1\underset{0}{\overset{\mathrm{}}{}}𝑑x_2\underset{x_2}{\overset{\mathrm{}}{}}𝑑y_2h_{\alpha _2\beta _2}^{}(x_2,y_2)$$

$$\times h_{\alpha _1\beta _1}(x_1,y_1)K_{\beta _1\beta _2}^{\alpha _1\alpha _2}(x_1,x_2,y_1,y_2)0.$$
(7)
The diagrammatic interpretation of this inequality is described in Appendix A.

Actually inequality (7) holds only for functions $`h_{\alpha \beta }(x,y)`$ vanishing at $`x=y`$ because of the virtual terms proportional to $`\delta (x_1y_1)\delta (x_2y_2)`$ which give a negative contribution to the kernel $`K_{\beta _1\beta _2}^{\alpha _1\alpha _2}(x_1,x_2,y_1,y_2)`$. If one ignores these terms then the stability of the positivity bound (3) with respect to the evolution (5) upwards in $`\mu `$ is a consequence of the property (7) of the one-loop evolution kernel. The last step is to notice that after the inclusion of the virtual terms proportional to $`\delta (x_1y_1)\delta (x_2y_2)`$ in the evolution equation, the positivity is still preserved under the evolution upwards. One can use the same argument as in the case of the forward distributions: if at some “critical point” $`\mu `$ inequality (3) is saturated and becomes an equality for some function $`f_{i\alpha }(x,P)`$ then at this point the virtual terms proportional to $`\delta (x_1y_1)\delta (x_2y_2)`$ do not contribute to the evolution of the LHS of inequality (3) so that at higher $`\mu `$ the positivity is again restored by the positive part of the evolution kernel. The details can be found in Appendix B.

## III Positivity bounds in the impact parameter representation

The positivity bound (3) contains a multidimensional integral in the LHS and this form is not quite convenient for practical applications. Turning to the impact parameter representation Burkardt-01 ; Diehl-02 ; Burkardt-02-a ; Burkardt-02-b  one can simplify this inequality. Let us show how this can be done for quark GPDs. The generalization for gluons is straightforward.

Let us introduce a light-cone vector $`n`$ and define the light-cone coordinates so that $`(nX)=X^+`$ for any 4-vector $`X`$. Below for simplicity we shall restrict our analysis to the case of “unpolarized” GPDs defined in terms of the matrix elements of the operator $`\overline{\psi }(n\gamma )\psi `$

$$\frac{d\lambda }{2\pi }\mathrm{exp}\left[\frac{1}{2}i\lambda (k_1^++k_2^+)\right]$$

$$\times P_2,\sigma _2|\overline{\psi }\left(\frac{\lambda n}{2}\right)(n\gamma )\psi \left(\frac{\lambda n}{2}\right)|P_1,\sigma _1$$

$$=F_{\sigma _2\sigma _1}[\frac{k_1^++k_2^+}{P_1^++P_2^+},\frac{P_1^+P_2^+}{P_1^++P_2^+},\frac{2(P_1^+P_2^{}P_2^+P_1^{})}{P_1^++P_2^+}].$$
(8)
The structure of the momentum dependence in the RHS is fixed by the Lorentz invariance. In order to see the relation between the arguments of function $`F_{\sigma _2\sigma _1}`$ and the standard notation of X. Ji Ji-98

$$\mathrm{\Delta }=P_2P_1,t=\mathrm{\Delta }^2,\xi =\frac{(nP_1)(nP_2)}{(nP_2)+(nP_1)}$$
(9)
one can use the frame where $`P_2^{}=P_1^{}\mathrm{\Delta }^{}/2`$. In this frame the RHS of Eq. (8) is simply $`F_{\sigma _2\sigma _1}(x,\xi ,\mathrm{\Delta }^{})`$.The longitudinal momentum flow corresponding to variables $`x,\xi `$ is shown in Fig. 1.

Let us insert Eq. (8) into the positivity condition (3) and turn to the impact parameter representation

$$\stackrel{~}{F}_{\sigma _2\sigma _1}(x,\xi ,b^{})=\frac{d^2\mathrm{\Delta }^{}}{(2\pi )^2}\mathrm{exp}\left[i(\mathrm{\Delta }^{}b^{})\right]$$

$$\times F_{\sigma _2\sigma _1}(x,\xi ,\mathrm{\Delta }^{}).$$
(10)
Using a factorized ansatz for the function $`f_{i_1\alpha _1}(x_1,P_1)`$ one can reduce inequality (3) to the following relatively simple form (technical details can be found in Appendix C)

$$\underset{\sigma _1\sigma _2}{}\underset{1}{\overset{1}{}}𝑑\xi \underset{|\xi |}{\overset{1}{}}𝑑x\frac{1}{(1x)^5}p_{\sigma _2}^{}\left(\frac{1x}{1\xi }\right)p_{\sigma _1}\left(\frac{1x}{1+\xi }\right)$$

$$\times \stackrel{~}{F}_{\sigma _2\sigma _1}(x,\xi ,\frac{1x}{1\xi ^2}b^{})0.$$
(11)
This inequality should be valid for any functions $`p_\sigma `$ (originating from the factorized ansatz for $`f_{i_1\alpha _1}(x_1,P_1)`$) and for any value of parameter $`b^{}`$ if one wants the original inequality (3) to hold for arbitrary $`f_{i_1\alpha _1}(x_1,P_1)`$.

For hadrons with spin 0 (e.g. pions) we define the GPD as follows

$$q(x,\xi ,t)=\frac{d\lambda }{2\pi }\mathrm{exp}\left(i\lambda x\right)$$

$$\times P_2|\overline{\psi }\left(\frac{\lambda n}{2}\right)(n\gamma )\psi \left(\frac{\lambda n}{2}\right)|P_1$$
(12)
and for hadrons with spin 1/2 we use the standard notation of X. Ji Ji-98

$$\frac{d\lambda }{2\pi }\mathrm{exp}(i\lambda x)U(P_2)|\overline{\psi }\left(\frac{\lambda n}{2}\right)(n\gamma )\psi \left(\frac{\lambda n}{2}\right)|U(P_1)$$

$$=H(x,\xi ,t)\overline{U}(P_2)(n\gamma )U(P_1)$$

$$+\frac{1}{2M}E(x,\xi ,t)\overline{U}(P_2)i\sigma ^{\mu \nu }n_\mu \mathrm{\Delta }_\nu U(P_1)$$
(13)
Here we assume the normalization condition $`n(P_1+P_2)=2`$. Note that we could not impose this condition earlier since $`P_1`$ and $`P_2`$ were integration variables in Eq. (3). We also follow the standard convention that the transverse coordinates are orthogonal to both $`n`$ and $`P_1+P_2`$. The Lorentz invariant squared momentum transfer $`t\mathrm{\Delta }^2<0`$ can be expressed in terms of the transverse part $`|\mathrm{\Delta }^{}|^2`$

$$t=(1\xi ^2)^1\left(|\mathrm{\Delta }^{}|^2+4\xi ^2M^2\right).$$
(14)
In the case of spin-0 hadrons the function $`F_{\sigma _2\sigma _1}`$ appearing in inequality (11) can be expressed in terms of GPD $`q`$ (12) as follows

$$F(x,\xi ,\mathrm{\Delta }^{})=q(x,\xi ,\frac{|\mathrm{\Delta }^{}|^2+4\xi ^2M^2}{1\xi ^2}).$$
(15)
For spin-1/2 hadrons we have

$$F_{\sigma _2\sigma _1}(x,\xi ,\mathrm{\Delta }^{})=\sqrt{1\xi ^2}$$

$$\times \left(\begin{array}{cc}2H\frac{2\xi ^2}{1\xi ^2}E& \frac{\mathrm{\Delta }^1+i\mathrm{\Delta }^2}{M(1\xi ^2)}E\\ \frac{\mathrm{\Delta }^1+i\mathrm{\Delta }^2}{M(1\xi ^2)}E& 2H\frac{2\xi ^2}{1\xi ^2}E\end{array}\right)_{\sigma _2\sigma _1}$$
(16)
where the arguments $`x,\xi `$ and $`t`$ (14) are assumed for $`H`$ and $`E`$ in the RHS.

Inequality (11) is our main result. Combining the expressions for $`F_{\sigma _2\sigma _1}`$ in terms of GPDs (15), (16) with the impact parameter representation (10) and inserting the result into inequality (11) one can obtain the explicit form of the inequalities for GPDs $`q,H,E`$.

## IV Positivity bounds on nucleon GPDs

From the practical point of view the most interesting case corresponds to the positivity bounds on nucleon GPDs. Let us rewrite our general inequality (11) for spin-1/2 hadrons in terms of the standard notations for nucleon GPDs $`H,E`$ (13).

Inserting Eqs. (10) and (16) into the general inequality (11) we find

$$\underset{\sigma _1\sigma _2}{}\underset{1}{\overset{1}{}}𝑑\xi \underset{|\xi |}{\overset{1}{}}𝑑x\frac{\sqrt{1\xi ^2}}{(1x)^5}p_{\sigma _2}^{}\left(\frac{1x}{1\xi }\right)p_{\sigma _1}\left(\frac{1x}{1+\xi }\right)$$

$$\times \frac{d^2\mathrm{\Delta }^{}}{(2\pi )^2}\mathrm{exp}\left[i\frac{1x}{1\xi ^2}(\mathrm{\Delta }^{}b^{})\right]$$

$$\times \left(\begin{array}{cc}2H\frac{2\xi ^2}{1\xi ^2}E& \frac{\mathrm{\Delta }^1+i\mathrm{\Delta }^2}{M(1\xi ^2)}E\\ \frac{\mathrm{\Delta }^1+i\mathrm{\Delta }^2}{M(1\xi ^2)}E& 2H\frac{2\xi ^2}{1\xi ^2}E\end{array}\right)_{\sigma _2\sigma _1}0.$$
(17)
The integration over $`\mathrm{\Delta }^{}`$ can be expressed in terms of the integration over variable $`t`$ (14). Since functions $`p_\sigma (z)`$ are arbitrary we can rescale them: $`p_\sigma (z)z^{3/2}p_\sigma (z)`$. Then inequality (17) takes the form

$$\underset{1}{\overset{1}{}}𝑑\xi \underset{|\xi |}{\overset{1}{}}\frac{dx}{(1x)^2}\underset{\mathrm{}}{\overset{t_0}{}}𝑑tJ_0\left(\frac{1x}{\sqrt{1\xi ^2}}|b^{}|\sqrt{t_0t}\right)$$

$$\times \{(p_{+,}^{}p_{+,+}+p_,^{}p_{,+})(H\frac{\xi ^2}{1\xi ^2}E)$$

$$+\frac{p_{+,}^{}p_{,+}+p_,^{}p_{+,+}}{|b^{}|M(1x)}\frac{}{t}\left[(t_0t)E\right]\}0.$$
(18)
Here $`J_0`$ is the Bessel function, parameter

$$t_0=\frac{4\xi ^2M^2}{1\xi ^2}$$
(19)
corresponds to the maximal kinematically allowed value of $`t`$ (14), functions $`H,E`$ are taken at their standard arguments $`x,\xi ,t,`$ and $`p_{\sigma ,\pm }`$ are the values of arbitrary functions $`p_\sigma (z)`$ at points $`z=(1x)/(1\pm \xi )`$

$$p_{\sigma ,\pm }=p_\sigma \left(\frac{1x}{1\pm \xi }\right).$$
(20)
## V Special cases

The general positivity bound (11) imposes rather serious constraints on GPDs since it should hold for any functions $`p_\sigma `$ and for any values of the impact parameter $`b^{}`$. Let us show how choosing various functions $`p_\sigma `$ one can reproduce most of the old positivity bounds and obtain new interesting results. For simplicity we consider the case of spinless hadrons (the generalization for spin-1/2 hadrons is straightforward).

1. Integrating inequality (11) over $`b^{}`$ with the weight $`|b^{}|^{2n}`$ and replacing $`p(z)z^{n+2}p(z)`$ we obtain the following inequality for GPD $`q`$ (15)

$$\underset{1}{\overset{1}{}}𝑑\xi \underset{|\xi |}{\overset{1}{}}𝑑x\frac{1}{(1x)^3}p^{}\left(\frac{1x}{1\xi }\right)p\left(\frac{1x}{1+\xi }\right)$$

$$\times \frac{^n}{t^n}q(x,\xi ,t)|_{t=4\xi ^2M^2/(1\xi ^2)}0$$
(21)
which should hold for any integer $`n=0,1,2,\mathrm{}`$ and for any function $`p`$.

2. Taking $`p(z)=\delta (zc)`$ in (11) one arrives at

$$\frac{d^2\mathrm{\Delta }^{}}{(2\pi )^2}\mathrm{exp}\left[i(\mathrm{\Delta }^{}b^{})\right]q(x,0,t)|_{t=|\mathrm{\Delta }^{}|^2}0.$$
(22)
This coincides with the inequality derived earlier in Ref. Burkardt-01 .

3. Using function

$$p\left(z\right)=c_1\delta \left(z\frac{1x}{1+\xi }\right)+c_2\delta \left(z\frac{1x}{1\xi }\right)$$
(23)
in inequality (11) we obtain

$$|c_1|^2\left(\frac{1+\xi }{1\xi }\right)^3\stackrel{~}{F}(\frac{x+\xi }{1+\xi },0,\frac{1x}{1+\xi }b^{})$$

$$+|c_2|^2\left(\frac{1\xi }{1+\xi }\right)^3\stackrel{~}{F}(\frac{x\xi }{1\xi },0,\frac{1x}{1\xi }b^{})$$

$$+\frac{c_1^{}c_2+c_2^{}c_1}{1\xi ^2}\stackrel{~}{F}(x,\xi ,\frac{1x}{1\xi ^2}b^{})0.$$
(24)
Taking into account that this inequality should hold for any $`c_1,c_2`$ and rescaling $`b^{}(1\xi ^2)(1x)^1b^{}`$ we find

$$(1\xi ^2)^1\stackrel{~}{F}(x,\xi ,b^{})$$

$$\sqrt{\stackrel{~}{F}(x_1,0,(1\xi )b^{})\stackrel{~}{F}(x_2,0,(1+\xi )b^{})}$$
(25)
where

$$x_{1,2}=(x\pm \xi )/(1\pm \xi ).$$
(26)
This inequality was derived in Ref. Diehl-02  (with a different normalization used in the definitions of GPD $`\stackrel{~}{F}`$ and of the impact parameter $`b^{}`$).

4. Taking $`c_k=d_k\mathrm{exp}\left(iP_k^{}b^{}\right)`$ in inequality (24), integrating over $`b^{}`$ and optimizing the resulting inequality with respect to arbitrary coefficients $`d_k`$ we obtain

$$q(x,\xi ,t)\sqrt{q\left(x_1\right)q\left(x_2\right)}$$
(27)
where $`x_1,x_2`$ are given by Eq. (26) and $`q\left(x_{1,2}\right)q(x_{1,2},0,0)`$ is the usual forward parton distribution. This inequality was obtained earlier in Ref. PST-99  (the authors of PST-99  use a different normalization of GPD $`q`$).

5. Using the modification of the ansatz (23) corresponding to spin-1/2 hadrons one can derive from the general inequality (11)

$$\left(H\frac{\xi ^2}{1\xi ^2}E\right)^2+\left(\frac{E}{2M(1\xi ^2)}|\mathrm{\Delta }^{}|\right)^2|_{x,\xi ,t}$$

$$\frac{q(x_1)q(x_2)}{1\xi ^2}.$$
(28)
Again relation (26) is assumed between the arguments $`x,\xi `$ in the LHS and the variables $`x_{1,2}`$ in the RHS. This inequality was derived earlier in Ref. Pobylitsa-01 .

This short list of particular cases is only a small part of the bounds that can be extracted from the general inequality (11).

## VI Positivity bounds and renormalization

GPDs are defined in terms of matrix elements (12), (13) of parton fields separated by a light-cone interval. Similarly to the case of the forward parton distributions (FPDs) this formal definition makes sense only in a combination with some renormalization procedure. Generally speaking the renormalization includes subtractions which can violate naive positivity bounds.

Working with the regularizations preserving the positivity of the norm in the Hilbert space of states (with all necessary comments concerning the color singlet sector, the insertion of the $`P`$ exponent between parton fields etc.) one seems to be on the safe ground, which gives an argument in favor of the validity of the positivity bounds at high normalization points. On the other hand, in the case of FPDs it is well known that only the cross sections associated with FPDs must be positive whereas the naive positivity bounds may be violated for FPDs at low normalization points in nonphysical renormalization schemes AFR-98 . One should keep in mind that starting from the general inequality (3) one can reproduce the standard positivity properties of FPDs. Therefore a violation of the positivity properties of FPDs at low normalization points would lead to the breakdown of the positivity bounds on GPDs (this can be directly seen in inequalities (27), (28) where the GPDs are constrained by FPDs).

In this paper we have shown that the validity of the positivity bounds on GPDs at high normalization points is compatible with the one-loop evolution. This self-consistency check is encouraging but certainly more serious analysis is needed in order to clarify the status of inequality (3) in the context of the renormalization.

## VII Conclusions

Using the impact parameter representation we have derived positivity bound (11) on GPDs. This inequality should hold for any function $`p_\sigma `$ so that actually we deal with an infinite set of inequalities. These positivity bounds impose certain constraints on models of GPDs used in the phenomenological analysis of hard exclusive processes.

Acknowledgements. I appreciate discussions with A. Belitsky, V. Braun, J.C. Collins, M. Diehl, L. Frankfurt, D.S. Hwang, X. Ji, M. Kirch, N. Kivel, L.N. Lipatov, A. Manashov, D. Müller, V.Yu. Petrov, M.V. Polyakov, A.V. Radyushkin, M. Strikman and O. Teryaev. This work was supported by DFG and BMBF.

## Appendix A Positivity properties of evolution kernels

In this Appendix we show how the positivity property (7) of the evolution kernel $`K_{\beta _1\beta _2}^{\alpha _1\alpha _2}(x_1,x_2;y_1,y_2)`$ can be seen in the direct diagrammatic calculation of these kernels. Generally speaking, the leading order evolution kernel gets contributions from several diagrams BFKL-85  but in the region $`0<x_k<y_k`$ (6) which is interesting for us this kernel is given by the single cut diagram of Fig. 2. This diagram leads to the following structure of the evolution kernel

$$K_{\beta _1\beta _2}^{\alpha _1\alpha _2}(x_1,x_2;y_1,y_2)=\underset{\gamma }{}\delta (y_1x_1y_2+x_2)$$

$$\times V_{\alpha _1\beta _1\gamma }(x_1,y_1)V_{\alpha _2\beta _2\gamma }^{}(x_2,y_2).$$
(29)
In the framework of Ref. BFKL-85  this contribution naturally appears from the pole integration over the light-like component of the momentum. In the case $`x_1,x_2>0`$ the poles of the propagators associated with $`x_1`$ and $`x_2`$ lie on the same side of the real axis. Shifting the integration contour to the opposite side one can get a nonzero contribution only from the pole of the propagator corresponding to the intermediate parton in the diagram of Fig. 2. Functions $`V_{\alpha _k\beta _k\gamma }(x_k,y_k)`$ correspond to the two vertices of the diagram of Fig. 2. The kinematical and normalization factors associated with the cut propagator are obviously positive (this positivity is explicit in the light-cone gauge) and can be included into factors $`V_{\alpha _k\beta _k\gamma }(x_k,y_k)`$. In the RHS of Eq. (29) one sums over the polarization $`\gamma `$ of the intermediate parton. Using the light-cone gauge and the techniques of Ref. BFKL-85  it is easy to see that the summation over $`\gamma `$ is restricted to the physical polarizations of the intermediate parton. Note that we can rewrite (29) in the form

$$K_{\beta _1\beta _2}^{\alpha _1\alpha _2}(x_1,x_2;y_1,y_2)$$

$$=\underset{0}{\overset{\mathrm{}}{}}𝑑z\underset{\gamma }{}W_{\alpha _1\beta _1\gamma ,z}(x_1,y_1)W_{\alpha _2\beta _2\gamma ,z}^{}(x_2,y_2)$$
(30)
where

$$W_{\alpha \beta \gamma ,z}(x,y)=\delta (yxz)V_{\alpha \beta \gamma }(x,y).$$
(31)
Obviously the form (30) of the evolution kernel $`K`$ automatically leads to the positivity property (7). Actually the general decomposition (30) is essentially equivalent to the inequality (7). It is decomposition (30) that is needed for the proof of the stability of the positivity bounds on GPDs with respect to the evolution to higher normalization points.

In order to illustrate the above general formulas with an explicit example let us consider the evolution kernel for the “helicity-independent” quark GPDs

$$\underset{\alpha }{}K_{\beta _1\beta _2}^{\alpha \alpha }(x_1,x_2;y_1,y_2)\delta _{\beta _1\beta _2}K(x_1,x_2;y_1,y_2).$$
(32)
In the region (6), $`0<x_k<y_k`$, the kernel $`K`$ is given by the following expression BFKL-85

$$K(x_1,x_2;y_1,y_2)=\frac{C_F}{2\pi ^2}$$

$$\times \frac{\delta (y_1x_1y_2+x_2)}{y_1x_1}\left(1+\frac{x_1x_2}{y_1y_2}\right).$$
(33)
Obviously this kernel has the form (29)

$$K(x_1,x_2;y_1,y_2)=\frac{C_F}{2\pi ^2}\delta (y_1x_1y_2+x_2)$$

$$\times \left(\frac{1}{\sqrt{y_1x_1}}\frac{1}{\sqrt{y_2x_2}}+\frac{x_1}{y_1\sqrt{y_1x_1}}\frac{x_2}{y_2\sqrt{y_2x_2}}\right).$$
(34)
This contribution comes from the diagram of Fig. 2 with a gluon playing the role of the intermediate parton. The two terms in the brackets in the RHS correspond to two possible polarizations of this gluon.

The above expression for the kernel $`K`$ is singular at the points $`y_k=x_k`$. The proper treatment of these singularities is described by the following expression for the convolution of the kernel $`K`$ with arbitrary functions $`h(x,y)`$ Ji-97 ; Radyushkin-97 ; BM-98

$$\underset{0}{\overset{\mathrm{}}{}}𝑑x_1\underset{x_1}{\overset{\mathrm{}}{}}𝑑y_1\underset{0}{\overset{\mathrm{}}{}}𝑑x_2\underset{x_2}{\overset{\mathrm{}}{}}𝑑y_2h(x_1,y_1)h^{}(x_1,y_1)$$

$$\times K(x_1,x_2;y_1,y_2)=\frac{C_F}{2\pi ^2}\underset{\epsilon 0}{lim}\underset{0}{\overset{\mathrm{}}{}}dx_1\underset{0}{\overset{\mathrm{}}{}}dx_2\{\underset{\epsilon }{\overset{\mathrm{}}{}}\frac{dz}{z}$$

$$\times \left[1+\frac{x_1x_2}{(x_1+z)(x_2+z)}\right]h(x_1,x_1+z)h^{}(x_2,x_2+z)$$

$$+(\frac{3}{2}\mathrm{ln}\frac{4x_1x_2}{\epsilon ^2})h(x_1,x_1)h^{}(x_2,x_2)\}.$$
(35)
The last term in the RHS proportional to $`h(x_1,x_1)h^{}(x_2,x_2)`$ corresponds to the virtual $`\delta (y_1x_1)\delta (y_2x_2)`$ contributions to the evolution kernel $`K`$. In the limit $`\epsilon 0`$ both integral and contact terms are divergent but these two divergences cancel each other.

## Appendix B Stability of positivity bounds under evolution

In this Appendix we present a detailed derivation of the stability of the positivity bounds (3) under the one-loop evolution (5) to higher normalization points. We shall use compact Dirac notation for the integral appearing in the LHS of inequality (3)

$$f|G|f\underset{i_1\alpha _1i_2\alpha _2}{}\frac{dP_1^+d^2P_1^{}dx_1}{2P_1^+(2\pi )^3}\frac{dP_2^+d^2P_2^{}dx_2}{2P_2^+(2\pi )^3}$$

$$\times \theta (x_1)\theta (x_2)f_{i_2\alpha _2}^{}(x_2,P_2)f_{i_1\alpha _1}(x_1,P_1)G_{P_1i_1,P_2i_2}^{\alpha _1\alpha _2}(x_1,x_2)$$
(36)
so that the positivity bound (3) can be written as follows

$$f|G|f0.$$
(37)
Let us assume that at some normalization point inequality (3) holds for all functions $`f_{i\alpha }(x,P)`$. Imagine that during the evolution to higher $`\mu `$ this inequality breaks down at some point $`\mu _0`$ for some function $`f^{(0)}`$ so that at $`\mu >\mu _0`$

$$f^{(0)}|G_\mu |f^{(0)}<0(\mu >\mu _0)$$
(38)
but at $`\mu <\mu _0`$ we still have for all functions $`f`$

$$f|G_\mu |f0(\mu \mu _0).$$
(39)
Then at point $`\mu _0`$

$$f^{(0)}|G_{\mu _0}|f^{(0)}=0.$$
(40)
The fate of this “matrix element” at $`\mu >\mu _0`$ is determined by the evolution equation (5)

$$\mu \frac{}{\mu }f^{(0)}|G_\mu |f^{(0)}=f^{(0)}|(KG_\mu )|f^{(0)}.$$
(41)
Here we use short notation $`KG_\mu `$ for the convolution of the evolution kernel $`K`$ with GPD $`G`$. If one could show that at the point $`\mu _0`$ we have

$$f^{(0)}|(KG_{\mu _0})|f^{(0)}0$$
(42)
then this would guarantee that the evolution from $`\mu _0`$ to higher $`\mu `$ would restore the positivity of $`f^{(0)}|G_\mu |f^{(0)}`$. This would invalidate our assumption (38). Thus in order to prove the stability of the positivity bounds on GPD under the evolution to higher normalization points we need the following

Statement. If at some normalization point $`\mu _0`$ for all functions $`f`$

$$f|G_{\mu _0}|f0$$
(43)
and for some $`f^{(0)}`$

$$f^{(0)}|G_{\mu _0}|f^{(0)}=0$$
(44)
then

$$f^{(0)}|(KG_{\mu _0})|f^{(0)}0.$$
(45)
To prove this statement we make use of the positivity property of the evolution kernels (7). As it is explained in Appendix A this positivity property is equivalent to the decomposition (30) of the evolution kernel. In our short Dirac-like notations Eq. (30) takes the following form

$$K=𝑑z\underset{\gamma }{}|W_{\gamma ,z}W_{\gamma ,z}|$$
(46)
Inserting representation (46) for $`K`$ into the LHS of inequality (45) we find

$$f^{(0)}|(KG_{\mu _0})|f^{(0)}$$

$$=𝑑z\underset{\gamma }{}f^{(0)}W_{\gamma ,z}|G_{\mu _0}|f^{(0)}W_{\gamma ,z}0.$$
(47)
The RHS is positive according to Eq. (43). This completes the derivation of the inequality (45).

The above proof of the stability of the positivity bounds (3) under the evolution (5) ignored the problem of the negative terms proportional to $`\delta (x_1y_1)\delta (x_2y_2)`$ which are present in the evolution kernel $`K`$. Because of these terms it is allowed to use inequality (7) only for functions $`h(x,y)`$ which vanish at $`x=y`$. As a result instead of Eq. (46) one has to work with the following representation for the evolution kernel

$$K=\underset{\epsilon +0}{lim}\left(K_\epsilon ^{(1)}+K_\epsilon ^{(2)}\right)$$
(48)
where the first term is given by the $`\epsilon `$ regularized integral decomposition (46)

$$K_\epsilon ^{(1)}(x_1,x_2,y_1,y_2)=\underset{\epsilon }{\overset{\mathrm{}}{}}𝑑z\underset{\gamma }{}W_{\gamma ,z}(x_1,y_1)W_{\gamma ,z}^{}(x_2,y_2)$$
(49)
and the second piece

$$K_\epsilon ^{(2)}(x_1,x_2,y_1,y_2)=C(\epsilon ,x_1,x_2)\delta (x_1y_1)\delta (x_2y_2)$$
(50)
comes from the virtual contribution to the evolution kernel which is proportional to $`\delta (x_1y_1)\delta (x_2y_2)`$ with a divergent coefficient $`C(\epsilon ,x_1,x_2)`$. This coefficient is regularized by the same small parameter $`\epsilon `$ so that

$$C(\epsilon ,x_1,x_2)=c_1+c_2\mathrm{ln}\epsilon +F_1(x_1)+F_2(x_2).$$
(51)
The singular contribution $`\mathrm{ln}\epsilon `$ compensates the divergence of the $`z`$ integral in the RHS of Eq. (49). Obviously representation (48) guarantees that $`K`$ obeys inequality (7) for arbitrary functions $`h(x,y)`$ vanishing at $`x=y`$. An example of the general structure (51) can be seen in the explicit expression (35) for the quark evolution kernel where the term $`\mathrm{ln}((4x_1x_2)/\epsilon ^2)`$ represents all terms appearing in the rhs of (51).

Now we have to show that the kernel $`K`$ (48) still has the property (45). For the first term in the RHS of Eq. (48) this can be done in the same way as in the inequality (47):

$$f^{(0)}|(K_\epsilon ^{(1)}G_{\mu _0})|f^{(0)}0.$$
(52)
Next let us show that

$$f^{(0)}|(K_\epsilon ^{(2)}G_{\mu _0})|f^{(0)}=0.$$
(53)
According to Eqs. (50), (51) we have

$$f^{(0)}|(K_\epsilon ^{(2)}G_{\mu _0})|f^{(0)}=f^{(0)}|[(c_1+c_2\mathrm{ln}\epsilon )G_{\mu _0}$$

$$+G_{\mu _0}F_1+F_2G_{\mu _0}]|f^{(0)}.$$
(54)
Note that $`|f^{(0)}`$ is associated with $`f^{(0)}(x_1,y_1)`$ whereas $`f^{(0)}|`$ corresponds to $`f^{(0)}(x_2,y_2)`$ in Eq. (36). This explains the ordering of “operators” $`F_1`$ and $`F_2`$ associated with functions $`F_1(x_1)`$ and $`F_2(x_2)`$ respectively in the RHS of Eq. (54). Taking into account Eq. (44) we get rid of the $`c_1+c_2\mathrm{ln}\epsilon `$ contribution:

$$f^{(0)}|(K_\epsilon ^{(2)}G_{\mu _0})|f^{(0)}$$

$$=f^{(0)}|(G_{\mu _0}F_1+F_2G_{\mu _0})|f^{(0)}.$$
(55)
Since property of $`G_{\mu _0}`$ (43) holds for any $`f`$ we conclude from Eq. (44) that

$$G_{\mu _0}|f^{(0)}=0,f^{(0)}|G_{\mu _0}=0.$$
(56)
This means that

$$f^{(0)}|(G_{\mu _0}F_1+F_2G_{\mu _0})|f^{(0)}=0.$$
(57)
Inserting this result into Eq. (55) we obtain Eq. (53). Combining Eqs. (52), (53) we derive inequality (45), which means that the virtual terms in the evolution kernel do not violate the stability of the positivity bounds on parton distributions under the evolution to higher normalization points.

## Appendix C Derivation of positivity bounds in impact parameter representation

In this Appendix we show how inequality (3) containing a multidimensional integral can be reduced to the relatively simple form (11).

Using representation (8) for GPD $`G`$ we find from inequality (3)

$$\underset{0}{\overset{\mathrm{}}{}}𝑑k_1^+\underset{0}{\overset{\mathrm{}}{}}\frac{dP_1^+}{P_1^+}d^2P_1^{}\underset{0}{\overset{\mathrm{}}{}}𝑑k_2^+\underset{0}{\overset{\mathrm{}}{}}\frac{dP_2^+}{P_2^+}d^2P_2^{}$$

$$\times \delta \left[(P_2^+k_2^+)(P_1^+k_1^+)\right]$$

$$\times F_{\sigma _2\sigma _1}[\frac{k_1^++k_2^+}{P_1^++P_2^+},\frac{P_1^+P_2^+}{P_1^++P_2^+},\frac{2(P_1^+P_2^{}P_2^+P_1^{})}{P_1^++P_2^+}]$$

$$\times f_{\sigma _2}^{}(k_2^+,P_2^+,P_2^{})f_{\sigma _1}(k_1^+,P_1^+,P_1^{})0.$$
(58)
Taking functions $`f_\sigma `$ in the factorized form

$$f_\sigma (k^+,P^+,P^{})=d_\sigma (k^+,P^+)\eta \left(\frac{P^{}}{P^+}\right)$$
(59)
we derive from inequality (58)

$$\underset{0}{\overset{\mathrm{}}{}}𝑑k_1^+\underset{0}{\overset{\mathrm{}}{}}𝑑P_1^+\underset{0}{\overset{\mathrm{}}{}}𝑑k_2^+\underset{0}{\overset{\mathrm{}}{}}𝑑P_2^+d_{\sigma _2}^{}(k_2^+,P_2^+)d_{\sigma _1}(k_1^+,P_1^+)$$

$$\times \stackrel{~}{F}_{\sigma _2\sigma _1}[\frac{k_1^++k_2^+}{P_1^++P_2^+},\frac{P_1^+P_2^+}{P_1^++P_2^+},\frac{P_1^++P_2^+}{2P_1^+P_2^+}b^{}]$$

$$\times \delta [(P_2^+k_2^+)(P_1^+k_1^+)]\left(\frac{P_1^++P_2^+}{2P_1^+P_2^+}\right)^20$$
(60)
where $`\stackrel{~}{F}_{\sigma _2\sigma _1}`$ is defined by Eq. (10). This inequality should hold for any value of $`b^{}`$ and for any function $`d_\sigma (k^+,P^+)`$. The next step is to take

$$d_\sigma (k^+,P^+)=g(P^+k^+)h_\sigma (P^+).$$
(61)
In the limit

$$\left|g(u)\right|^2\delta (u\nu )$$
(62)
we find from inequality (60)

$$\underset{1}{\overset{1}{}}𝑑\xi \underset{|\xi |}{\overset{1}{}}𝑑x\frac{1}{(1x)(1\xi ^2)^2}h_{\sigma _2}^{}\left(\nu \frac{1\xi }{1x}\right)$$

$$\times h_{\sigma _1}\left(\nu \frac{1+\xi }{1x}\right)\stackrel{~}{F}_{\sigma _2\sigma _1}(x,\xi ,\frac{(1x)b^{}}{\nu (1\xi ^2)})0.$$
(63)
Since function $`h_\sigma `$ and parameter $`b^{}`$ are arbitrary we can set $`\nu =1`$ without losing generality. Taking $`h_\sigma (z)=z^2p_\sigma (z^1)`$ we obtain inequality (11).