text:
# The Pinch Technique to All Orders

(August 21, 2002)

## Abstract

The generalization of the pinch technique to all orders in perturbation theory is presented. The effective Green’s functions constructed with this procedure are singled out in a unique way through the full exploitation of the underlying Becchi-Rouet-Stora-Tyutin symmetry. A simple all-order correspondence between the pinch technique and the background field method in the Feynman gauge is established.

preprint: FTUV-02-0821preprint: IFIC-02-19

It is well-known that, to any finite order, the conventional perturbative expansion gives rise to expressions for physical amplitudes which are endowed with crucial properties. $`S`$-matrix elements, for example, are independent of the gauge-fixing scheme and parameters chosen to quantize the theory, they are gauge-invariant (current conservation), they are unitary (conservation of probability), and well behaved at high energies. However, the above properties are in general not reflected by the individual off-shell Green’s functions, which are the building blocks of the aforementioned perturbative expansion. The latter depend on the gauge-fixing parameters in a complicated way , grow much faster than physical amplitudes at high energies, and display unphysical thresholds. Evidently, when combining unphysical Green’s functions to form a physical amplitude, subtle field-theoretical mechanisms are at work, which enforce non-trivial cancellations among them at any given order.

There are considerable conceptual and phenomenological advantages in reformulating the perturbative expansion in terms of off-shell Green’s functions which display manifestly the same properties as the physical amplitudes. To begin with, the sharp difference between observables and Green’s functions suggests a great deal of redundancy in the conventional diagrammatic formulation of gauge theories, in the sense that extensive underlying cancellations beg to be made manifest and be explicitly exploited as early within a calculation as possible. Implementing these cancellations at an early stage not only renders the book-keeping aspects more tractable Feng:1995vg , but allows for theoretically safe reorganizations or resummations of the perturbative series. For example, identifying and Dyson-resuming the correct sub-set of propagator-like corrections gives rise to physically meaningful Born-improved amplitudes Papavassiliou:1995fq . In addition, the generalization into a non-Abelian context of the characteristic properties of the QED effective charge, has a wide range of phenomenological applications Brodsky:2001wx . Finally, $`n`$-point functions free of unphysical artifacts could serve, at least in principle, as the new building blocks of manifestly gauge-invariant Schwinger-Dyson equations Cornwall:1982zr .

It would clearly be preferable to enforce the relevant cancellations already at the level of the functional path-integral defining the theory, and obtain directly from it the desired Green’s functions; this is however beyond our powers at the moment. On the other hand, there exists a diagrammatic method, called the pinch technique (PT) Cornwall:1982zr ; Cornwall:1989gv , which systematically exploits the symmetries built into physical observables, such as $`S`$-matrix elements, in order to construct off-shell sub-amplitudes that are kinematically akin to conventional Green’s functions, but, unlike the latter, are also endowed with desirable properties. The basic observation, which essentially defines the PT, is that there exists a fundamental cancellation between sets of diagrams with different kinematic properties, such as self-energies, vertices, and boxes. This cancellation is driven by the underlying Becchi-Rouet-Stora-Tyutin symmetry Becchi:1976nq , and is triggered when longitudinal momenta circulating inside vertex and box diagrams generate (by “pinching” out internal fermion lines) propagator-like terms. The latter are reassigned to conventional self-energy graphs in order to give rise to the aforementioned gauge-invariant effective Green’s functions. In its original one- Cornwall:1982zr ; Cornwall:1989gv  and two-loop Papavassiliou:2000az  application, the PT boils down to the study of the kinematic rearrangements produced into individual Feynman diagrams when elementary tree-level Ward identities (WIs) are triggered.

One of the most pressing questions in this context is whether one can extend the PT algorithm to all orders in perturbation theory, thus achieving the systematic construction of effective $`n`$-point functions displaying the aforementioned characteristic features. To accomplish this it is clear that one needs to abandon algebraic operations inside individual Feynman graphs, and resort to a more formal procedure. In this Letter we will show that the PT algorithm can be successfully generalized to all orders in perturbation theory, through the collective treatment of entire sets of diagrams. This is accomplished through the judicious use of the Slavnov-Taylor identity (STI) Slavnov:1972fg  satisfied by a special Green’s function, which serves as a common kernel to all higher order self-energy and vertex diagrams.

We will consider for concreteness the $`S`$-matrix element for the quark–anti-quark elastic scattering process $`q(r_1)\overline{q}(r_2)q(p_1)\overline{q}(p_2)`$ in QCD. We set $`q=r_2r_1=p_2p_1`$, with $`s=q^2`$ the square of the momentum transfer. The longitudinal momenta responsible for the aforementioned kinematical rearrangements stem either from the bare gluon propagators or from the pinching part $`\mathrm{\Gamma }_{\alpha \mu \nu }^\mathrm{P}(q,k_1,k_2)`$ appearing in the characteristic decomposition of the elementary tree-level three-gluon vertex $`\mathrm{\Gamma }_{\alpha \mu \nu }^{eab,[0]}=gf^{eab}\mathrm{\Gamma }_{\alpha \mu \nu }^{[0]}`$ into Cornwall:1982zr

$`\mathrm{\Gamma }_{\alpha \mu \nu }^{[0]}(q,k_1,k_2)`$ $`=`$ $`\mathrm{\Gamma }_{\alpha \mu \nu }^\mathrm{F}(q,k_1,k_2)+\mathrm{\Gamma }_{\alpha \mu \nu }^\mathrm{P}(q,k_1,k_2),`$
$`\mathrm{\Gamma }_{\alpha \mu \nu }^\mathrm{F}(q,k_1,k_2)`$ $`=`$ $`(k_1k_2)_\alpha g_{\mu \nu }+2q_\nu g_{\alpha \mu }2q_\mu g_{\alpha \nu },`$
$`\mathrm{\Gamma }_{\alpha \mu \nu }^\mathrm{P}(q,k_1,k_2)`$ $`=`$ $`k_{2\nu }g_{\alpha \mu }k_{1\mu }g_{\alpha \nu }.`$ (1)
The above decomposition is to be carried out to “external” three-gluon vertices only, i.e., the vertices where the physical momentum transfer $`q`$ is entering Papavassiliou:2000az . In what follows we work in the renormalizable Feynman gauge (RFG); this choice eliminates the longitudinal momenta from the bare propagators, and allows us to focus our attention on the all-order study of the longitudinal momenta originating from $`\mathrm{\Gamma }_{\alpha \mu \nu }^\mathrm{P}`$.

We will denote by $`𝒜`$ the subset of the graphs which will receive the action of the longitudinal momenta stemming from $`\mathrm{\Gamma }_{\alpha \mu \nu }^\mathrm{P}(q,k_1,k_2)`$ (see Fig.1). We have that

$$𝒜=ig^2\overline{u}(r_1)\frac{\lambda ^e}{2}\gamma _\alpha v(r_2)f^{eab}\mathrm{\Gamma }^{\mathrm{P},\alpha \mu \nu }(q,k_1,k_2)𝒯_{\mu \nu }^{ab}(k_1,k_2,p_1,p_2),$$
(2)
where $`\lambda ^e`$ are the Gell-Mann matrices, and $`𝒯_{\mu \nu }^{ab}`$ is the sub-amplitude $`g_\mu ^a(k_1)g_\nu ^b(k_2)q(p_1)\overline{q}(p_2)`$, with the gluons off-shell and the fermions on-shell; for the latter $`\overline{v}(p_2)S^1(p_2)|_{p/_2=m}=S^1(p_1)u(p_1)|_{p/_1=m}=0`$, where $`S(p)`$ is the (full) quark propagator. In terms of Green’s functions we have

$$𝒯_{\mu \nu }^{ab}=\overline{v}(p_2)[𝒞_{\rho \sigma }^{ab}(k_1,k_2,p_1,p_2)\mathrm{\Delta }_\mu ^\rho (k_1)\mathrm{\Delta }_\nu ^\sigma (k_2)]u(p_1).$$
(3)
Clearly, there is an equal contribution from the $`\mathrm{\Gamma }^\mathrm{P}`$ situated on the right hand-side of $`𝒯`$.

Let us focus on the STI satisfied by the amplitude $`𝒯_{\mu \nu }^{ab}`$; it reads

$`k_1^\mu C_{\mu \nu }^{ab}+k_{2\nu }G_1^{ab}igf^{bcd}Q_{1\nu }^{acd}gX_{1\nu }^{ab}+g\overline{X}_{1\nu }^{ab}=0,`$ (4)
where the Green’s function appearing in it are defined in Fig.2. The terms $`X_{1\nu }`$ and $`\overline{X}_{1\nu }`$ die on-shell, since they are missing one fermion propagator. Thus, we arrive at the on-shell STI for $`𝒯_{\mu \nu }^{ab}`$

$$k_1^\mu 𝒯_{\mu \nu }^{ab}=𝒮_{1\nu }^{ab},$$
(5)
with

$`𝒮_{1\nu }^{ab}`$ $`=`$ $`\overline{v}(p_2)[igf^{bcd}𝒬_{1\nu }^{acd}(k_1,k_2,p_1,p_2)D(k_1)`$ (6)
$``$ $`k_{2\nu }𝒢_1^{ab}(k_1,k_2,p_1,p_2)D(k_1)D(k_2)]u(p_1),`$
where we have defined $`Q_{1\nu }^{acd}(k_1,k_2,p_1,p_2)=𝒬_{1\nu }^{acd}(k_1,k_2,p_1,p_2)D(k_1)S(p_1)S(p_2)`$.

In perturbation theory both $`𝒯_{\mu \nu }^{ab}`$ and $`𝒮_{1\nu }^{ab}`$ are given by Feynman diagrams, which can be separated into distinct classes, depending on their kinematic dependence and their geometrical properties. Graphs which do not contain information about the kinematical details of the incoming test-quarks are self-energy graphs, whereas those which display a dependence on the test quarks are vertex graphs. The former depend only on the variable $`s`$, whereas the latter on both $`s`$ and the mass $`m`$ of the test quarks; equivalently, we will refer to them as $`s`$-channel or $`t`$-channel graphs, respectively. In addition to the $`s`$-$`t`$ decomposition, Feynman diagrams can be separated into one-particle irreducible (1PI) and one-particle reducible (1PR) ones. The crucial point is that the action of the momenta $`k_1^\mu `$ or $`k_2^\nu `$ on $`𝒯_{\mu \nu }^{ab}`$ does not respect, in general, the original $`s`$-$`t`$ and 1PI-1PR separation furnished by the Feynman diagrams (see third paper of Papavassiliou:1995fq ).

In other words, even though Eq.(5) holds for the entire amplitude, it is not true for the individual sub-amplitudes, i.e.,

$$k_1^\mu \left[𝒯_{\mu \nu }^{ab}\right]_{x,\mathrm{Y}}\left[𝒮_{1\nu }^{ab}\right]_{x,\mathrm{Y}},x=s,t;\mathrm{Y}=\mathrm{I},\mathrm{R},$$
(7)
where I (respectively R) indicates the one-particle irreducible (respectively reducible) parts of the amplitude involved. Evidently, whereas the characterization of graphs as propagator- and vertex-like is unambiguous in the absence of longitudinal momenta (e.g., in a scalar theory), their presence tends to mix propagator- and vertex-like graphs. Similarly, 1PR graphs are effectively converted into 1PI ones (the opposite cannot happen). The reason for the inequality of Eq.(7) are precisely the propagator-like terms, such as those encountered in the one- and two-loop calculations; they have the characteristic feature that, when depicted by means of Feynman diagrams contain unphysical vertices, i.e., vertices which do not exist in the original Lagrangian (Fig.3). All such diagrams cancel diagrammatically against each other. Thus, after the aforementioned rearrangements have taken place, for the $`t`$-channel irreducible part of the amplitude we will have the equality

$$\left[k_1^\mu 𝒯_{\mu \nu }^{ab}\right]_{t,\mathrm{I}}^{\mathrm{PT}}\left[𝒮_{1\nu }^{ab}\right]_{t,\mathrm{I}}.$$
(8)
Eq.(8) merits particular attention, because it is of central importance for the generalization of the PT to all orders. The superscript “PT” on the left hand-side denotes that the corresponding amplitude must be rearranged following the well-defined PT algorithm, as it has been explained in the literature Papavassiliou:2000az . In particular, one tracks down the rearrangments induced when the action of (virtual) longitudinal momenta ($`k`$) on the bare vertices of diagrams trigger elementary WIs. Eventually a WI of the form $`k_\mu \gamma ^\mu =S^1(k/+p/)S^1(p/)`$ will give rise to propagator-like parts, by removing (pinching out) the internal bare fermion propagator $`S(k/+p/)`$. Depending on the topology of the diagram under consideration this last WI may be activated immediately, or as the final outcome of a sequential triggering of intermediate WIs. We emphasize that, in order to preserve the special unitarity and analyticity properties of the PT Green’s functions, “internal” three-gluon vertices should not pinch, nor should one carry out sub-integrations Papavassiliou:2000az .

The non-trivial step for generalizing the PT to all orders is then the following: Instead of going through the arduous task of manipulating the left hand-side of Eq.(8), following the aforementioned rules, in order to determine the pinching parts and explicitly enforce their cancellation, use directly the right-hand side, which already contains the answer! Indeed, the right-hand side involves only conventional (ghost) Green’s functions, expressed in terms of normal Feynman rules, with no reference to unphysical vertices. That this must be so follows from the same PT rules mentioned above: due to the absence of external three-gluon vertices the right-hand side cannot be pinched further, i.e. its separation into propagator- and vertex-like graphs is unambiguous, since there is no possibility (without violating the PT rules) to obtain further mixing. Thus, the right-hand side of Eq.(8) serves as a practical definition of the PT to all orders.

After these observations, we proceed to the PT construction to all orders. Once the effective Green’s functions have been derived, they will be compared to the corresponding Green’s functions obtained in the Feynman gauge of the background field method (BFG for short) in order to establish whether the known correspondence persists to all orders; as we will see, this is indeed the case (for an extended list of related references see Binosi:2002ez ).

To begin with, it is immediate to recognize that in the RFG box diagrams of arbitrary order $`n`$, to be denoted by $`B^{[n]}`$, coincide with the PT boxes $`\widehat{B}^{[n]}`$, since all three-gluon vertices are “internal”, i.e., they do not provide longitudinal momenta. Thus, they coincide with the BFG boxes, $`\stackrel{~}{B}^{[n]}`$, i.e., $`\widehat{B}^{[n]}=B^{[n]}=\stackrel{~}{B}^{[n]}`$ for every $`n`$.

We then continue with the construction of the 1PI PT gluon-quark–anti-quark vertex $`\widehat{\mathrm{\Gamma }}_\alpha ^e`$. We start from the corresponding vertex in the RFG, to be denoted by $`\mathrm{\Gamma }_\alpha ^e`$, and focus only on the class of vertex diagrams containing an external bare three-gluon vertex; we will denote this subset by $`\mathrm{\Gamma }_{A^3,\alpha }^e`$ \[Fig.4(a)\]. All other types of graphs contributing to $`\mathrm{\Gamma }_\alpha ^e`$ are inert as far as the PT procedure is concerned, because they do not furnish pinching momenta Papavassiliou:2000az . The next step is to carry out the vertex decomposition of Eq.(1) to the external three-gluon vertex $`\mathrm{\Gamma }_{\alpha \mu \nu }^{eab,[0]}`$ appearing in $`\mathrm{\Gamma }_{A^3,\alpha }^e`$. This will result in the obvious separation $`\mathrm{\Gamma }_{A^3,\alpha }^e=\mathrm{\Gamma }_{A^3,\alpha }^{\mathrm{F},e}+\mathrm{\Gamma }_{A^3,\alpha }^{\mathrm{P},e}`$. The part $`\mathrm{\Gamma }_{A^3,\alpha }^{\mathrm{F},e}`$ is also inert, and will be left untouched. Thus, the only quantity to be further manipulated is $`\mathrm{\Gamma }_{A^3,\alpha }^{\mathrm{P},e}`$; it reads

$$\mathrm{\Gamma }_{A^3,\alpha }^{\mathrm{P},e}=gf^{eba}\left[(kq)^\mu g_\alpha ^\nu +k^\nu g_\alpha ^\mu \right]\left[𝒯_{\mu \nu }^{ab}\right]_{t,\mathrm{I}},$$
(9)
where $`\mu ^{2\epsilon }d^dk/(2\pi )^d`$, $`d=D2\epsilon `$, $`D`$ is the space-time dimension, and $`\mu `$ is the ’t Hooft mass. Following the discussion presented above, the pinching action amounts to the replacement $`k^\nu [𝒯_{\mu \nu }^{ab}]_{t,\mathrm{I}}[k^\nu 𝒯_{\mu \nu }^{ab}]_{t,\mathrm{I}}=\left[𝒮_{2\mu }^{ab}(k+q,k)\right]_{t,\mathrm{I}}`$ and similarly for the term coming from the momentum $`(kq)^\mu `$, i.e., $`[(kq)^\mu 𝒯_{\mu \nu }^{ab}]_{t,\mathrm{I}}=\left[𝒮_{1\nu }^{ab}(k+q,k)\right]_{t,\mathrm{I}}`$, or, equivalently,

$$\mathrm{\Gamma }_{A^3,\alpha }^{\mathrm{P},e}(q)gf^{eba}\left([𝒮_{2\alpha }^{ab}]_{t,\mathrm{I}}[𝒮_{1\alpha }^{ab}]_{t,\mathrm{I}}\right).$$
(10)
At this point the construction of the effective PT vertex $`\widehat{\mathrm{\Gamma }}_\alpha ^e`$ has been completed. The next important point is to study the connection between $`\widehat{\mathrm{\Gamma }}_\alpha ^e`$ and the vertex $`\stackrel{~}{\mathrm{\Gamma }}_\alpha ^e`$ in the BFG. To begin with, all “inert” terms contained in the original $`\mathrm{\Gamma }_\alpha ^e`$ carry over to the same sub-groups of graphs obtained in the BFG; most notably, the $`\mathrm{\Gamma }_{A^3,\alpha }^{\mathrm{F},e}`$ is precisely the $`\stackrel{~}{\mathrm{\Gamma }}_{\stackrel{~}{A}A^2,\alpha }^e`$ part of $`\stackrel{~}{\mathrm{\Gamma }}_\alpha ^e`$, where $`\stackrel{~}{A}`$ is the background gluon. The only exception are the ghost-diagrams contributing to $`\mathrm{\Gamma }_\alpha ^e`$ \[Fig.4(b)\]; the latter do not coincide with the corresponding ghost contributions in the BFG.

The important step is to recognize that the BFG ghost sector is provided precisely by combining the RFG ghosts with the right-hand side of Eq.(8). Specifically, one arrives at both the symmetric vertex $`\stackrel{~}{\mathrm{\Gamma }}_{\stackrel{~}{A}\overline{c}c}^e`$, characteristic of the BFG, as well as at the four-particle ghost vertex $`\stackrel{~}{\mathrm{\Gamma }}_{\stackrel{~}{A}A\overline{c}c}^e`$, which is totally absent in the conventional formalism \[Fig.4(c)\]. Indeed we find (omitting the spinors)

$`{\displaystyle \left[𝒮_{1\alpha }^{ab}\right]_{t,\mathrm{I}}}`$ $`=`$ $`{\displaystyle }D(k+q)\{k_\alpha \left[𝒢_1^{ab}(k+q,k)\right]_{t,\mathrm{I}}D(k)`$ (11)
$`+`$ $`igf^{bcd}\left[𝒬_{1\alpha }^{acd}(k+q,k)\right]_{t,\mathrm{I}}\}.`$
A similar equation, in which we have to trade the $`𝒢_1^{ab}`$ and $`𝒬_{1\alpha }^{acd}`$ Green’s functions for their Bose symmetric ones $`𝒢_2^{ab}`$ and $`𝒬_{2\alpha }^{acd}`$, holds for the $`𝒮_{2\alpha }`$ term. It is then easy to show that

$`\stackrel{~}{\mathrm{\Gamma }}_{\stackrel{~}{A}\overline{c}c,\alpha }^e(q)`$ $``$ $`\mathrm{\Gamma }_{A\overline{c}c,\alpha }^e(q)+gf^{eba}{\displaystyle }\{k_\alpha \left[𝒢_1^{ab}(k+q,k)\right]_{t,\mathrm{I}}`$
$`+`$ $`(kq)_\alpha \left[𝒢_2^{ab}(k+q,k)\right]_{t,\mathrm{I}}\}D(k+q)D(k),`$
$`\stackrel{~}{\mathrm{\Gamma }}_{\stackrel{~}{A}A\overline{c}c,\alpha }^e(q)`$ $``$ $`ig^2f^{eba}{\displaystyle }\{f^{acd}\left[𝒬_{2\alpha }^{cdb}(k+q,k)\right]_{t,\mathrm{I}}D(k)`$ (12)
$``$ $`f^{bcd}\left[𝒬_{1\alpha }^{acd}(k+q,k)\right]_{t,\mathrm{I}}D(k+q)\}.`$
This concludes the proof that $`\widehat{\mathrm{\Gamma }}_\alpha ^e\stackrel{~}{\mathrm{\Gamma }}_\alpha ^e`$. We emphasize that the sole ingredient in the above construction has been the STI of Eq.(5); in particular, at no point have we employed a priori the background formalism. Instead, its special ghost sector has arisen dynamically, once the PT rearrangement has taken place.

The final step is to construct the (all orders) PT gluon self-energy $`\widehat{\mathrm{\Pi }}_{\mu \nu }^{ab}`$. Notice that at this point one would expect that it too coincides with the BFG gluon self-energy $`\stackrel{~}{\mathrm{\Pi }}_{\mu \nu }^{ab}`$, since both the boxes as well as the vertex do coincide with the corresponding quantities in BFG, and the $`S`$-matrix is unique (renormalization may be carried out order-by-order without any complications, see second paper in Papavassiliou:2000az ). We will carry out a proof based on the strong induction principle, which states that a given predicate $`P(n)`$ on $``$ is true $`n`$, if $`P(k)`$ is true whenever $`P(j)`$ is true $`j`$ with $`j<k`$. We will use a schematic notation, suppressing Lorentz, color, and momentum indices. At one- and two-loop, we know that the result is true Cornwall:1982zr ; Papavassiliou:2000az . Assuming then that the PT construction has been successfully carried out up to the order $`n1`$, we will show that the PT gluon self-energy is equal to the BFG gluon self-energy at order $`n`$, hence proving that this equality holds true at any given $`n`$. From the inductive hypothesis, we know that $`\widehat{\mathrm{\Pi }}^{[\mathrm{}]}\stackrel{~}{\mathrm{\Pi }}^{[\mathrm{}]}`$, $`\widehat{\mathrm{\Gamma }}^{[\mathrm{}]}\stackrel{~}{\mathrm{\Gamma }}^{[\mathrm{}]}`$, and $`\widehat{B}^{[\mathrm{}]}\stackrel{~}{B}^{[\mathrm{}]}B^{[\mathrm{}]}`$, with $`\mathrm{}=1,\mathrm{},n1`$. Now, the $`S`$-matrix element of order $`n`$, to be denoted as $`S^{[n]}`$, assumes the form $`S^{[n]}=\left\{\mathrm{\Gamma }\mathrm{\Delta }\mathrm{\Gamma }\right\}^{[n]}+B^{[n]}`$. Moreover, since it is unique, regardless if it is written in the RFG, in the BFG, as well as before and after the PT rearrangement, we have that $`S^{[n]}\widehat{S}^{[n]}\stackrel{~}{S}^{[n]}`$. Using then the fact that $`\widehat{B}^{[\mathrm{}]}\stackrel{~}{B}^{[\mathrm{}]}`$ holds true even when $`\mathrm{}=n`$, we find that $`\left\{\mathrm{\Gamma }\mathrm{\Delta }\mathrm{\Gamma }\right\}^{[n]}\{\widehat{\mathrm{\Gamma }}\widehat{\mathrm{\Delta }}\widehat{\mathrm{\Gamma }}\}^{[n]}\{\stackrel{~}{\mathrm{\Gamma }}\stackrel{~}{\mathrm{\Delta }}\stackrel{~}{\mathrm{\Gamma }}\}^{[n]}`$. These amplitudes can then be split into 1PR and 1PI parts; in particular, the 1PR part after the PT rearrangement coincides with the 1PR part written in the BFG, since $`\left\{\mathrm{\Gamma }\mathrm{\Delta }\mathrm{\Gamma }\right\}_\mathrm{R}^{[n]}=\mathrm{\Gamma }^{[n_1]}\mathrm{\Delta }^{[n_2]}\mathrm{\Gamma }^{[n_3]}`$ with $`n_1,n_2,n_3<n`$, and $`n_1+n_2+n_3=n`$. This implies in turn the equivalence of the 1PI parts, i.e.,

$$\left(\widehat{\mathrm{\Gamma }}^{[n]}\stackrel{~}{\mathrm{\Gamma }}^{[n]}\right)\mathrm{\Delta }^{[0]}\mathrm{\Gamma }^{[0]}+\mathrm{\Gamma }^{[0]}\mathrm{\Delta }^{[0]}\left(\widehat{\mathrm{\Gamma }}^{[n]}\stackrel{~}{\mathrm{\Gamma }}^{[n]}\right)+\mathrm{\Gamma }^{[0]}\mathrm{\Delta }^{[0]}\left(\widehat{\mathrm{\Pi }}^{[n]}\stackrel{~}{\mathrm{\Pi }}^{[n]}\right)\mathrm{\Delta }^{[0]}\mathrm{\Gamma }^{[0]}0.$$
(13)
At this point, by means of the explicit construction presented for the vertex, we have that $`\widehat{\mathrm{\Gamma }}^{[n]}\stackrel{~}{\mathrm{\Gamma }}^{[n]}`$, so that one immediately gets $`\widehat{\mathrm{\Pi }}^{[n]}\stackrel{~}{\mathrm{\Pi }}^{[n]}`$. Hence, by strong induction, the above relation is true for any given perturbative order $`n`$, i.e., we have $`\widehat{\mathrm{\Pi }}_{\mu \nu }^{ab}\stackrel{~}{\mathrm{\Pi }}_{\mu \nu }^{ab}`$, q.e.d.

In conclusion, we have shown that the use of the STI (5) satisfied by the special Green’s function (3), allows for the generalization of the PT procedure to all orders. It would be interesting to further explore the physical meaning of the $`n`$-point functions obtained Bernabeu:2002nw , and establish possible connections with related formalisms Rebhan:1984bg .

Acknowledgments: This work has been supported by the CICYT Grants AEN-99/0692 and BFM2001-0262. D. B. thanks the Theoretical Physics Department of the University of Trento, where part of this work has been carried out.


format: text

text:
# Exotic shapes in 40Ca and 36Ar studied with Antisymmetrized Molecular Dynamics

## I Introduction

Owing to the experiments of gamma-ray measurements, many excited bands in the nuclei near 40Ca have been recently observed . The rotational bands with large moments of inertia are hot subjects relating with super deformations in $`sd`$-shell nuclei. In 40Ca, there exist many low-lying bands, which imply the coexistence of various shapes. It is interesting problem why the shape coexistence occurs even in the doubly magic nucleus. The shape coexistence problem and the mechanism of deformations in these $`sd`$-shell nuclei is one of the attractive subjects.

Our interests are possible exotic shapes and cluster aspects in the deformed excited states. It has been already known that clustering is one of the essential features in light nuclei . Exotic shapes are formed due to cluster structures in some light nuclei, For example, a famous cluster structure in 12C is three $`\alpha `$ clusters which form a triangle shape. Parity doublets of the rotational bands in 20Ne are described by a parity asymmetry shape with a 16O+$`\alpha `$ cluster structure. It is a long problem whether or not cluster aspects appear in heavier nuclei. Many of the familier cluster structures are those in light stable nuclei. In the recent experimental and theoretical research , cluster structures were suggested in light unstable nuclei such as 10Be and 12Be, where new type clusters such as 6He and 8He are proposed. Our aim in this paper is systematic study of the $`sd`$-shell nuclei near 40Ca to find possible exotic states with cluster aspects.

For the systematic study of ground and excited states of $`sd`$-shell nuclei, theoretical difficulties exist in the coexistence of cluster and mean-field aspects. Needless to say, cluster aspects are important in light nuclei, while the mean-field is an essential feature in heavy nuclei. Since the $`sd`$-shell is an inter-mediate mass-number region, where both aspects should be significant, we need a miscroscopic method beyond the traditional clsuter models. Namely, we adopt a method of antisymmetrized molecular dynamics(AMD) . The method AMD was applied for light nuclei and has been proved to be a powerful approach to describe the ground and excited states of unstable nuclei as well as stable nuclei. In the pioneering study with the AMD method, the structure changes on the ylast line of 20Ne are explained in terms of alteration of the cluster structure . In the recent studies of light unstable nuclei with AMD, a new cluster aspect in unstable nuclei was suggested that cluster cores develop in the deformed neutron mean-field in neutron-rich nuclei.

In the present work, we adopt a new effective force with a finite-rage three-body term, with which the binding energies and the radii of 4He, 16O and 40Ca are systematically reproduced, because the usual effective forces used in the AMD study, such as Minnesota forces, Volkov forces and Modified Volkov forces, have serious problems in reproducing these basic properties in the wide mass-number region.

In this paper, the structures of the ground and excited states of 40Ca and 36Ar are studied by the AMD method. We investigate the shape coexistence problem of these nuclei. The mechanism of deformations are discussed, focusing on cluster aspects. We also introduce the AMD study of neuton-rich $`sd`$-shell nuclei near 32Mg.

## II Formulation

In this section, the formulation of AMD for the nuclear structure study of ground and excited states is briefly explained. For more detailed descriptions of the AMD framework, the reader is referred to Refs .

The wave function of a system is written by AMD wave functions,

$$\mathrm{\Phi }=c\mathrm{\Phi }_{AMD}+c^{}\mathrm{\Phi }_{AMD}^{}+\mathrm{}.$$
(1)
An AMD wave function of a nucleus with a mass number $`A`$ is a Slater determinant of Gaussian wave packets;

$`\mathrm{\Phi }_{AMD}(𝐙)={\displaystyle \frac{1}{\sqrt{A!}}}𝒜\{\phi _1,\phi _2,\mathrm{},\phi _A\},`$ (2)
$`\phi _i=\varphi _{𝐗_i}\chi _{\xi _i}\tau _i:\{\begin{array}{c}\varphi _{𝐗_i}(𝐫_j)\mathrm{exp}\left[\nu \left(𝐫_j\frac{𝐗_i}{\sqrt{\nu }}\right)^2\right],\hfill \\ \chi _{\xi _i}=\left(\begin{array}{c}\frac{1}{2}+\xi _i\hfill \\ \frac{1}{2}\xi _i\hfill \end{array}\right),\hfill \end{array}`$ (7)
where the $`i`$th single-particle wave function $`\phi _i`$ is a product of the spatial wave function $`\varphi _{\mathrm{X}_i}`$, the intrinsic spin function $`\chi _{\xi _i}`$ and the iso-spin function $`\tau _i`$. The spatial part $`\varphi _{\mathrm{X}_i}`$ is presented by variational complex parameters $`X_{1i}`$, $`X_{2i}`$, $`X_{3i}`$. $`\chi _{\xi _i}`$ is the intrinsic spin function defined by $`\xi _i`$, and $`\tau _i`$ is the iso-spin function which is fixed to be up(proton) or down(neutron) in the present calculations. The values $`𝐙\{X_{ni},\xi _i\}(n=1,2,3\text{and }i=1,\mathrm{},A)`$ are the variational parameters which express an AMD wave function.

In order to obtain the wave functions of ground and excited states, we perform a generator coordinate method in the framework of AMD. First we vary the energy of the parity eigen state projected from an AMD wave function under a constraint that the total oscillator quanta must equal to a given number $`𝒩`$ as $`aa^{}=𝒩`$. The energy variation is numerically calculated by a frictional cooling method . An energy curve is obtained as a function of the coordinate $`𝒩`$ in the constraint. In the second step, the spin-parity eigen states projected from the obtained AMD wave functions are superposed by diagonalizing Hamiltonian and Norm matrices, $`P_{MK^{}}^{J\pm }\mathrm{\Phi }_{AMD}(𝐙_i)|H|P_{MK^{\prime \prime }}^{J\pm }\mathrm{\Phi }_{AMD}(𝐙_j)`$ and $`P_{MK^{}}^{J\pm }\mathrm{\Phi }_{AMD}(𝐙_i)|P_{MK^{\prime \prime }}^{J\pm }\mathrm{\Phi }_{AMD}(𝐙_j)`$, We call the present calculations as VBP(variation before projection) because the variation is performed before the spin projection.

## III Interactions

The ordinary effective forces with no three-body term such as the Volkov force or with a zero-range three-body term like the MV1 force are not appropriate to describe the binding energies and radii of nuclei covering wide mass number region from $`\alpha `$ to 40Ca. Hence we use an interaction containing a finite-range three-body term in the present calculations. The central part of the interaction is explained by a conbination of the two-body and three-body terms. The interaction parameters used in the present paper are as follows,

$`V_{central}={\displaystyle \underset{i<j}{}}V^{(2)}+{\displaystyle \underset{i<j<k}{}}V^{(3)},`$ (8)
$`V^{(2)}=(1m+bP_\sigma hP_\tau mP_\sigma P_\tau )\left\{V_a\mathrm{exp}[({\displaystyle \frac{r_{12}}{r_a}})^2]V_b\mathrm{exp}[({\displaystyle \frac{r_{12}}{r_b}})^2]\right\}`$ (9)
$`+V_c\mathrm{exp}[({\displaystyle \frac{r_{12}}{r_c}})^2],`$ (10)
$`V^{(3)}=V_d\mathrm{exp}[d(r_{12}^2+r_{23}^2+r_{31}^2)^2],`$ (11)
$`V_a=198.34\mathrm{MeV},V_b=300.86\mathrm{MeV},V_c=22.5\mathrm{MeV},`$ (12)
$`r_a=1.2\mathrm{fm},r_b=0.7\mathrm{fm},r_c=0.9\mathrm{fm},V_d=600\mathrm{MeV},d=0.8\mathrm{fm}^2`$ (13)
$`m=0.193,b=0.185,h=0.37,`$ (14)
where the width and strength parameters are chosen so as to reproduce reasonably the sizes of $`\alpha `$ and 40Ca, and the binding energies of $`\alpha `$, 16O and 40Ca. In choosing parameters, $`\alpha `$+$`\alpha `$ phase shifts and the saturation property of symmetric nuclear matter have been also taken into consideration. The adopted interaction is a sum of the central force, the G3RS spin-orbit force  with the strength $`u_{ls}=2500`$ MeV and the Coulomb force.

## IV Results

The excited deformed bands in the 40Ca and 36Ar have been recently observed in the gamma-ray measurements , which reveal the existence of many low-lying bands in these nuclei. In the experimental level scheme of 40Ca, the rotational band with a large moment of inertia (No.1 in Fig.1) is a hot topic relating with a possible super deformation.

We study the ground and excited states of 40Ca and 36Ar, with the method of a generator coordinate method in the framework of AMD with respect to the total oscillator quanta. In the energy curves obtained after spin-projection as a function of the oscillator quanta, we can not necessarily find local minima except for the absolute minimum. However, after diagonalization of the Hamiltonian and Norm matrices, the energy curves for the excited states have other local minima due to the conditions orthogonal to the lower states. As a result, we obtain many low-liying rotational bands (Fig.2), each of which is dominated by the spin-parity projected states from an intrinsic AMD wave function.

By analyzing deformations of these dominant intrinsic states, it is found that various shapes coexist in 40Ca. The interesting point is that each of the intrinsic states in 40Ca is well described by either of $`0p0h`$, $`4p4h`$ or $`8p8h`$ configurations. First we obtain the normal spherical ground state with the doubly closed $`sd`$-shell configuration. With the increase of the oscillator quanta, the deformed excited states appear in the low-energy region. In the region concerning $`4p4h`$ states, there exist opposite shapes as two local minimum states. One is an oblate shape with the deformation parameter $`\beta =0.2`$, while the other is a prolate deformation with $`\beta =0.2`$. The oblate state provides the excited $`K^\pi =0^+`$ band (2), and the bands (3) and (4) with $`K^\pi =0^+,2^+`$ are dominated by the states projected from the prolate intrinsic state. The reason for the higher $`0^+`$ state than the $`2^+`$ state in the $`K^\pi =0^+`$ band (3) is the mixing with the ground state. With a larger oscillator quanta, we find another largely prolate deformation which is dominated by $`8p8h`$ configurations. The energy levels of the band (5) obtained from this large deformation indicate a large moment of inertia of this rotational band. It is very surprising that the excited band with such highly excited configurations as $`8p8h`$ exists in low-energy region. We identify this band as the experimentally observed super-deformed band which starts from the $`0^+`$ state at 5.2 MeV.

As shown in Table I, the theoretical values of the $`E2`$ transition strength in the band (5) are extremely large comparing with those in other excited band. According to a quantitative comparison of the strength in the band (5) with the experimental data, the theoretical $`E2`$ strength for $`4^+2^+`$ is smaller than the corresponding experimental data $`B(E2;6.54\mathrm{MeV}5.63\mathrm{MeV})=100(W.u.)`$. Therefore, it is conjectured that the intrinsic deformation $`\beta =0.4`$ of the band (5) in the present results must be still small compared with the value $`\beta =0.6`$ for a typical super deformation. The underestimation of the deformations in the present theory will be improved by the superposition along generator coordinates and also the extension of the single-particle wave functions.

As mentioned above, the point in the present results of 40Ca is the coexistence of various shapes: the spherical state, the oblate shape, the normal prolate deformation and the large prolate deformation. By analysing the intrinsic structures of these deformed states, exotic shapes and cluster aspects are revealed. The density distributions of the dominant intrinsic states of the excited bands are presented in Fig.3. As shown in Fig.3, the oblate state (2) has a hexagon structure like a three-leaf clover which consists of three 12C cores surrounding an $`\alpha `$ at the center. In the super deformation, the system has a parity-asymmetry shape like a pear because of a 28Si+12C-like clustering. We focus the characteristics of the 28Si+12C-like cluster structures in the super deformation. The first point is that the clusters are not weak coupling but are strong coupling ones. It is interesting that the 12C cluster consists of four nucleons in the $`sd`$-shell and eight nucleons in the $`pf`$-shell. Namely, it is a inter-shell cluster which lies over different shells. Because of the strong coupling feature, it is expected that the deformed mean-field effect can be important as well as cluster aspects in the state. In second, if such a cluster structure indeed develops, it forms a parity asymmetric deformation, which may provide the parity doublets made of positive and negative-parity bands. Thirdly, it is reasonable that both clusters are the sub-shell closed nuclei. In other words, the stability of these cluster cores is based on the shell effects of the sub-shell closure. For further analysis of the cluster aspects, we need advanced calculations where the residual effects such as the coupling of other cluster channels, the tail of inter-cluster motion, and the core deformations are more carefully taken into consideration.

It should be pointed that the present results imply the possible exotic shapes due to the cluster effects such as 12C cores, though there are few experimental evidences of the existence of clusters in 40Ca. The investigation of the level structures of the negative-parity states is one of the probes, because the parity-doublets of positive and negative-parity bands support a cluster structure as the origin of a parity-asymmetry intrinsic state. The present calculations predict a $`K^\pi =0^{}`$ band with the band-head $`1^{}`$ state at 11.7 MeV , which is about 3 MeV higher than theoretical band-head energy of the positive-parity super deformation. The deformation of the dominant intrinsic state is rather large as $`\beta =0.3`$ and the intrinsic structure is similar to that of the super deformation in the positive-parity states. Therefore we consider this negative-parity band as a candidate of the parity doublet. Strictly speaking, in negative parity states, the rotational levels fragments because the state mixing between bands is rather strong. As a result, the intra-band $`E2`$ transitions in the $`K^\pi =0^{}`$ band are not so strong as those in the positive-parity band of the super deformation.

In the theoretical results of 36Ar, we find excited bands with a large prolate deformation $`\beta =0.3`$ which are dominated by $`4p8h`$ configurations. This band starts from a 4 MeV excitation energy, and well corresponds to the newly observed rotational band with a large moment of inertia in the experimental data .

At the end of this section, I briefly report the AMD study on neutron-rich $`sd`$-shell nuclei near 32Mg calculated by one of the authors(M. K.). The structures of the ground and excited states of 32Mg are hot subjects concerning the island of inversion of a neutron magic number $`N=20`$. The ground and excited states of Mg isotopes are calculated by an extended version of AMD with the Gogny force. The calculations well reproduce the energy levels of 32Me and the features of the island of inversion. In the neutron-rich $`sd`$-shell nuclei, the deformed neutron mean-field plays an important role. The unique predictions in this work are the low-lying negative parity states in 32Mg. Since the nuclei 32Mg and 40Ca have the same neutron number, it is interesting future problem to compare the intruder states of 32Mg with the excited states of 40Ca.

## V Summary

The structures of the ground and excited states of $`sd`$-shell nuclei: 40Ca and 36Ar were studied with the method of antisymmetrized molecular dynamics. The recently observed excited bands were well described with the deformed intrinsic states in the present calculations. It was found that various kinds of shape coexist in the low-energy region. In 40Ca, it was suggested in the present results that the oblate shape, the normal prolate deformation and the large prolate deformation coexist in the low-lying excited bands besides the spherical ground state. The oblate, normal prolate, and the large prolate states are explained by the excited configurations with $`4p4h`$, $`4p4h`$, and $`8p8h`$, respectively. It is a mystery that the rotational band with the highly excited configurations($`8p8h`$) starts from the low-energy region. The exotic shapes due to the cluster effects were implied in the intrinsic structures of these deformed states. The present results suggested possible cluster aspects in heavy $`sd`$-shell nuclei.

###### Acknowledgements.

Authors would like to thank Professors A. Tohsaki, Y. Akaishi and K. Ikeda for helpful discussions and comments. They are also thankful to Dr. N. Itagaki and E. Ideguchi for many discussions. This work was partially performed in the “Research Project for Study of Unstable Nuclei from Nuclear Cluster Aspects” sponsored by Institute of Physical and Chemical Research (RIKEN). The computational calculations of this work are supported by the Supercomputer Project No.58, No.70 of High Energy Accelerator Research Organization(KEK), and Research Center for Nuclear Physics(RCNP) in Osaka University.


text:
# A stochastic model for heart rate fluctuations

## I Introduction

Various methods and models have been used in attempts to characterize the dynamics of the heart rate control mechanism. For short time periods and under stationary conditions there are successful models of heart rate and blood pressure regulation Seidel (1998); ten Voorde (1992), but the characterization of long-term behavior has been a very difficult problem. Some models have been introduced in order to explain long-term fluctuations, but usually they can only describe well-controlled in vitro experiments, or the models depend on large number of parameters, which cannot be easily determined from experimental data Glass and Mackey (1988). Furthermore, these models can predict only global statistical features like scaling properties of power spectrum and correlations Ivanov et al. (1998) and tell us very little about the details of the time evolution.

Many features can be extracted from long time series of heart rate measurements, quantities like entropy measures Bettermann and van Leeuwen (1998); Kaspar and Schuster (1987); Pincus and Goldberger (1995); Pincus (1995); Rezek and Roberts (1998); Richman and Moorman (2000); Zhang and Roy (1999) correlation dimension Farmer et al. (1983); Fell et al. (1996); Grassberger and Procaccia (1983); Kantz and Schreiber (1995); Mayer-Kress et al. (1988); Yum et al. (1999), detrended fluctuations Peng et al. (1993, 1995); Iyengar et al. (1996), fractal dimensions Bassingthwaighte and Raymond (1995); Chau et al. (1993); Gough (1993); Zhang and Roy (1999), spectrum power-law exponents Bigger et al. (1996); Iyengar et al. (1996)) and symbolic dynamics complexity Palazzolo et al. (1998); Voss et al. (1995, 1996), but these are all purely statistical characterizations and as such cannot provide us a mathematical model of heart rate dynamics, not even a simple one. However, some of these statistical methods do characterize the complexity of the dynamics underlying the time series Kuusela et al. (2002), or are directly related to their fractal or chaotic features. Mathematical analysis of many physiological rhythms, including long-term heart rate fluctuations, has revealed that they are generated by processes which must be nonlinear, since linear systems can not produce such a complex behaviour Glass (2001). Nonlinear purely deterministic models can display chaotic dynamics and generate apparently unpredictable oscillations, but in practice it has not yet been possible to extract such models from real noisy experimental data. It is also possible that the underlying system is stochastic, i.e., the time evolution of the system is subject to a noise source. \[This kind of dynamical noise is different from measurement noise, which is mostly generated in the experimental apparatus.\] In any case, there is increasing evidence that noise, originated either from the system itself or as a reflection of external influences, is actually an integral part of the dynamics of biological systems Collins et al. (1996); Hidaka et al. (2000); Mar et al. (1999).

A typical R-R interval recording is shown in Fig. 1. The time series is generated by recording a 24-hour electrocardiogram and detecting the R-peak from each heartbeat, the R-R interval is the time difference between two consecutive R-peaks. In the upper panel of Fig. 1 we have the R-R interval time series for 6 hours. We can see sections where the oscillations are rather regular but there are also abrupt changes. In the lower panel of Fig. 1 we have zoomed into a part of the time series, about $`50`$ minutes, and also on this time range we can see apparently random oscillations with rapid changes.

It is well known that most short-time fluctuations of heart-rate are generated by respiration (periods typically in the couple of seconds range) and blood pressure regulation (so called Meyer waves with periods of about 10 seconds Guyton and Hall (1996)). In the following we are not interested in these fast rhythms (which can be analyzed quite well using linear or semi-linear models) but rather in time scales from minutes to hours. We will show that in this time range the dynamics of the heart-rate fluctuations can be well described by a one-dimensional Langevin type difference equation. This equation contains a deterministic part and additive Gaussian noise, and we have found that it works well when the delay parameter in the equation is in the range of 2–20 minutes.

## II The Model

An important and wide class of dynamic systems can be described by the Langevin differential equation van Kampen (1981); Risken (1984)

$$\frac{dX(t)}{dt}=g(X(t),t)+h(X(t),t)\mathrm{\Gamma }(t).$$
(1)
Here $`X(t)`$ represent the state of the system at time $`t`$, the function $`g`$ gives the nonlinear deterministic change, and the last term $`h`$ is the amplitude of the stochastic contribution and $`\mathrm{\Gamma }(t)`$ stands for uncorrelated white noise with vanishing mean. These kinds of stochastic differential equations always need an interpretation rule for the noise term, normally one uses the Ito interpretation Ito (1950). In general the functions $`g`$ and $`h`$ could depend explicitly on time $`t`$. The equation (1) can be easily generalized to higher dimensions. We will now show that long-term behaviour of heart rate can be modeled using a difference version of the Langevin equation van Kampen (1981)

$$X(t+\tau )=X(t)+g(X(t);\tau )+h(X(t);\tau )\mathrm{\Gamma }(t).$$
(2)
Here $`X(t)`$ again represents the state of the system, in this case the R-R interval, at time $`t`$, and $`\tau `$ is the time delay. If arbitrary small delays $`\tau `$ are possible then one can take the limit $`\tau 0`$ and get the differential equation (1) \[if the $`\tau `$ dependence is given by $`g(X(t);\tau )\tau g(X(t))`$\], but in the present case it will turn out that there is a minimum $`\tau `$ for which the model (2) seems to be valid. We assume that $`g`$ and $`h`$ do not have explicit time dependence but they may depend on the delay $`\tau `$. It is convenient to extract the term $`X(t)`$ in the deterministic part, as is done in (2), then a nonzero $`g(X(t);\tau )`$ stands for changes in the state of the system. An essential feature of models of the above type is that for time evolution we only need to know the state at one given moment and not its evolution in the past, i.e., they are Markovian Hänggi and Thomas (1982); van Kampen (1981).

The computational problem is now to determine the functions $`g`$ and $`h`$ from measured time series and to verify that the description using (2) is accurate. The principle of the method is very simple Gradišek et al. (2000); Siegert et al. (1998): at every time $`t_i`$ when the trajectory of the system meets an arbitrary but fixed point $`x`$ in state space, we look at the future state of the system at time $`t_i+\tau `$. The set of these future values (for a chosen $`x`$ and $`\tau `$) has a distribution in the state space and from this distribution we can determine the deterministic part $`g(x)`$ and the stochastic part $`h(x)`$, see Fig. 2A. In practice we first divide the range of the dynamical variable $`X`$ into equal boxes. By scanning the whole measured time series we check when $`X`$ is inside a given box $`x`$, i.e., $`|X(t_i)x|\mathrm{\Delta }x`$ , where $`x`$ is the middle value of the box and $`\mathrm{\Delta }x`$ is the half width of the box. When $`X`$ is found on the box, we look at the future value of the variable, $`X(t_i+\tau )`$, where $`\tau `$ is the fixed delay parameter. Since the trajectory of the system passes each box several times, we can calculate the distribution of the future values $`X(t_i+\tau )`$ for each box $`x`$. If we assume that the noise is Gaussian, we can fit a Gaussian function on each distribution, and as a result we get the mean and the deviation parameters for each $`x`$; the mean of this distribution is equal to $`x+g(x)`$ and the deviation is equal to $`h(x)`$ Friedrich et al. (2000); Timmer (2000). A typical case is given in Fig. 2B, and it shows that the distribution is actually very well described by Gaussian noise (the correlation is better than 0.95; the correlation is calculated as $`\sqrt{1S_{res}/S_{tot}}`$, where $`S_{res}`$ is the sum of the squared residuals and $`S_{tot}`$ is the variance). From the given data we can in this way determine the functions $`g(X)`$ and $`h(X)`$ needed in the stochastic model (2). It should be noted that we can calculate only the absolute value of $`h(X)`$ since the deviation parameter found from the fitted Gaussian function is in squared form.

In our analysis we have used R-R interval time series of $`22`$$`24`$ hours, corresponding to $`80.000`$$`100.000`$ data points. Our data is actually interval data, i.e., it consists of a sequence of R-R interval values. It is then convenient to count the delay in our analysis in terms of heart-beats rather than seconds, i.e., we have not used cumulative time as time variable but the beat index. However, since the R-R interval values vary a lot within the used delay range, the beat index actually gives a delay as if computed with the average beat-rate. We have tested both methods and found only minor differences between them (in the details of the functions $`g`$ and $`h`$). We will show later that the functional forms of $`g`$ and $`h`$ are quite insensitive on the time delay, and since this holds for both methods we will use the more convenient beat index.

## III Results

### III.1 A typical case

In Fig. 3 we have presented results obtained for a particular case using the method described earlier. The value of the delay parameter $`\tau `$ was $`500`$ beats, and the number of boxes used to construct local distributions was $`150`$. Distributions were fitted using Gaussian function. The $`g(X`$) function, the deterministic part of the system, is displayed on the left panel in Fig 3. It has a very clear and simple functional form (between the vertical lines) which is typical for systems exhibiting bistable behaviour Bergé et al. (1986); van Kampen (1981). The function crosses the zero line three times, these crossings are the fixed points of the system. The fixed points marked with arrows are stable: without any noise term these points attract all nearby states, because the control function $`g(X)`$ is locally decreasing. The middle fixed point is repulsive. Due to the stochastic part the system has a tendency to jump between the stable points if the amplitude of the noise is high enough. Far away from the stable points $`g(X)`$ increases or decreases strongly forcing the system rapidly back to oscillate around the stable points. The amplitude of the stochastic part of the system, function $`h(X)`$, is almost constant except between the stable points where it has a clear maximum (the middle panel in Fig. 3). One interpretation is that the system has a larger inherent freedom to oscillate randomly when the trajectory is between the stable points but outside this range the character of the system is more deterministic. From the physiological point of view this kind of dynamics can be useful since it lets the R-R interval to wander most of the time but prevents it from escaping too far away from the normal range. On the right panel in Fig. 3 we have shown the correlation coefficient of each local distribution. Most of the time the correlation is remarkably high, about $`0.85`$$`0.95`$, but near the largest and smallest $`X`$ values there are only rather few data points and therefore the corresponding distributions do not have clear Gaussian shape resulting with lower correlation. The high average correlation value is a clear indication that the noise in this system is really Gaussian type. We have used the value of $`0.8`$ as a threshold level, and the corresponding range is marked with the vertical lines in Fig. 3.

What is remarkable in this description is that the functional forms of $`g(X)`$ and $`h(X)`$ are fairly independent of the delay parameter $`\tau `$ in a rather extensive delay range, typically $`100`$$`1000`$ beats (corresponding to $`2`$$`20`$ minutes). In Fig. 4 we have plotted the functions $`g(X)`$ and $`h(X)`$ for a range of $`\tau `$ values. The $`g`$-function is practically $`\tau `$ independent, except for the shortest R-R intervals, where some cumulative effects show up. The $`h`$-function seems to grow very slowly as $`\tau `$ increases. For still smaller delay values $`g(X)`$ is more flat and $`h(X)`$ is more scattered, and for longer delays $`g(X)`$ is typically a straight line and $`h(X)`$ is constant. Behavior at these extremes can be easily understood by recalling that when the time scale is small, the heart rate system is clearly multidimensional depending directly on blood pressure, respiration and other rapidly changing physiological variables and our 1-dimensional description is no longer valid. On the other hand, if the delay parameter is very large, we cannot reconstruct the local dynamics in terms of local distributions, we just get the global distribution that is independent of dynamics and no longer Gaussian. In the right panel of Fig. 4 we have given the values of the $`g(X)`$ function at $`X=780`$ ms (marked with a vertical dashed line in the left panel) computed with delays of $`5`$$`10240`$ beats. We can see a plateau in the delay range of $`100`$$`1000`$ beats which means that the $`g(X)`$ curves for these delays are bundled. In principle the curves for a delay of $`2\tau `$ should be obtainable by iterating (2) with delay $`\tau `$. Direct numerical calculations of joint probabilities using experimentally determined $`g(X)`$ (within $`100`$$`1000`$ beats delay range) indicate that $`g(X)`$ and $`h(X)`$ do not change significantly in one iteration, mostly because in our case the Gaussian distribution is not so narrow. In general iterations tend to sharpen the bends in $`g(X)`$ and this feature is indeed visible in Fig. 4. The small $`\tau `$-dependency of $`g`$ and $`h`$ in the range of short R-R intervals can then be interpreted either as the expected result from repeated iterations, or as a sign of higher order dynamics: possibly the heart rate regulation system is more complex when the system must readjust at a fast heart rate.

### III.2 Variation between subjects

In order to find whether different subjects have are any common features in the deterministic and stochastic parts $`g(X)`$ and $`h(X)`$ we analyzed the data from $`27`$ healthy subjects of various age and gender \[$`18`$ cases from PhysioBank Goldberger et al. (2000) and $`9`$ cases from Kuopio University Hospital\]. Analyses were done using the same parameter values as in Fig. 3. The deterministic part, the $`g(X)`$ function, is displayed in Fig. 5 for a set of $`9`$ typical cases. The most common form for this function is the bistable type, already shown in Fig. 3, where the $`g(X)`$ function has three zeroes, and $`60`$% of all cases can be classified to this group (cases 1–5 in Fig. 5). The next most common group, $`25`$% of all cases, has a $`g(X)`$ function with $`5`$ zeroes, a kind of double pitchfork system (cases 6 and 7 in Fig. 5). We also found $`3`$ cases where the $`g(X)`$ function seems to have even more zeroes (case 8 in Fig. 5). Only very few cases could not be clearly classified as bi- or multistable. In these cases it can be difficult to interpret the results. It is possible that the dynamical variable did not explore the whole state phase, and therefore we can see only part of the $`g(X)`$ function; for example case 9 in Fig. 5, where the system has only one stable fixed point and no unstable points at all, can be an example of this. The stochastic parts (function $`h(X)`$) are fairly similar: they are almost constant except that in all cases there are maxima on the R-R interval ranges between the stable fixed points of the deterministic part, as in the example in Fig. 3.

The description given by equation (2) contains both a deterministic and a stochastic component. It is an important to realize that the stochastic part is not a small perturbation but in fact forms an essential part of the description, furthermore it is $`10`$$`20`$ times higher than the measurement noise (uncertainty in detecting the position of the R-peak), which is typically only $`2`$$`5`$ ms. One way to compare the deterministic and stochastic components is to note that the size of the bend in the $`g(X)`$ function is of the order of $`30`$ to $`50`$ ms, while the average size of the $`h(X)`$ function is about $`70`$ to $`110`$ ms, as can be seen in Fig. 3. \[The extraction of small details in the $`g(X)`$ function under such noise is of course possible only because the noise is so cleanly Gaussian.\] On the other hand, the distance between the stable fixed points in the $`g(X)`$ function is of the order of $`50`$ to $`250`$ ms, and therefore the probability that the systems jumps between stable points is not extremely high, but nevertheless possible. It is also possible that external factors drive the system from one stable point to another, since during night-time the mean R-R interval is typically longer than during day-time \[although the R-R interval can abruptly jump to the faster rate also during the night, as can be seen on the lower panel in Fig. 1\].

### III.3 Same subject at different times

If the model (2) were to describe true heart rate dynamics, the functions $`g(X)`$ and $`h(X)`$ should have some constant features specific for each subject. In order to look at this aspect we made two recordings from the same subject within 4 days, the results are shown in Fig. 6. In general the deterministic and stochastic parts from different recordings are remarkably similar both having clear bistable character. In the R-R interval range of $`500`$$`800`$ ms the results are almost identical and the only difference seems to be a scaling towards the shorter R-R intervals in the $`800`$$`1100`$ ms range of the second recording. In the first recording the mean value of R-R interval calculated over the 24 hour period was $`781`$ ms and in the second one $`726`$ ms. Therefore in the second recording the shortest R-R intervals are significantly more frequent and this can affect on the analysis results. These deviations could also reflect true changes on the underlying control system: it is well know that there are daily variations on functions of the autonomic nervous system.

### III.4 Surrogate analysis

As a further validity check we also performed surrogate analysis Schreiber and Schmitz (2000); Theiler et al. (1992) in order to eliminate the possibility that the results are generated just from peculiar distribution of the R-R intervals imitating real dynamics. For this purpose the data was shuffled by dividing it into sections of equal size which were then repositioned randomly. As a result we get a new time series where the dynamical structure has been partially destroyed depending on the section size. Results of this surrogate analysis are shown in Fig. 7. The top panels display the deterministic $`g(X)`$ and stochastic $`h(X)`$ parts of the system and the correlation coefficient without any data shuffling (row A in Fig. 7). On the next row (row B in Fig. 7) we have used sections of $`800`$ data points for shuffling. There are only small changes in the deterministic part, but the correlation has decreased noticeably. When the section size is $`400`$ (row C in Fig. 7) we can no longer see the bistable character in the deterministic part, the stochastic part is flat with higher mean level, and the average level of the correlation coefficient has dropped well below our threshold value $`0.8`$. With still smaller section sizes the results do not change any further. In this analysis we have used the same delay of $`500`$ data points as previously and when the section size used in the shuffling process is less than this delay all dynamical properties disappear, as expected in the case of true time evolution. Therefore we conclude that our results are derived from the dynamical properties of the heart beat data, and not from their overall statistical characteristics.

## IV Conclusion

Our results indicate that the human heart-rate control dynamics can be accurately modeled with the 1-dimensional stochastic difference equation (2), where the time delay parameter is within $`2`$$`20`$ minutes. Stochasticity is an integral part of the dynamics, and in this delay range the effects of other variables are either embedded into the stochastic part of the system or averaged over time with no net effect. It is remarkable that the form of the control function $`g(X)`$ is similar from case to case. Their typically bistable character is also well justified on common physiological grounds. From this initial study we cannot yet identify what kind of dynamical structure is typical for healthy subjects (although our results already indicate that a simple bistable system is most common feature) and therefore the model cannot yet be used directly for clinical work, for that purpose one needs extensive demographic studies. We can nevertheless speculate that the form of the control function $`g(X)`$ should tell us something about the health of the subject. Also, some of the current knowledge based on statistical measures of heart rate time series can probably be explained within the framework of our model. Another interesting observation is the importance of the stochastic part, it could be the result of integrating the effects of a more detailed control mechanism over time, but it could also reflect some truly stochastic internal and external influences.

###### Acknowledgements.

We thank T. Laitinen from Kuopio University Hospital, Department of Clinical Physiology, for providing nine electrocardiogram recordings. This work was partially supported by the Academy of Finland.


text:
# Quark distributions in polarized 𝜌 meson and its comparison with those in pion

## 1 Introduction

Structure functions are one of the most significant characteristics of the inner structure of hadrons . For nucleon and pion there are experimental results (see, e.g. - for the nucleon, and -, for the pion). For other hadrons, one should use some models based on additional suppositions about inner structure. A very significant question about quark distribution dependence on hadron polarization has no model independent answer (moreover in various models it is usually simply supposed that there is no significant influence).

That is why determination of quark distribution functions in a model-independent way in QCD sum rules based only on QCD and the operator product expansion (OPE) seems to be very important task especially for polarized hadron (in this talk we will discuss polarized $`\rho `$ meson case).

QCD sum rules for valence quark distribution at intermediate $`x`$ was suggested in  and developed in -. The method was based on the fact, that the imaginary part (in $`s`$\- channel) of 4-point correlator, corresponding to the forward scattering of a virtual photon on the current with the quantum numbers of hadron of interest is dominated by contribution of small distances (at intermediate $`x`$) if the virtualities of the photon and hadronic current ($`q^2`$ and $`p^2`$, respectively) are large and negative $`|q^2||p^2|R_c^2`$ , where $`R_c`$ is the confinement radius. So the operator product expansion in this $`x`$ region is applicable. Then, comparing dispersion representation of the forward amplitude in terms of physical states with that, calculated in OPE and using Borel transformation to suppress higher resonance contributions, one can find quark distribution functions at intermediate $`x`$. Unfortunately, the accuracy of the sum rules for nucleon, obtained in , are bad, (especially for $`d`$-quark), moreover it was found to be impossible to calculate quark distributions in the $`\pi `$ and $`\rho `$ mesons in this way. The reason is that the sum rules in the form used in  have a serious drawback.

This drawback comes from the fact that contribution of nondiagonal transitions wasn’t suppressed in sum rules by borelization, and special additional procedure, which was used in  to supress them is incorrect for such sum rules. There is no time to discuss this in details, the full analisys can be found in our papers  and . In this papers the modified method of calculating of the hadron structure functions (quark distributions in hadrons) was suggested. The problem of supressing the nondiagonal terms is eliminated and valence quark distributions in pion and polarizied $`\rho `$ meson were calculated.

## 2 The method

Let me briefly present the method. Consider the non-forward 4-point correlator

$$\mathrm{\Pi }(p_1,p_2;q,q^{})=id^4xd^4yd^4ze^{ip_1x+iqyip_2z}$$

$$\times 0|T\{j^h(x),j^{el}(y),j^{el}(0),j^h(z)\}|0$$
(1)
Here $`p_1`$ and $`p_2`$ are the initial and final momenta carried by hadronic the current $`j^h`$, $`q`$ and $`q^{}=q+p_1p_2`$ are the initial and final momenta carried by virtual photons and Lorentz indices are omitted. It will be very essential for us to consider unequal $`p_1`$ and $`p_2`$ and treat $`p_1^2`$ and $`p_2^2`$ as two independent variables. However, we may put $`q^2=q^2`$ and $`t=(p_1p_2)^2=0`$. The general form of the double dispersion relation (in $`p_1^2,p_2^2`$) of the imaginary part of $`\mathrm{\Pi }(p_1^2,p_2^2,q^2,s)`$ in the $`s`$ channel has the form

$$Im\mathrm{\Pi }(p_1^2,p_2^2,q^2,s)=a(q^2,s)+\underset{0}{\overset{\mathrm{}}{}}\frac{\phi (q^2,s,u)}{up_1^2}𝑑u+\underset{0}{\overset{\mathrm{}}{}}\frac{\phi (q^2,s,u)}{up_2^2}𝑑u$$

$$+\underset{0}{\overset{\mathrm{}}{}}𝑑u_1\underset{0}{\overset{\mathrm{}}{}}𝑑u_2\frac{\rho (q^2,s,u_1,u_2)}{(u_1p_1^2)(u_2p_2^2)}$$
(2)
Double Borel transformation in $`p_1^2`$ and $`p_2^2`$ to (2) eliminates three first terms and we have

$$_{M_1^2}_{M_2^2}Im\mathrm{\Pi }(p_1^2,p_2^2,q^2,s)=\underset{0}{\overset{\mathrm{}}{}}𝑑u_1\underset{0}{\overset{\mathrm{}}{}}𝑑u_2\rho (q^2,s,u_1,u_2)exp\left[\frac{u_1}{M_1^2}\frac{u_2}{M_2^2}\right]$$
(3)
where $`M_1^2`$ and $`M_2^2`$ are the squared Borel mass. The integration region with respect to $`u_1,u_2`$ may be divided into four areas:

I. $`u_1<s_0,u_2<s_0;`$

II.$`u_1<s_0,u_2>s_0;`$

III. $`u_1>s_0,u_2<s_0;`$

IV. $`u_1,u_2>s_0.`$

Here $`s_0`$ is the continuum threshold in the standard QCD sum rule model of the hadronic spectrum with one lowest resonance plus continuum. Area I obviously corresponds to the resonance contribution and the spectral density in this area can be written as

$$\rho (u_1,u_2,x,Q^2)=g_h^22\pi F_2(x,Q^2)\delta (u_1m_h^2)\delta (u_2m_h^2),$$
(4)
where $`g_h`$ is defined as

$$0|j_h|h=g_h$$
(5)
(For simplicity we consider the case of the Lorentz scalar hadronic current. The necessary modifications for the $`\pi `$ and $`\rho `$ mesons will be presented below). If the structure, proportional to $`P_\mu P_\nu `$ $`[P_\mu =(p_1+p_2)_\mu /2]`$ is considered in $`Im\mathrm{\Pi }(p_1,p_2,q,q^{})`$, then in the lowest twist approximation $`F_2(x,Q^2)`$ is the structure function of interest.

In area IV, where both variables $`u_{1,2}`$ are far from the resonance region, the nonperturbative effects may be neglected and, as usual in the sum rules, the spectral function of a hadron state is described by the bare loop spectral function $`\rho ^0`$ in the same region:

$$\rho (u_1,u_2,x)=\rho ^0(u_1,u_2,x)$$
(6)
In areas II and III one of the variables is far from the resonance region, but another is in the resonance region, and the spectral function in this region is some unknown function $`\rho =\psi (u_1,u_2,x)`$, which corresponds to the transitions like $`hcontinuum`$.

Taking all these facts into account, the physical side of the sum rule (3) can be rewritten as

$$\widehat{B}_1\widehat{B}_2[Im\mathrm{\Pi }]=2\pi F_2(x,Q^2)g_h^2e^{m_h^2(\frac{1}{M_1^2}+\frac{1}{M_2^2})}+\underset{0}{\overset{s_0}{}}𝑑u_1\underset{s_0}{\overset{\mathrm{}}{}}𝑑u_2\psi (u_1,u_2,x)e^{(\frac{u_1}{M_1^2}+\frac{u_2}{M_2^2})}$$

$$+\underset{s_0}{\overset{\mathrm{}}{}}𝑑u_1\underset{0}{\overset{s_0}{}}𝑑u_2\psi (u_1,u_2,x)e^{(\frac{u_1}{M_1^2}+\frac{u_2}{M_2^2})}+\underset{s_0}{\overset{\mathrm{}}{}}\underset{s_0}{\overset{\mathrm{}}{}}𝑑u_1𝑑u_2\rho ^0(u_1,u_2,x)e^{(\frac{u_1}{M_1^2}+\frac{u_2}{M_2^2})}$$
(7)
In what follows, we put $`M_1^2=M_2^22M^2`$. (As was shown in , the values of the Borel parameters $`M_1^2,M_2^2`$ in the double Borel transformation are about twice as large as those in the ordinary ones).

One of the advantages of this method is that after the double Borel transformation the unknown contributions of II and III areas \[the second and the third terms in (7)\] are exponentially suppressed. So, and we would like to emphasise this, we do not need in any additional artificial procedure such as the differentiation with respect to the Borel mass. Using standard duality arguments, we estimate the contribution of all nonresonance region (i.e., areas II, III, IV) as the contribution of the bare loop in the same region and demand their value to be small (less than 30%). Finally, equating the physical and QCD representations of $`\mathrm{\Pi }`$ one can write the following sum rules:

$$Im\mathrm{\Pi }_{QCD}^0+\text{Power  correction}=2\pi F_2(x,Q^2)g_h^2e^{m_h^2(\frac{1}{M_1^2}+\frac{1}{M_2^2})}$$

$$Im\mathrm{\Pi }_{QCD}^0=\underset{0}{\overset{s_0}{}}\underset{0}{\overset{s_0}{}}\rho ^0(u_1,u_2,x)e^{\frac{u_1+u_2}{2M^2}}$$
(8)
It is worth mentioning that if we consider the forward scattering amplitude from the very beginning, (i.e. put $`p_1=p_2=p`$ as in -) and perform the Borel transformation in $`p^2`$, then the contributions of the second and third terms in (2),in contrast to (7), would not be suppressed comparing with the lowest resonance contribution we are interesting in. They correspond to the nondiagonal transition matrix elements discussed in the Introduction.

## 3  Quark distributions in pion

Let me briefly discuss the main points of the calculation and the results for the pion structure function, which can be treated as a check of the accuracy of the method due to fact that for pion the experimental results are avaible. To find the pion structure function by the method, described in the previous section, one should consider the imaginary part of 4-point correlator (1) with two axial and two electromagnetic currents. Since $`\overline{d}(x)=u(x)`$, it is enough to find the distribution of the valence $`u`$ quark in $`\pi ^+`$. The most suitable choice of the axial current is

$$j_{\mu 5}=\overline{u}\gamma _\mu \gamma _5d$$
(9)
The electromagnetic current is chosen as $`u`$-quark current with the unit charge

$$j_\mu ^{el}=\overline{u}\gamma _\mu u$$
(10)
The most convenient tensor structure, which is chosen to construct the sum rule, is a structure proportional to $`P_\mu P_\nu P_\lambda P_\sigma /\nu `$.

The sum rule for the valence $`u`$-quark distribution in the pion in the bare loop approximation is :

$$u_\pi (x)=\frac{3}{2\pi ^2}\frac{M^2}{f_\pi ^2}x(1x)(1e^{s_0/M^2})e^{m_\pi ^2/M^2},$$
(11)
where $`s_0`$ is the continuum threshold. In  the following corrections to Eq. (11) were taken into account:

1. Leading order (LO) perturbative corrections proportional to $`ln(Q^2/\mu ^2)`$ , where $`\mu ^2`$ is the normalization point. In what follows, the normalization point will be chosen to be equal to the Borel parameter $`\mu ^2=M^2`$.

2. Power corrections – higher order terms of OPE. Among them, the dimension-4 correction, proportional to the gluon condensate $`0|\frac{\alpha _s}{\pi }G_{\mu \nu }^nG_{\mu \nu }^n|0`$ was first taken into account, but it was found that the gluon condensate contribution to the sum rule vanishes after double borelization. There are two types of vacuum expectation values of dimension 6. One, involving only gluonic fields:

$$\frac{g_s}{\pi }\alpha _sf^{abc}0|G_{\mu \nu }^aG_{\nu \lambda }^bG_{\lambda \mu }^c|0$$
(12)
and the other, proportional to the four-quark operators

$$0|\overline{\psi }\mathrm{\Gamma }\psi \overline{\psi }\mathrm{\Gamma }\psi |0$$
(13)
It was shown in  that terms of the first type cancel in the sum rule and only terms of the second type survive. For the latter, one may use the factorization hypothesis which reduces all the terms of this type to the square of the quark condensate.

A remark is in order here. As was mentioneed in the Introduction, the present approach is invalid at small and large $`x`$. No-loop 4-quark condensate contributions are proportional to $`\delta (1x)`$ and, being outside of the applicability domain of the approach, cannot be taken into account. In the same way, the diagrams, which can be considered as radiative corrections to those, proportional to $`\delta (1x)`$ must be also omitted.

All dimension-6 power corrections to the sum rule were calculated in  and the final result is given by (the pion mass is neglected):

$$xu_\pi (x)=\frac{3}{2\pi ^2}\frac{M^2}{f_\pi ^2}x^2(1x)[(1+\left(\frac{a_s(M^2)ln(Q_0^2/M^2)}{3\pi }\right)$$

$$\times (\frac{1+4xln(1x)}{x}\frac{2(12x)lnx}{1x}))(1e^{s_0/M^2})$$

$$\frac{4\pi \alpha _s(M^2)4\pi \alpha _sa^2}{(2\pi )^43^72^6M^6}\frac{\omega (x)}{x^3(1x)^3}],$$
(14)
where $`\omega (x)`$ is the fourth degree polynomial in $`x`$,

$$a=(2\pi )^20|\overline{\psi }\psi |0$$
(15)
$$\omega (x)=5784x^41140x^320196x^2$$

$$+20628x8292)ln(2)+4740x^4+8847x^3$$

$$+2066x^22553x+1416$$
(16)
The function $`u_\pi (x)`$ may be used as an initial condition at $`Q^2=Q_0^2`$ for solution of QCD evolution equations (Dokshitzer-Gribov-Lipatov-Altarelli-Parisi equations).

In the numerical calculations we choose: the effective $`\mathrm{\Lambda }_{QCD}^{LO}=200MeV,Q_0^2=2GeV^2,\alpha _sa^2(1GeV^2)=0.13GeV^6`$ . The continuum threshold was varied in the interval $`0.8<s_0<1.2GeV^2`$ and it was found, that the results depend only slightly on it’s variation. The analysis of the sum rule (14) shows, that it is fulfilled in the region $`0.15<x<0.7`$; the power corrections are less than 30% and the continuum contribution is small ($`<25`$%). The stability in the Borel mass parameter $`M^2`$ dependence in the region $`0.4GeV^2<M^2<0.6GeV^2`$ is good. The result of our calculation of the valence quark distribution $`xu_\pi (x,Q_0^2)`$ in the pion is shown in Fig. 1. In Fig. 1, we plot also the valence $`u`$-quark distribution found in  by fitting the data on the production of $`\mu ^+\mu ^{}`$ and $`e^+e^{}`$ pairs in pion-nucleon collisions (Drell-Yan process). When comparing with the distribution found here it should be remembered, that the accuracy of our curve is of order of $`1020\%`$, (see discussion in Section 5). The $`u`$-quark distribution found from the experiment is also not free from uncertainties (at least 10-20%, see  \- ). For all these reasons, we consider the agreement of two curves as being good.

Assume, that $`u_\pi (x)1/\sqrt{x}`$ at small $`x\begin{array}{c}<\hfill \\ \hfill \end{array}0.15`$ according to the Regge behaviour and $`u_\pi (x)(1x)^2`$ at large $`x\begin{array}{c}>\hfill \\ \hfill \end{array}0.7`$ according to quark counting rules. Then, matching these functions with (14), one may find the numerical values of the first and the second moments of the $`u`$-quark distribution:

$$_1=\underset{0}{\overset{1}{}}u_\pi (x)𝑑x0.84(0.85)$$
(17)
$$_2=\underset{0}{\overset{1}{}}xu_\pi (x)𝑑x0.21(0.23)$$
(18)
where the values in the parentheses correspond to behavior $`u_\pi (x)(1x)`$ at large $`x`$. The results depend only slightly on the matching points (not more than 5%, when the lower matching point is varied in the region 0.15 - 0.2 and the upper one in the region 0.65 - 0.75). The moment $`_1`$ has the meaning of the number of $`u`$ quarks in $`\pi ^+`$ and it should be $`_1=1`$. The deviation of (17) from 1 characterizes the accuracy of our calculation. The moment $`_2`$ has the meaning of the pion momentum fraction carried by the valence $`u`$ quark. Therefore, the valence $`u`$ and $`\overline{d}`$ quarks carry about 40% of the total momentum.

The valence $`u`$-quark distribution in the pion was calculated recently in the instanton model . At intermediate $`x`$, the values of $`xu_\pi (x)`$ found in  are not more than 20% higher that our results. Recently, the pions valence quark momentum distribution using a model, based on the Dyson-Schwinger equation was also calculated . Our results are in a reasonable agreement with the results of . Our estimation (18) of the second moment of the valence quark distribution can be also compared with the calculation of the second moment of the total (valence plus sea) quark distribution in the pion . The value, obtained in , is 0.6 for the total momentum carried by all quarks (valence plus sea) with the accuracy about 10%. Taking into account, that sea quarks usually supposed to carry 15% of the total momentum, one can estimate from the result of , that the second moment of one valence quark distribution should be about 20-22%, which is in a good agreement with our result (18). The quark distribution in the pion was calculated also in  by using sum rules with nonlocal condensates. Unfortunately it is impossible to perform the comparison directly, because the quark distribution is calculated in  only at very low normalization point. But, comparig his result with different models, the author of this paper arrived at the conclusion, that the result is in agreement with experimental data at $`Q^2`$=20 $`GeV^2`$ within the accuracy about 20%. Our result is close also to experimental fit, so we can believe that our results are in an agreement with those of .

## 4 Quark distributions in $`\rho `$ meson

Let us calculate valence $`u`$-quark distribution in the $`\rho ^+`$ meson. The choice of hadronic current is evident

$$j_\mu ^\rho =\overline{u}\gamma _\mu d$$
(19)
The matrix element $`\rho ^+j_\mu ^\rho 0`$ is given by

$$\rho ^+j_\mu ^\rho 0=\frac{m_\rho ^2}{g_\rho }e_\mu $$
(20)
where $`m_\rho `$ is the $`\rho `$-meson mass, $`g_\rho `$ is the $`\rho \gamma `$ coupling constant, $`g_\rho ^2/4\pi =1.27`$ and $`e_\mu `$ is the $`\rho `$ meson polarization vector. Consider the coordinate system, where the collision of the $`\rho `$-meson with momentum $`p`$ and the virtual photon with momentum $`q`$ proceeds along $`z`$-axes. Averaging over $`\rho `$ polarizations is given by the formulas:

$$e_\mu ^Le_\nu ^L=\left(q_\mu \frac{\nu p_\mu }{m_\rho ^2}\right)\left(q_\nu \frac{\nu p_\nu }{m_\rho ^2}\right)\frac{m_\rho ^2}{\nu ^2q^2m_\rho ^2}$$
(21)
for longitudinally polarized $`\rho `$ and

$$\underset{T,r}{}e_\mu ^re_\nu ^r=\left(\delta _{\mu \nu }\frac{p_\mu p_\nu }{m_\rho ^2}\right)\frac{m_\rho ^2}{\nu ^2q^2m_\rho ^2}\left(q_\mu \frac{\nu p_\mu }{m_\rho ^2}\right)\left(q_\nu \frac{\nu p_\nu }{m_\rho ^2}\right)$$
(22)
for transversally polarized $`\rho `$.

The imaginary part of the forward $`\rho \gamma `$ scattering amplitude $`W_{\mu \nu \lambda \sigma }`$ (before multiplication by $`\rho `$ polarizations) satisfies the equations $`W_{\mu \nu \lambda \sigma }q_\mu =`$ $`W_{\mu \nu \lambda \sigma }q_\nu =`$ $`W_{\mu \nu \lambda \sigma }p_\lambda =`$ $`W_{\mu \nu \lambda \sigma }p_\sigma =0`$, which follow from current conservation. The indices $`\mu ,\nu `$ refer to the initial and final photon; $`\lambda ,\sigma `$ – to the initial and final $`\rho `$. The general form of $`W_{\mu \nu \lambda \sigma }`$ is:

$$W_{\mu \nu \lambda \sigma }=[(\delta _{\mu \nu }\frac{q_\mu q_\nu }{q^2})(\delta _{\lambda \sigma }\frac{p_\lambda p_\sigma }{m_\rho ^2})A(\delta _{\mu \nu }\frac{q_\mu q_\nu }{q^2})(q_\lambda \frac{\nu p_\lambda }{m_\rho ^2})(q_\sigma \frac{\nu p_\sigma }{m_\rho ^2})B$$

$$\left(p_\mu \frac{\nu q_\mu }{q^2}\right)\left(p_\nu \frac{\nu q_\nu }{q^2}\right)\left(\delta _{\lambda \sigma }\frac{p_\lambda p_\sigma }{m_\rho ^2}\right)C+\left(p_\mu \frac{\nu q_\mu }{q^2}\right)\left(p_\nu \frac{\nu q_\nu }{q^2}\right)$$

$$\times (q_\lambda \frac{\nu p_\lambda }{m_\rho ^2})(q_\sigma \frac{\nu p_\sigma }{m_\rho ^2})D]$$
(23)
where A,B,C and D are the invariant functions. By averaging Eq. (23) over polarizations for longitudinal and transverse $`\rho `$ mesons, one can find that the structure function $`F_2(x)`$ proportional to $`p_\mu p_\nu `$ is given in the scaling limit $`(\nu ^2q^2m_\rho ^2)`$ by the contribution of invariants $`C+(\nu ^2/m_\rho ^2)D`$ and $`C`$ in the cases of the longitudinal and the transversal $`\rho `$ meson, respectively. This means that, in the forward scattering amplitude $`W_{\mu \nu \lambda \sigma }`$ (23), one must separate the structure proportional to $`p_\mu p_\nu p_\lambda p_\sigma `$ in the first case and the structure $`p_\mu p_\nu \delta _{\lambda \sigma }`$ in the second case.

Consider now the non-forward 4-point correlator

$$\mathrm{\Pi }_{\mu \nu \lambda \sigma }(p_1,p_2;q,q^{})=id^4xd^4yd^4ze^{ip_1x+iqyip_2z}$$

$$\times 0T\{j_\lambda ^\rho (x),j_\mu ^{el}(y),j_\nu ^{el}(0),j_\sigma ^\rho (z)\}0,$$
(24)
where the currents $`j_\mu ^{el}(x)`$ and $`j_\lambda ^\rho (x)`$ are given by Eq.(10) and (20). It is evident from the consideration, that in the non-forward amplitude the most suitable tensor structure for determination of $`u`$-quark distribution in the longitudinal $`\rho `$ meson is that proportional to $`P_\mu P_\nu P_\sigma P_\lambda `$, while $`u`$-quark distribution in the transverse $`\rho `$ can be found by considering the invariant function at the structure ($`P_\mu P_\nu \delta _{\lambda \sigma }`$).

In the case of longitudinal $`\rho `$ meson the tensor structure, that is separated is the same as in the case of the pion. Since bare loop contributions for vector and axial hadronic currents are equal at $`m_q=0`$, the only difference from the pion case is the normalization. It can be shown, that $`u`$-quark distribution in the longitudinal $`\rho `$ meson can be found from Eq.(14) by substituting $`m_\pi m_\rho `$, $`f_\pi m_\rho /g_\rho `$ and, therefore, one can easily write down sum rules for this distribution:

$$xu_\rho ^L(x)=\frac{3}{2\pi ^2}M^2\frac{g_\rho ^2}{m_\rho ^2}e^{m_\rho ^2/M^2}x^2(1x)[(1+\left(\frac{a_s(M^2)ln(Q_0^2/M^2)}{3\pi }\right)$$

$$\times (\frac{1+4xln(1x)}{x}\frac{2(12x)lnx}{1x}))(1e^{s_0/M^2})\frac{\alpha _s(M^2)\alpha _sa^2}{\pi ^23^72^6M^6}\frac{\omega (x)}{x^3(1x)^3}],$$
(25)
where $`a`$ and $`\omega (x)`$ are given by Eqs.(15) and (16), respectively. Sum rules for $`u_\rho ^L(x)`$ are satisfied in the wide $`x`$ region: $`0.1<x<0.85`$. The Borel mass $`M^2`$ dependence is weak in the whole range of $`x`$, except $`x0.15`$ and $`x0.7`$. Figure. 2 presents $`xu_\rho ^L(x)`$ as a function of $`x`$. The values $`M^2=1`$ GeV<sup>2</sup> and $`s_0=1.5`$ GeV<sup>2</sup>, $`Q_0^2=4`$ GeV<sup>2</sup> were chosen, the parameters $`\mathrm{\Lambda }_{QCD}^{LO}`$ and $`\alpha _sa^2`$ are the same as in the calculation of $`xu_\pi (x)`$.

Let us now consider the case of transversally polarizied $`\rho `$-meson, i.e., the term proportional to the structure $`P_\mu P_\nu \delta _{\lambda \sigma }`$. The bare loop contribution (with leading order perturbative corrections) is

$$u^T(x)=\frac{3}{8\pi ^2}\frac{g_\rho ^2}{m_\rho ^4}e^{m_\rho ^2/M^2}M^4E_1\left(\frac{s_0}{M^2}\right)\phi _0(x)$$

$$\left[1+\frac{ln(Q_0^2/\mu ^2)\alpha _s(\mu ^2)}{3\pi }\left((4x1)/\phi _0(x)+4ln(1x)\frac{2(12x+4x^2)lnx}{\phi _0(x)}\right)\right]$$
(26)
where

$$\phi _0(x)=12x(1x)$$
(27)
We consider now the power correction contribution to the sum rules. The power correction of the lowest dimension is proportional to the gluon condensate $`G_{\mu \nu }^qG_{\mu \nu }^q`$ with $`d=4`$. The $`G_{\mu \nu }^qG_{\mu \nu }^q`$ correction was calculated in the standard way in the Fock-Schwinger gauge $`x_\mu A_\mu =0`$ .

The quark propagator $`iS(x,y)=\psi (x)\psi (y)`$ in the external field $`A_\mu `$ has the well-known form . In contrast to the pion case, the $`G_{\mu \nu }^aG_{\mu \nu }^a`$ correction for transversally polarized $`\rho (\rho _T)`$, does not vanish.

$$Im\mathrm{\Pi }_T^{(d=4)}=\frac{\pi }{8x}0\frac{\alpha _s}{\pi }G_{\mu \nu }^20$$
(28)
There are a great number of loop diagrams for $`d=6`$ correction. It is convenient to divide them into two types and discuss these types separately. Type I diagrams are those, in which only the interaction with the external gluon field is taken into account and type II diagrams are those in which expansion of the quark field is also taken into account.

Discuss briefly the special features of calculating diagrams of these two types. The type-I diagrams are obviously proportional to $`0g^3f^{abc}G_{\mu \nu }^aG_{\alpha \beta }^bG_{\rho \sigma }^c0`$, $`0D_\rho G_{\mu \nu }^aD_\tau G_{\alpha \beta }^a0`$ and $`0G_{\mu \nu }^aD_\rho D_\tau G_{\alpha \beta }^a0`$.

One may demonstrate  that these tensor structures are proportional to two vacuum averages

$$0g^2j_\mu ^20\text{and}0g^3G_{\mu \nu }^aG_{\nu \rho }^bG_{\rho \mu }^cf^{abc}0.$$

By using the factorization hypothesis the first of them, $`0g^2j_\mu ^20`$ easily reduces to $`g\overline{\psi }\psi ^2`$, which is well known:

$$0g^2j_\mu ^20=(4/3)[0g\overline{\psi }\psi 0]^2.$$
(29)
But $`0g^3G_{\mu \nu }^aG_{\nu \rho }^bG_{\rho \mu }^cf^{abc}0`$ is not well known; there are only some estimates based on the instanton model . In contrast to $`\pi `$ and $`\rho _L`$-meson cases, the terms proportional to
$`0g^3f^{abc}G_{\mu \nu }^aG_{\nu \rho }^b`$ $`G_{\rho \mu }^c0`$ are not cancelled for $`\rho _T`$ and one should estimate it. The estimation based on the instanton model  gives

$$g^3f^{abc}G_{\mu \nu }^aG_{\nu \rho }^bG_{\rho \mu }^c=\frac{48\pi ^2}{5}\frac{1}{\rho _c^2}0(\alpha _s/\pi )G_{\mu \nu }^20,$$
(30)
where $`\rho _c`$ is the effective instanton radius.

Among type-II diagrams only those, in which the interaction with the vacuum takes place inside the loop are considered. Such diagrams cannot be treated as the evolution of any non-loop diagrams and are pure power corrections of dimension 6.

The total number of $`d=6`$ diagrams is enormous – about 500. Collecting the results we get finally the following sum rules for the valence $`u`$-quark distribution in the transversally polarized $`\rho `$ meson.

$$xu_\rho ^T(x)=\frac{3}{8\pi ^2}g_\rho ^2e^{m_\rho ^2/M^2}\frac{M^4}{m_\rho ^4}x\{\phi _0(x)E_1\left(\frac{s_0}{M^2}\right)[1+\frac{1}{3\pi }ln\left(\frac{Q_0^2}{M^2}\right)\alpha _s(M^2)(\frac{(4x1)}{\phi _0(x)}+$$

$$+4ln(1x)\frac{2(12x+4x^2)lnx}{\phi _0(x)})]\frac{\pi ^2}{6}\frac{0|(\alpha _s/\pi )G^2|0}{M^4x^2}$$

$$+\frac{1}{2^83^5M^6x^3(1x)^3}0|g^3f^{abc}G_{\mu \nu }^aG_{\nu \lambda }^bG_{\lambda \mu }^c|0\xi (x)$$

$$+\frac{\alpha _s(M^2)(\alpha _sa^2)}{2^53^8\pi ^2M^6x^3(1x)^3}\chi (x)\}$$
(31)
where

$$\xi (x)=1639+8039x15233x^2+10055x^3624x^4974x^5$$
(32)
$$\chi (x)=851341692x+64589x^260154x^3+99948x^4$$

$$112516x^5+45792x^6+(1808604x+53532x^2$$

$$75492x^328872x^4+109296x^555440x^6)ln2$$
(33)
The standard value of the gluonic condensate $`0(\alpha _s/\pi )G^20=0.012GeV^4`$ was taken in numerical calculations. Equation (30) was used and the effective instanton radius $`\rho _c`$ was chosen as $`\rho _c=0.5fm`$. This value is between the estimations of  ($`\rho _c=1fm`$) and  ($`\rho _c=0.33fm`$). (In the recent paper  it was argued that the liquid gas instanton model overestimates higher order gluonic condensates and, in order to correct this effect, larger values of $`\rho _c`$ comparing with  should be used). The Borel mass dependence of $`xu_T(x)`$ in the interval $`0.2<x<0.65`$ is weak at $`0.8<M^2<1.2GeV^2`$. Figure 3 shows $`xu_T(x)`$ at $`M^2=1GeV^2`$ and $`Q_0^2=4GeV^2`$. Dashed and thin solid lines demonstrate the influence of the variation in $`\rho _c`$ on the final result: the lower line corresponds to $`\rho _c=0.6fm`$ and the upper – to $`\rho _c=0.4fm`$. Our results are reliable at $`0.2<x<0.65`$, where $`d=4`$ and $`d=6`$ power corrections each comprise less than 30% of the bare loop contribution. (The contributions $`0|(\alpha _s/\pi )G^2|0`$ and $`0|g^3f^{abc}G_{\mu \nu }^aG_{\nu \lambda }^bG_{\lambda \mu }^c|0`$ are of the opposite sign and compensate one another, $`\alpha _s(M^2)\alpha _sa^2`$ contribution is negligible.) At $`\rho _c=0.4fm`$ the applicability domain shrinks to $`0.25<x<0.6`$.

The moments of quark distributions in the longitudinal $`\rho `$ meson are calculated in the same way, as it was done in for the case of pion: by matching with Regge behavior $`u(x)1/\sqrt{x}`$ at low $`x`$ and with quark counting rule $`u(x)(1x)^2`$ at large $`x`$. The matching points were chosen as $`x=0.10`$ at low $`x`$ and $`x=0.80`$ at large $`x`$. The numerical values of the moments for the longitudinally polarized $`\rho `$ are

$$_1^L=\underset{0}{\overset{1}{}}𝑑xu_\rho ^L(x)=1.06(1.05)$$

$$_2^L=\underset{0}{\overset{1}{}}x𝑑xu_\rho ^L(x)=0.39(0.37)$$
(34)
The values of moments, obtained by assuming that $`u(x)(1x)`$ at large $`x`$ are given in the parentheses.

Reliable calculation of moments for the case of transversally polarized $`\rho `$ meson is impossible, because of a narrow applicability domain in $`x`$ and expected double humps shape of the $`u`$-quark distribution, which does not allow soft matching with expected behaviour $`xu_\rho ^T(x)`$ at small and large $`x`$.

Now let us discuss the nonpolarizied $`\rho `$-meson case. The quark distribution function $`u(x)`$ in this case is equal to

$$u_\rho (x)=(u_\rho ^L(x)+2u_\rho ^T(x))/3$$

and we can determine $`u(x)`$ only in the region, where sum rules for $`u_\rho ^L(x)`$ and $`u_\rho ^T`$ are satisfied, i.e. in $`0.2\begin{array}{c}<\hfill \\ \hfill \end{array}x0.65`$. In this region, $`u_\rho (x)`$ is found to be very close to $`u_\pi (x)`$ (the difference in whole range of $`x`$ is no more than 10-15%).

## 5 Summary and discussion

Figure 4 gives the comparison of the valence $`u`$-quark distributions in the pion and longitudinally and transversally polarized $`\rho `$ mesons. The shapes of the curves are quite different, especially of $`xu_\rho ^T(x)`$ in comparison with $`xu_\rho ^L(x)`$ and $`xu_\pi (x)`$. Strongly different are also the second moments in the pion and longitudinal $`\rho `$-meson: the momentum fraction, carried by valence quarks and antiquarks $`(u+\overline{d})`$ in the longitudinal $`\rho `$ meson is about 0.8, while in the pion it is much less – about 0.4-0.5. All these differences are very large and many times larger than estimated uncertainties of our results. In the case of $`u`$-quark distribution in the pion the main source of them is the value of the renorminvariant quantity $`(2\pi )^4\alpha _s\overline{\psi }\psi 0^2`$. In our calculations, we took it to be equal to 0.13 GeV<sup>6</sup>. In fact, however, it is uncertain by a factor of 2. (Recent determination  of this quantity from $`\tau `$-decay data indicates that it may be two times larger). The perturbative corrections also introduce some uncertainties, especially at large $`x(x>0.6)`$ where the LO correction, which is taken into account, is large. The estimation of both effects shows that they may result in 10-20% variation (increase at $`x<0.3`$ and decrease at $`x>0.3`$) of $`xu_\pi (x)`$

For the $`u`$-quark distribution in the longitudinally polarized $`\rho `$ meson, the uncertainties in $`\alpha _s0\overline{\psi }\psi 0^2`$ do not play any role, because of higher $`M^2`$ values, so the expected accuracy is even better. One should note, that the accuracy of prediction of moments for longitudally polarizied $`\rho `$ meson are very high, because, as one can see from Fig.4, sum rules predict quark distribution almost in whole region of $`x`$ (from 0.1 up to 0.9). The fact that the first moment is so close to 1, also confirm that accuracy is high. Large value of the second moment mean, that valence quarks in longitudally polarizied $`\rho `$ meson carried about 0.8 of the total momentum, so one can conclude that the total (gluon and quark) sea is strongly suppressed. Moreover, this prediction for the second moment is close to those obtained in  from quite another sum rules, and also agreed with the lattice calculation results , so the sea supression in longitudally polarizied $`\rho `$ meson can be treated as a theoretically well established fact.

The accuracy of our results for the $`u`$-quark distribution in the transversally polarized $`\rho `$ meson is lower, because of a large role of $`d=4`$ and $`d=6`$ gluonic condensate contributions. For the latter, as was discussed before, there are only estimations, based on the instanton model, and the gluon condensate $`0(\alpha _s/\pi )G^20`$ is also uncertain by a factor 1.5. Estimations show that they result in not more than 30-40% variation in $`xu_\rho ^T`$ at $`x0.30.4`$, but in much less variation at $`x0.50.6`$. (The LO perturbative corrections are not more than 20% at small $`x`$ and negligible at large $`x`$). So one can see, that the difference, obtained in the quark distribution in the pion, longitudinally and transversally polarizied $`\rho `$-meson in fact are much more larger than any possible uncertainties of the results.

In summary the main physical conclusion are the following:

(1) The sea is strongly supressed in the longitudinally polarizied $`\rho `$-meson.

(2) The quark distribution in the polarizied $`\rho `$-meson are significantly dependend on polarization.

(3) The quark distributions in the the pion and $`\rho `$-meson have a not too much in common. The specific properties of the pion, as a Goldstone boson, manifest themselves in the quark distributions different from those in the $`\rho `$ meson. $`SU(6)`$ symmetry may, probably, take place for the static properties of $`\pi `$ and $`\rho `$, but not for their inner structure. We have no explanation, why $`u`$-quark distributions in the pion and unpolarized $`\rho `$-meson at $`0.2<x<0.65`$ are close to one another – whether it is a pure accident or there are some deep reasons for it.

This study was supported in part by Award no. RP2-2247 of U.S. Civilian Research and Development Foundation for Independent States of Former Soviet Union (CRDF), by the Russian Foundation of Basic Research, project no. 00-02-17808 and by INTAS Call 2000, project 587.