text:
# Relevance Vector Machines for classifying points and regions in biological sequences.

## 1 Background

Machine learning technologies have played an important role in bioinformatics for some years now, both for extracting information from existing data sets and for making predictions of interesting features in previously unannotated data . Important examples include Hidden Markov Models for annotating protein domains  and neural networks for predicting protein secondary structure . Here, we describe a family of new machine learning methods for classifying biological sequences based on the occurrence of particular motifs within the sequence. All our methods share some common features: they can all be represented as generalized linear models and trained using the Sampling Relevance Vector Machine strategy, and they all use position-weight matrices (PWMs) as the basic sensors to detect motifs in the sequence data. However, differences in the way that that PWM scores are presented to the linear model provide us with several learning methods suited to different purposes.

In this paper, we first review the basic Relevance Vector Machine machine method, and describe some extensions which make the RVM easier to apply to complex problems. We then describe three different RVM-sequence models, one for analyzing individual points in sequences (e.g. transcription start sites, splice sites) and two for analyzing larger regions of sequence. These methods have been applied to a range of interesting sequence analysis questions, and we highlight some applications throughout the paper.

### 1.1 Generalized Linear Models

Generalized linear models (GLMs) are a commonly used form of model for both classification (separate data into two or more classes) and regression (estimate the value of a continuous function) problems. GLMs take the form:

$$\eta (𝐱)=\underset{m=1}{\overset{M}{}}\beta _m\varphi _m(𝐱)+K$$
(1.1)
where $`\varphi `$ is a set of $`M`$ basis functions (which can be arbitrary real-valued functions) and $`\beta `$ is a vector of weights. While, in the machine learning literature, the data items are usually vectors in some data space, since the model only refers to the data via the basis functions, any type of data can be used so long as appropriate basis function are provided. The classical view of ‘training’ a GLM is simply to find some weights which offer a good fit to the provided training data. Since in this case we are considering classification problems, we apply the logistic link function to the model output:

$$\pi (𝐱)=\frac{1}{1+e^{\eta (𝐱)}}$$
(1.2)
This renormalizes the model output such that $`0\pi (𝐱)1`$, and can be interpreted as a probability that $`𝐱`$ is a member of the ‘positive’ class for the classification problem in hand. The classical problem of learning a set of weights, $`\beta `$, such that the model optimally classifies some set of labeled training data has been called logistic regression, and can be solved using techniques such as Iterated Reweighted Least Squares .

One way of looking at generalized linear models is that the basis functions define a projection of the data into a high-dimensional space (called feature space) where the data is either linear (for regression problems) or linearly separable (for classification problems). Clearly, the choice of projection is important. An important step in GLM learning is to find a feature space which allows a linear model to fit the training data, while not being of such high dimensionality that overfitting becomes a problem. Processes for solving this problem are described as feature selection.

### 1.2 The Relevance Vector Machine

The Relevance Vector Machine (RVM)  is a sparse method for training generalized linear models. This means that it will generally select only a subset (often a small subset) of the provided basis functions to use in the final model. Sparsity is generally considered a desirable feature in a machine learning system . This is consistent with the intuition that a simple model is more likely to make useful generalizations which can be applied to unseen data, rather than solving the problem ‘trivially’ by memorizing the training data. Sparse training methods also offer a partial solution to the question of feature selection. The RVM is named by analogy to the better-known Support Vector Machine method, which is also a kind of sparse GLM trainer, and indeed the RVM was initially presented as an alternative and direct competitor to the SVM. However, SVMs can only be applied to training a restricted subset of generalized linear models – those that can be defined by a suitable kernel function – while the RVM can train a GLM with any collection of basis functions. The ability to perform sparse training with arbitrary basis functions opens up many new possibilities, and is crucial for all the sequence analysis methods described below.

Briefly, the Relevance Vector Machine is a Bayesian approach for training a linear model. In a binary classification problem, where each training datum $`𝐱_n`$ has a label $`t_n`$ (either $`0`$ or $`1`$), we can represent the probability that the data set is correctly labeled given some classifier model $`\pi (𝐱)`$ as:

$$P(t|𝐱,\beta )=\underset{n=1}{\overset{N}{}}\pi (𝐱_n)^{t_n}(1\pi (𝐱_n))^{1t_n}$$
(1.3)
Where $`\beta `$ is the weights vector of our classification GLM. If we assume that the training data is correctly labeled, Bayes’ theorem allows us to turn this expression around and infer likely values of the weights given some labeled data.

$$P(\beta |𝐱,t)P(\beta )P(t|𝐱,\beta )$$
(1.4)
In this expression, we have introduced an extra probability distribution, $`P(\beta )`$, which is our prior belief in the values of the weights. If we merely wished to perform classical GLM training, we would have no preference for any particular value of the weights, and would encode this by providing a very broad (“non-informative”) prior. However, as discussed previously, we believe that simple models are more likely to make useful generalizations. A preference for simplicity can be encoded using an Automatic Relevance Determination prior . In this case, we introduce an additional vector of parameters, $`\alpha `$. Each element of the $`\alpha `$ vector controls the width of the prior over the corresponding weight:

$$P(\beta )=\underset{m=1}{\overset{M}{}}𝒢(\beta _m|0,\alpha _m^1)$$
(1.5)
To include these new $`\alpha `$ parameter in the inference process, we also need to specify a hyperprior over values of $`\alpha `$. For the RVM, a very broad gamma distribution is used – a standard choice for non-informative priors over Gaussian width parameters – but the precise choice is of little importance so long as it is sufficiently broad.

Considering just a single basis function, there are two possibilities:

* The basis function provides additional information about the specified classification problem. When its weight is set to some non-zero value, the amount of misclassified training data is reduced. This increases the value of equation 1.3, and therefore the probability of that model given the data.
* If the basis function provides no information because it is irrelevant to the problem, there is no value of the weight that will lead to a significant increase in the likelihood. At this point, the prior term in the model comes into play: by setting the $`\alpha _m`$ parameter to a large value, the prior distribution $`P(\beta _m)`$ becomes sharply peaked around zero. By then setting $`\beta _m`$ to zero, the posterior probability of the model is maximized. Similarly, when two basis functions offer redundant information, the posterior is maximized by using only one of them in the model.

When a basis function has a sufficiently high $`\alpha `$, it can be marked as irrelevant, and removed from the model. As a result, the RVM will learn simple models even when presented with a large starting set of basis functions.

In this discussion, we have neglected the practical issues of evaluating equation 1.4 to find good values for the model parameters. A number of strategies can be used, including Monte-Carlo simulation, maximum likelihood estimation , and variational inference . While we have experimented with all these methods, the training library which we use used when developing all the models described here is based on the variational approximation given in . This training method involves iterative update of the model parameters, and it is possible to remove irrelevant basis functions during the course of the training, so the training process runs progressively faster as it converges towards the optimal solution (and corresponding minimal set of basis functions).

### 1.3 Sampling-RVM

In principle, we would like to provide as little prior knowledge as possible about which feature spaces to consider when solving a problem, instead allowing a sparse trainer such as the RVM to pick freely from a wide range of possible basis functions. Unfortunately, computational costs rapidly become significant. The progressive simplification of the problem means that it is difficult to quantify the full computational complexity since the size of the matrices varies from cycle to cycle. However, scaling is quite substantially worse than linear.

Here, we suggest a pragmatic approach to handling large sets of candidate basis function. A working set is first initialized with a subset of basis functions, picked at random from the complete pool of candidates. The trainer is run as previously described, and as before some $`\alpha `$ values increase to a level such that it is possible to remove the associated basis functions from further consideration. Once the size of the working set drops below a designated ‘low water mark’, additional basis functions are added from the pool. At this point, all $`\alpha `$ and $`\beta `$ values are reinitialized and training continues. In this way, the working set fluctuates between high and low water marks until the pool is exhausted, at which point the trainer continues to run until the weights and priors no longer change significantly between cycles (i.e. convergence) to give a complete model.

To evaluate the performance benefits of this method, we considered a simple classification problem in Cartesian 2-space. This involved two equal-sized classes of points sampled from two Gaussians, with a substantial space between them. Basis functions were generated in an “SVM-like” manner, with a radial basis function centred on each point in the training data. Thus, the size of the basis-value matrix of the RVM trainer increases with the square of dataset size. As shown in figure 1, an optimal model requires only two basis functions, and gives 100% classification accuracy on the training data.

A range of problem sizes were tested, from 50 points in each class (i.e. 100 basis functions) to 300 points. We compared the full-set training method with the incremental method, using a high water mark of 20 basis functions and a low water mark of 15. Timings are shown in figure 2. For small problems, the incremental training approach is in fact significantly slower, due to the increased number of cycles and the need to periodically restart the training process as more basis functions are added. However, the scalability of the method to large problems is substantially better.

It is interesting to note that, while in all cases all the training data was correctly classified, the full-set training runs with 300 points did not learn the simple 2-basis function model seen with smaller training sets, but instead appeared to converge before all unnecessary basis functions had been removed. However, the incremental trainer consistently learned a fully sparse model for the same dataset with no complications. We believe that this indicates the limits of either the linear algebra routines used by the trainer, or the precision of the standard ‘double’ floating point mathematics available on most computers. Errors caused by lack of precision (often described as numerical instability) are a common problem in numerical computation, and large linear algebra computations are a frequent source of trouble. By keeping the size of individual computations down, using the small-working-set variant of the RVM neatly sidesteps these problems and makes it possible to tackle very large problems efficiently.

Having made the step toward the approach of gradually sifting through a large set of basis functions, it is only a small additional step to consider an implementation which generates new basis functions on the fly during the training process. In principle, this can allow exploration of infinite sets of basis functions. This is particularly feasible in cases where there are clear correlations between basis functions. For instance, a radial basis function with centre (10.0, 10.2) and variance 10 will have outputs that are highly correlated to a second radial basis function with centre (10.0, 10.0) and equal variance. This means that if the first function is found to be informative for modeling a particular problem, it is likely that the second will also be effective. Therefore, after some period of training, it is likely to be worthwhile proposing more candidate basis functions which are correlated to those currently in the working set.

In many practical applications, including the sequence models described here, a significant proportion of total run time may be spent evaluating the basis functions. However, since the evaluation of each basis function for each item in the training set is an independent operation, this part of the process is easily parallelized on both individual multiprocessor computers and clusters.

### 1.4 Position-weight matrices

When studying biological sequences, it is often natural to look for significant short motifs. In many cases – for example, gene promoters and exonic splice enhancers – these motifs act as preferred binding sites for specific proteins. It is rare, however, that there is only one sequence which will bind a particular protein: normally, there are at least some positions in the motif where more than one base is acceptable. Given a large number of examples, a good way to represent knowledge about a motif is a Position Weight Matrix (PWM) . This is a matrix where each element represents the probability of a given nucleotide occurring at a particular position in the sequence. In other words, each column of the matrix is a probability distribution over the DNA alphabet. A PWM can be viewed as a probabilistic model of a fixed-length sequence:

$$W(s)=\underset{i=1}{\overset{|W|}{}}W_i(s_i)$$
(1.6)
Where $`|W|`$ is the length of the PWM, and $`W_i`$ is the probability distribution represented by the $`i`$’th column of the matrix. The higher the probability of a given sequence under this model, the more similar the sequence is to the “ideal” motif, and therefore the higher the probability that it will function as a binding site. This is the zeroth order model: each position in the sequence is assumed to be independent of all the others.

It should be pointed out that the assumption of bases being independent in a motif is unlikely to be true in most cases. When proteins bind to DNA, they often significantly deform it, so mechanical properties of the DNA other than the actual base sequence will be significant – in particular, the flexibility of the double helix. This property is determined largely by interactions between neighbouring base pairs. Models which take these, and perhaps also longer range interactions, into account can be expected to better predict the binding of a protein to a given sequence . However, since non-independent models have many more parameters than simple PWMs, they require more example sequences to learn effectively.

## 2 Eponine RVM models

### 2.1 Eponine Anchored Sequence

The Eponine Anchored Sequence (EAS) model is a system for analyzing the sequence around some particular interesting point (the anchor point). A number of important sequence analysis questions, such as prediction of transcription start sites and splice sites, are best expressed in terms of classifying individual points in the sequence. In practice, the data is presented to this model as a large piece of sequence data and an integer defining the anchor point under consideration. When searching for features in bulk sequence, the same sequence is presented many times while scanning the anchor point along its length. The basic element of an EAS model is the positioned constraint (hereafter, PC). This consists of:

* A preferred sequence motif, defined as a position-weight matrix.
* A probability distribution over integer offsets relative to the anchor point, which defines the expected localization of the motif. In the work presented here, these distributions are always discretized Gaussians (i.e. the result of integrating the Gaussian probability density function over unit intervals). Gaussians were chosen because of their familiarity, and a smooth shape that made Gaussian-based models less prone to overfitting than functions with abrupt changes, such as square waves. However, any distribution over integers could in principle be used here.

To obtain a score for a PC on a given piece of data, the program scans over all positions in the sequence which are assigned a non-infinitesimal probability by the chosen position distribution. For each position, the probability of the sequence motif starting at that position being emitted by the chosen weight matrix is evaluated. The final score is given by:

$$\varphi (C)=\frac{1}{|W|}\mathrm{log}\underset{i=\mathrm{}}{\overset{\mathrm{}}{}}P(i)W(C,i)$$
(2.1)
where $`C`$ is a genomic context, $`P`$ is a position distribution, and $`W(C,i)`$ is a DNA weight matrix probability for offset $`i`$ relative to the anchor point of $`C`$. Note the division by $`|W|`$, the number of columns in the weight matrix (i.e. the length of the sequence motif which it defines). This is important, since this method allow motifs with a wide range of lengths – with the trainer implementation described here, the length varies between 2 and 20 columns. The RVM trainer has a weak bias towards selecting basis function with higher absolute magnitudes. Normalizing the scores allows unbiased selection between motifs of different lengths. It is somewhat analogous to the whitening process often used to pre-process data for SVM classifiers, where all the training vectors are normalized to constant length .

A single PC describes an individual sequence motif and its relationship to a point in a sequence, but a set of them can be combined to describe more complex structures. Figure 3 shows a schematic of a model combining three positioned constraints. If the final output score is defined as a weighted sum of individual PC scores, the combined model is a generalized linear model over genomic contexts, with the PCs as basis functions. Therefore, it is possible use the Relevance Vector Machine method to reduce a large set of candidate PCs down to a sparse model containing a small, informative subset.

While this description of the model architecture emphasizes detection of single, well-defined “words” in the sequence, since the overall PC score is based on a sum of weight matrix scores across a region, it is also possible to represent general compositional biases of a region by picking a PC with a short weight matrix and a very broad position distribution. It is also quite acceptable for the distributions of two or more PCs to overlap. Since arbitrary weights are assigned in the training process to fit the model to the training data, any issues with double-counting of a particular piece of information are corrected automatically

The space of potentially interesting PCs is extremely large. Even taking a highly simplified view, restricting the constraints to simple motifs rather than weight matrices and the position distributions to Gaussians of a constant width, there are over one million PCs representing six-base motifs with distributions centred at positions in the range \[-250:50\] relative to the anchor point. This space it too large to search exhaustively. Of course, in practice the EAS framework allows for an infinite (limited only by numerical precision on the probability values) number of PCs. Fortunately, as in the case of radial basis functions, the space of possible basis functions for this model is highly correlated. Making a small change to, say, one of the probabilities in a weight matrix will give a second PC whose output on a given sequence is correlated with the first. For this reason, exploring regions of a conceptual “PC-space” in the neighborhood of constraints which have already proved to be informative is likely to reveal even more informative constraints.

Taking advantage of these correlations, EAS models can be trained using an RVM with a sampling strategy to create new basis functions. When the size of the working set falls below the low water mark, the trainer selects a sampling strategy at random from the following set:

* Constructing a new PC, not based on the current set. This is performed by the following algorithm:

1. First, select a context at random from the training set (either positive or negative, without bias).
    2. Pick some point relative to that context’s anchor point.
    3. From that point, take a sequence motif of between 3 and 6 bases in length, and construct a weight matrix which optimally matches that consensus sequence, but includes some degree of uncertainty.
    4. Construct a PC using the newly selected weight matrix with a Gaussian position distribution of random width, centred at the position at which the motif was originally found.

Obviously, PCs selected in this way will strongly match the training example from which they were originally derived. This is closely analogous to selecting radial basis functions centred on points in the training data set.
* Selecting an existing PC and adjusting the emission spectrum of one column of its weight matrix, by sampling from a tightly-focused Dirichlet distribution centered on the current values.
* Adding an extra column to either the start or the end of an existing weight matrix, up to a maximum number of columns (in this case 20).
* Removing the start or end column from an existing weight matrix, down to a minimum of two columns
* Adjusting the width parameter of a Gaussian position distribution
* Adjust the centre position for a Gaussian position distribution

For an initial period of 200 cycles, only the first sampling rule (creation of novel basis functions) was used. After this point, the full range of sampling strategies were available, and training proceeded by a mixture of sampling and introduction of novel PCs.

This model was originally developed for the purposes of transcription start site prediction, and has been shown to perform well for that purpose, at least in mammalian genomes . Subsequently, it has been used for predicting transcription termination sites \[A. Ramadass, personal communication\].

### 2.2 Eponine Windowed Sequence

The Eponine Windowed Sequence (EWS) model is designed by analogy to the EAS model, but rather than targeting individual points in the sequence, it is designed to classify small regions or windows of a sequence, based purely on their own sequence content. A naturally windowed classification problem would be to distinguish exons from intergenic sequence. Once again, the sensors are position-weight matrices, but each weight matrix is scanned across the entire region.

Normalizing for the length of the sequence being inspected and the size of the PWM, the basis functions of the model take the form:

$$\varphi (S)=\frac{4^{|W|}}{|S||W|+1}\underset{i=1}{\overset{|S||W|+1}{}}W(S_i^{i+|W|1})$$
(2.2)
where $`W(s)`$ is the probability that sequence $`s`$ was emitted by weight matrix $`W`$, and $`S_i^j`$ denotes a subsequence from $`i`$ to $`j`$.

An initial set of basis functions is proposed by taking all possible DNA motifs of a specified length (typically 5) and generate weight matrices which recognize these motifs. When new basis functions are needed, one of the following sampling moves is used to generate a variant of a basis function currently in the working set.

* Generate a new weight matrix in which each column is a sample from a Dirichlet distribution with its mode equal to the weights in the corresponding column of the parent weight matrix.
* Generate a new weight matrix one column shorter than the parent by removing either the first of the last column.
* Generate a new weight matrix with an extra column at either the start or the end, biased in favour of a random base.

We have used the EWS model to explore common motifs in regions of non-coding sequence conserved between the mouse and human genomes, as discussed in . We have also successfully trained this model on windows of sequence upstream of annotated genes, in which case it works successfully as a promoter predictor, although the EAS model is generally superior for this purpose since it is able to make predictions of the exact locations of transcription start sites, rather than merely highlighting broad regions likely to have promoter activity.

### 2.3 Convolved Eponine Windowed Sequence (C-EWS)

The C-EWS model is an extension of EWS to capture larger-scale patterns in sequences. In this case, each basis function defines a ‘scaffold’ consisting of one or more PWMs, each with an associated position distribution relative to a scaffold anchor. In principle, distributions such as the discretized Gaussian extend to infinity, but in practice it is reasonably to apply some cut-off: for instance, only considering the portion of the distribution which includes 99% of the total probability mass. The probabilities of all points outside this region are assumed to be infinitesimal and ignored. Now that the distributions have finite size, for a given scaffold there is a pair of integers, $`n`$ and $`m`$, such that when the scaffold anchor is placed in the interval \[$`n`$ : $`m`$\], the non-infinitesimal parts of all the position distributions fall entirely within the length of a particular target sequence. A score for scanning this scaffold across the sequence can be given by

$$\varphi (S)=Z\underset{i=n}{\overset{m}{}}\left(\underset{k=1}{\overset{K}{}}\left(\underset{j=\mathrm{}}{\overset{\mathrm{}}{}}P_k(j)W_k(S_{i+j}^{i+j+|W_k|})\right)\right)$$
(2.3)
where $`P_k`$ is the $`k`$’th position distribution and $`W_k`$ is the $`k`$’th weight matrix in the scaffold. $`Z`$ is the normalizing constant:

$$Z=\frac{4^{_{k=1}^K|W_k|}}{mn+1}$$
(2.4)
For the case when the scaffold only contains one PWM with a narrow distribution, the results will be the same as those from equation 2.2. So this can be considered a direct extension to the original EWS model that can capture information about sets of motifs with correlated positions. Since the scaffold scores are only evaluated in the regions where the whole scaffold fits onto the sequence window, there is a risk of introducing edge effects. A possible future solution to this would be to use windows that are a little larger than the actual region, of interest, and use “soft boundaries” where scores from the edges of the window are given less weight than those at the centre. To train scaffold-based models, all the sampling rules for the EWS model are used, together with some additional rules for constructing and adjusting scaffolds:

* Combine the sets of motifs from two scaffolds, with randomly chosen offsets between the two (up to some maximum number of weight matrices per scaffold).
* Take a scaffold with two or more PWMs and return the scaffold with one of those PWMs (picked at random) removed
* Alter the position or width of one of the relative position distributions in a scaffold.

To date, we have only investigated scaffolds with small numbers of PWMs (up to a maximum of 3). If more complex scaffolds were used, it might be appropriate to consider additional sampling rules which combined a subset of PWMs from two parent scaffolds: this would be analogous to the use of ‘cross-over’ rules in the field of genetic algorithms.

The C-EWS model can be applied in any situation where EWS can be used. It should never give significantly worse results, but is better able to capture longer-range features in a sequence: for instance binding sites for two separate proteins which interact with one another (and therefore have binding sites close to one another). This model scheme has recently been applied successfully to the discovery of exonic splice enhancer motifs in coding exon sequence \[B. Leong, manuscript in preparation\].

## 3 Conclusions

The Sampling-RVM learning systems presented here are powerful strategies for analyzing sequence data. All these methods have been successfully applied in the past and we hope they will be applied more widely and on new systems in the future. We continue to develop and experiment with these methods, with a particular interest in promoter and transcription start site research.

It should be noted that in all cases, it is in principle possible to use any of these models to distinguish between more than two classes of object – for example, a set of promoter training data could be split into several classes based on the tissues in which they are active. However, the variational RVM trainer which we have used here relies on an approximation which is only valid when considering two-class problems. We have developed a prototype RVM implementation which uses Metropolis-Hastings sampling to perform the Bayesian inference, and this has successfully completed three-way classification problems, albeit at the cost of increased computational requirements.

A key part of the strength of these approaches is the use of the Relevance Vector Machine as an extremely flexible method for training sparse models. The RVM has already been applied in the classification of microarray gene expression data pocock.thesis, and we encourage its use in other areas of bioinformatics and beyond.

## 4 Availability

The RVM training library and Eponine sequence models are implemented in the Java programming language, using the COLT and BioJava libraries. Java source code is available to all interested researchers on request from TD.

## 5 Acknowledgments

Thanks to Matthew Pocock for many interesting discussions on machine learning in bioinformatics. TD would like to thank the Wellcome trust for funding.


format: text

text:
# 1 Introduction: the fate of spin-orbit interaction for heavy nuclei and the mechanism of SCHNC

## 1 Introduction: the fate of spin-orbit interaction for heavy nuclei and the mechanism of SCHNC

Heavy nuclei are strongly absorbing targets. Whether the spin-flip effects in high energy scattering are washed out by this absorption or not is not an obvious issue which we address in this communication on an example of diffractive vector mesons.

The standard argument for vanishing of the spin-flip in elastic scattering of spin $`\frac{1}{2}`$ particles off strongly absorbing nuclei goes as follows: In the presence of the spin-orbit interaction the scattering amplitude $`f=A+2B\widehat{𝐬}𝐧`$, where $`\widehat{𝐬}`$ is the spin operator, $`𝐧`$ is the normal to the scattering plane and the partial wave expansion of the helicity-non-flip, $`f_0`$, and the helicity-flip, $`f_1`$, amplitudes reads

$`f_0`$ $`=`$ $`{\displaystyle \frac{1}{2ip}}{\displaystyle \underset{l}{}}\{(l+1)[\mathrm{exp}(2i\delta _l^+)1]+l[\mathrm{exp}(2i\delta _l^{})1]\}P_l(\mathrm{cos}\theta ),`$
$`f_1`$ $`=`$ $`{\displaystyle \frac{1}{2ip}}{\displaystyle \underset{l}{}}[\mathrm{exp}(2i\delta _l^+)\mathrm{exp}(2i\delta _l^{})]P_l^1(\mathrm{cos}\theta ),`$ (1)
where $`\delta _l^\pm `$ are the scattering phases for $`j=l\pm \frac{1}{2}`$. In the presence of strong absorption the scattering phases acquire large imaginary parts, $`\mathrm{exp}(2i\delta _l^\pm )0`$, and, consequently, for the momentum transfer $`𝚫`$ within the diffraction cone $`f_1/f_00`$. More detailed treatment for elastic scattering of protons off nuclei is found in Glauber’s lectures , the net result is that the small contribution to the spin-flip amplitudes comes only from the periphery of the nucleus, so that $`f_1/f_0A^{1/3}`$, where $`A`$ is the mass number of a nucleus.

In high-energy QCD the diffractive production of vector mesons,

$$\gamma ^{}pVp^{},$$

proceeds via the exchange of colorless system of gluons in the $`t`$-channel.

The fundamental property of such (multiple) gluon exchange is an exact conservation of the s-channel helicity of quarks, recall the similar property in the QED scattering of electrons in the Coulomb field . None the less, QCD predicts a non-vanishing helicity flip in diffractive production of vector mesons off unpolarized nucleons : the origin of this SCHNC is in the subtle possibility that the sum of helicities of the quark-antiquark pair in the diagrams of fig. (1) can be unequal to the helicity of photons and vector mesons. In the non-relativistic case the pure S-wave deuteron with spin up consists of the spin-up proton and neutron. However, in the relativistic case while the longitudinal virtual photon contains the $`q\overline{q}`$ pair with $`\lambda _q+\lambda _{\overline{q}}=0`$, the transverse photon with helicity $`\lambda _\gamma =\pm 1`$ besides the $`q\overline{q}`$ state with $`\lambda _q+\lambda _{\overline{q}}=\lambda _\gamma =\pm 1`$, also contains the state with $`\lambda _q+\lambda _{\overline{q}}=0`$, in which the helicity of the photon is carried by the orbital angular momentum in the $`q\overline{q}`$ system. Furthermore, it is precisely the state $`\lambda _q+\lambda _{\overline{q}}=0`$ which gives the dominant contribution to the absorption of transverse photons and the proton SF $`F_{2p}(x,Q^2)`$ in the Bjorken limit. From the point of view of the vector meson production, it is important that the transverse and longitudinal $`\gamma ^{}`$ and $`V`$ share the intermediate $`q\overline{q}`$ state with $`\lambda _q+\lambda _{\overline{q}}=0`$, which allows the s-channel helicity non-conserving (SCHNC) transitions between the transverse (longitudinal) $`\gamma ^{}`$ and longitudinal(transverse) vector meson V. As a matter of fact, this mechanism of SCHNC does not require an applicability of pQCD, B.G. Zakharov was the first to introduce it in application to nucleon-nucleon scattering . The theoretical predicition of energy-independent SCHNC in diffractive vector meson production has been confirmed experimentally at HERA , the detailed comparison of the theory and experiment is found in .

Now notice that the argument about exact SCHC of quarks and antiquarks applies equally to a one-pomeron exchange in the scattering off a free nucleon and to multiple pomeron exchange in the scattering off a nuclear target. Then for a sufficiently high energy such that the lifetime of the $`q\overline{q}`$ fluctuation of the photon, often referred to as the coherence time , and the formation time of the vector meson are larger than the radius of the nucleus $`R_A`$ (for instance, see ), the above described origin of SCHNC must be equally at work for the nuclear and free-nucleon targets. In this communication we expand on this point, from the practical point of view one speaks of the values of

$$x=\frac{Q^2+m_V^2}{2\nu }<x_A\text{ }<0.1A^{1/3}.$$
(2)
Another property of interactions with nuclei is the so-called saturation scale $`Q_A`$ which for partons with $`x\text{ }<x_A`$ defines the transverse momentum below which their density is lowered by parton fusion effects . It is interesting to see how the emergence of the saturation scale affects the $`Q^2`$-dependence of diffractive vector meson production, in particular, its SCHNC properties. Here one must compare the saturation scale $`Q_A^2`$ to the usual hard scale for diffractive vector meson production

$$\overline{Q}^2\frac{1}{4}(Q^2+m_V^2).$$
(3)
The further presentation is organized as follows. In section 2 we start with the formulation of the color-dipole approach to calculation of helicity amplitudes for diffractive vector meson production. In section 3 we introduce the scanning radius and comment on the potential sensitivity to the short-distance behaviour of the wave function of vector mesons. In section 4 we start a discussion of nuclear effects on an example of coherent diffraction $`\gamma ^{}AVA`$ when the recoil nucleus remains in the ground state. We demonstrate how the dependence $`\overline{Q}^4`$ which is common to all helicity amplituides changes to $`\overline{Q}^2Q_A^2`$ for $`\overline{Q}^2\text{ }<Q_A^2`$. Even stronger impact of saturation is found for incoherent diffractive vector mesons considered in section 5. In the Conclusions we summarize our principal findings and comment on the possibilities of the COMPASS experiment at CERN.

## 2 The free nucleon target

For the purposes of our discussion it is convenient to resort to the color dipole formalism: the production process depicted in fig. (1) factorizes into splitting of the photon into $`q\overline{q}`$ dipole way upstream the target, s-channel helicity conserving elastic scattering of the dipole off a target, and projection of the $`q\overline{q}`$ dipole onto the vector meson state. We restrict ourselves to the contribution from the $`q\overline{q}`$ Fock states of the vector meson which is a good approximation for $`xx_A`$. The momentum-space calculation of the helicity amplitudes has been worked out time ago in , the crucial ingredient in preserving the rotation invariance is the concept of the running polarization vector for the longitudinal vector mesons. Following , we make the Fourier transform to the color-dipole space and represent the helicity amplitudes $`𝒜_{fi}(x,𝚫)`$, where $`i=\lambda _\gamma `$ and $`f=\lambda _V`$ are helicities of the initial state photon and the final state vector meson, respectively, in the color dipole factorization form

$$𝒜_{fi}(x,𝚫)=V_f|𝒜_{q\overline{q}}(𝐫,𝚫)|\gamma _i^{}=i_0^1𝑑zd^2𝐫\sigma (𝐫,𝚫)\mathrm{exp}[\frac{i}{2}(12z)(𝐫𝚫)]I_{fi}(z,𝐫),$$
(4)
where $`𝚫`$ is the transverse momentum transfer in the $`\gamma ^{}V`$ transition, $`I_{fi}(z,𝐫)=\mathrm{\Psi }_{V,f}^{}(z,𝐫)\mathrm{\Psi }_{\gamma ^{},i}(z,𝐫)`$ and the summation over the helicities $`\lambda ,\overline{\lambda }`$ of the intermediate $`q\overline{q}`$ pair is understood. The wave function $`\mathrm{\Psi }_{V,f}(z,𝐫)`$ of the final state vector meson contains the spin-orbital part and the ”radial” wave functions, defined in terms of the vertex function $`\mathrm{\Gamma }_V(z,𝐤)\overline{q}S_\mu qV_\mu `$, where $`S_\mu `$ is the relevant Dirac structure and $`V_\mu `$ is the running polarization vector which must be so defined as to guarantee the rotational invariance for the fixed invariant mass $`M`$ of the lightcone $`q\overline{q}`$ Fock state of the vector meson ,

$$M^2=\frac{m_f^2+𝐤^2}{z(1z)},$$

where $`𝐤`$ is the transverse momentum of the quark in the vector meson and $`z`$ and $`(1z)`$ are fractions of the lightcone momentum of the vector meson carried by the quark and antiquark, respectively. In the momentum-space $`\psi _V(z,𝐤)\mathrm{\Gamma }_V(z,𝐤)/(M^2m_V^2)`$, the Fourier transform to the dipole space depends on the Dirac structure $`S_\mu `$. For the sake of simplicity, here we take $`S_\mu =\gamma _\mu `$, the exact form for the pure $`S`$ and $`D`$ wave states is found in , the major conclusions on the impact of nuclear absorption do not depend on the exact form of $`S_\mu `$. If one defines the radial wave functions for the transverse (T) and longitudinal (L) vector mesons as

$$\psi _T(z,𝐫)=d^2𝐤\psi (z,𝐤)\mathrm{exp}(i\mathrm{𝐤𝐫}),\psi _L(z,𝐫)=d^2𝐤M\psi (z,𝐤)\mathrm{exp}(i\mathrm{𝐤𝐫}),$$
(5)
then

$`I_{LL}`$ $`=`$ $`4Qz^2(1z)^2K_0(\epsilon r)\psi _L(z,r)`$ (6)
$`I_{TT}`$ $`=`$ $`m_f^2K_0(\epsilon r)\psi _T(z,r)[z^2+(1z)^2]\epsilon K_1(\epsilon r)\psi _T^{}(z,r)`$ (7)
$`I_{LT}`$ $`=`$ $`i2z(1z)(12z)\psi _L(z,r)\epsilon K_1(\epsilon r){\displaystyle \frac{(\mathrm{𝐞𝐫})}{r}}`$ (8)
$`I_{TL}`$ $`=`$ $`i2Qz(1z)(12z)K_0(\epsilon r)\psi _T^{}(z,r){\displaystyle \frac{(𝐕^{}𝐫)}{r}}`$ (9)
$`I_{TT^{}}`$ $`=`$ $`4z(1z)\epsilon K_1(\epsilon r)\psi _T^{}(z,r){\displaystyle \frac{(\mathrm{𝐞𝐫})^2}{r^2}}`$ (10)
where $`\psi _T^{}(z,r)=\psi _T(z,r)/r`$, the polarization vectors $`𝐞`$ and $`𝐕`$ are for the transverse photon and vector meson, respectively, $`\epsilon ^2=z(1z)Q^2+m_f^2`$, the Bessel functions $`K_0(x)`$ and $`K_1=K_0^{}(x)`$ describe the lightcone wave function of the photon, $`I_{TT}`$ and $`I_{TT^{}}`$ describe the helicity-non-flip and double-helicity-flip production of transverse vector mesons by transverse photons, in the latter case $`𝐕^{}=𝐞`$.

The off-forward generalization of the color-dipole scattering amplitude has been introduced in

$`\sigma (𝐫,𝚫)={\displaystyle \frac{2\pi }{3}}{\displaystyle \frac{d^2𝜿}{(𝜿\frac{1}{2}𝚫)^2(𝜿+\frac{1}{2}𝚫)^2}(x,𝜿,𝚫)\alpha _S(𝜿^2)}`$
$`\times \left\{[1\mathrm{exp}(i𝜿𝐫)][1\mathrm{exp}(i𝜿𝐫)][1\mathrm{exp}({\displaystyle \frac{1}{2}}i𝚫𝐫)][1\mathrm{exp}({\displaystyle \frac{1}{2}}i𝚫𝐫)]\right\}`$ (11)
Here $`(x,𝜿,𝚫)`$ is the off-forward unintegrated differential gluon structure function of the nucleon, the gross features of its $`𝚫`$-dependence are discussed in .

The analysis of  focused on the $`LL`$ and $`TT`$ amplitudes, in which case the dominant contribution comes from $`z\frac{1}{2}`$ and corrections to the $`𝚫`$-dependence from the factor $`\mathrm{exp}[\frac{i}{2}(12z)(𝐫𝚫)]`$ in (4) can be neglected. This factor is crucial, though, for the helicity-flip transitions. For small dipoles and within the diffraction cone the leading components of the $`LT`$ and $`TL`$ amplitudes come from the second term in the expansion

$$\mathrm{exp}\left(\frac{1}{2}i(12z)(𝚫𝐫)\right)=1+\frac{1}{2}i(12z)(𝚫𝐫)$$
(12)
so that upon the azimuthal averaging the effective integrands take the form

$`I_{LT}`$ $`=`$ $`{\displaystyle \frac{1}{2}}z(1z)(12z)^2\psi _L(z,r)\epsilon rK_1(\epsilon r)(𝐞𝚫)`$ (13)
$`I_{TL}`$ $`=`$ $`{\displaystyle \frac{1}{2}}Qz(1z)(12z)^2K_0(\epsilon r)r\psi _T^{}(z,r)(𝐕^{}𝚫)`$ (14)
In the case of the double-flip $`TT^{}`$ amplitude $`𝐞^2=0`$ and one needs to expand the integrand up to the terms $`(𝐫𝚫)^2`$:

$$\mathrm{exp}\left(\frac{1}{2}i(12z)(𝚫𝐫)\right)\frac{(\mathrm{𝐞𝐫})^2}{r^2}\sigma (x,𝐫,𝚫)\frac{1}{60}(𝐞𝚫)^2r^2[\sigma (\mathrm{},0)+(12z)^2\sigma (r,0)]$$
(15)
so that the corresponding integrand of the double-flip amplitude will be of the form

$`I_{TT^{}}={\displaystyle \frac{1}{15}}z(1z)\epsilon r^2K_1(\epsilon r)\psi _T^{}(z,r)(𝐞𝚫)^2[\sigma (\mathrm{},0)+(12z)^2\sigma (r,0)]`$ (16)
## 3 The scanning radius, hard scale expansion and sensitivity to the short distance wave function of the vector meson

Consider the pQCD regime of large $`Q^2`$. The useful representation for small color dipoles is

$`\sigma (𝐫,0)`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{3}}r^2\alpha _S(q^2)G(x,q^2),q^2{\displaystyle \frac{10}{r^2}}.`$ (17)
Then, because of the exponential decrease, $`K_{0,1}(\epsilon r)\mathrm{exp}(\epsilon r)`$, the amplitudes for the free nucleon target will be dominated the contribution from $`r=r_S`$, where the scanning radius

$$r_S\frac{a_S}{\epsilon }\frac{a_S}{\overline{Q}}=\frac{2a_S}{\sqrt{Q^2+m_V^2}},a_S3.$$
(18)
The simplest case is that of the $`LL`$ amplitude:

$`𝒜_{LL}`$ $``$ $`Qr_S^2\sigma (r_S,0)K_0(a_S){\displaystyle 𝑑zz^2(1z)^2\psi _L(z,r_S)}`$ (19)
$``$ $`{\displaystyle \frac{Q}{\overline{Q}^4}}\alpha _S(\overline{Q}^2)G(x,\overline{Q}^2){\displaystyle \frac{QG(x,\overline{Q}^2)\alpha _S(\overline{Q}^2)}{(Q^2+m_V^2)^2}}.`$
The expansion in powers of the scanning radius $`r_S`$ is an expansion in inverse powers of the hard scale $`\overline{Q}`$. Here one factor $`1/(Q^2+m_V^2)`$ is the same as in the Vector Dominance Model, in the color dipole language it can be identified with the overlap of the photon and vector meson wave functions, the second factor $`1/(Q^2+m_V^2)`$ derives from the pQCD form (17) of the dipole cross section.

The helicity-flip amplitudes will be of the form

$`𝒜_{LT}`$ $``$ $`(𝐞𝚫)r_S^2\sigma (r_S,0)a_SK_1(a_S){\displaystyle 𝑑zz(1z)(12z)^2\psi _L(z,r_S)}`$ (20)
$`𝒜_{TL}`$ $``$ $`(𝐕^{}𝚫)Qr_S^3\sigma (r_S,0)K_0(a_S){\displaystyle 𝑑zz(1z)(12z)^2\psi _T^{}(z,r_S)}`$ (21)
Notice the senisitivity to the short-distance behavior of the vector meson wave wave function in the last result. The soft, oscillator-like interaction would give the wave function which at short distances is a smooth function of $`𝐫^2`$, to that

$$\psi _T^{}(z,r)\frac{r}{R_V^2}\psi _T(z,0)$$
(22)
whereas the attractive Coulomb inetraction at short distances suggests  ”hard”, Coulomb-like $`\psi _T(z,r)\mathrm{exp}(r/R_C)`$ when

$$\psi _T^{}(z,r)\frac{1}{R_C}\psi _T(z,0)$$
(23)
In the case of the soft short-distance wave function the helicity-flip amplitude $`𝒜_{TL}`$ would acquire extra small factor $`r_S/R_V1/(R_V\overline{Q})`$. The similar discussion is relevant to the contribution to the $`𝒜_{TT}`$ from the term $`\epsilon K_1(\epsilon r)\psi _T^{}(z,r)`$ in $`I_{TT}`$, eq. (8), which for the hard short-distance wave function is larger and somewhat enhances the transverse cross section $`\sigma _T`$ and lowers $`\sigma _L/\sigma _T`$, see also the discussion in . Finally, the leading term of expansion in powers of $`r_S`$ of the double-helicity-flip amplitude is of the form

$$𝒜_{TT^{}}(𝐞𝚫)^2z(1z)r_S^3a_SK_1(a_S)\sigma (\mathrm{},0)𝑑zz(1z)\psi _T^{}(z,r_S).$$
(24)
It is proportional to the dipole cross section for the nonperturbative large color dipole

## 4 Nuclear saturation effects in coherent diffraction

In the coherent diffractive production of vector mesons the target nucleus remains in the ground state,

$$\gamma ^{}AVA$$

For heavy nuclei such that their radii are much larger than the dipole size and the diffraction slope in the dipole-nucleon scattering, only the forward dipole-nucleon scattering enters the calculation of the nuclear profile function. Compared to the conventional derivation of the Glauber formulas for the nuclear profile functions, there are little subtleties with the presence of the phase factor $`\mathrm{exp}[\frac{i}{2}(12z)(𝐫𝚫)]`$ in the color dipole factorization formula (4), but a careful rederivation gives the nuclear diffractive amplitudes of a form

$$𝒜_{fi}=2id^2𝐛\mathrm{exp}(i𝚫𝐛)V_f|\left\{1\mathrm{exp}[\frac{1}{2}\sigma (r,0)T(𝐛)]\right\}|\gamma _i^{}\mathrm{exp}[\frac{i}{2}(12z)(𝐫𝚫)],$$
(25)
where $`T(𝐛)=𝑑zn_A(𝐛,z)`$ is the standard nuclear optical thickness at an impact parameter $`𝐛`$. A comparison with the free nucleon amplitude (4) shows that the nuclear profile function

$$\mathrm{\Gamma }(𝐛,𝐫)=1\mathrm{exp}[\frac{1}{2}\sigma (r,0)T(𝐛)]$$
(26)
plays the rôle of the color dipole cross section per unit area in the impact parameter space.

First of all, we notice that the color dipole dependence of the ovelap of wave functions of the photon and vector meson does not change from the free nucleon to the nuclear case. The same is true of the calculation of expectation values over the orientation of color dipoles, see eqs. (14 \- 16). The $`r`$ dependence of the integrands changes, though.

Ref.  gives the detailed discussion of the reinterpretation of the nuclear profile function in terms of the saturating nuclear gluon density  and of the limitations of such an interpretation for observables more complex than the single particle spectra. For the purposes of our discussion it is sufficient to know that in terms of the so-called saturation scale

$$Q_A^2(𝐛)=\frac{4\pi ^2}{3}\alpha _S(Q_A^2)G(Q_A^2)T(𝐛)$$
(27)
the nuclear attenuation factor in (26) can be represented as

$$\mathrm{exp}[\frac{1}{2}\sigma (r,0)T(𝐛)]\mathrm{exp}[\frac{1}{8}Q_A^2r^2]$$
(28)
Let the nucleus be very heavy such that the saturation scale $`Q_A^2`$ is very large (for the estimates of $`Q_A^2`$ for realistic nuclei see ). For color dipoles with $`r^2>r_A^2=8/Q_A^2`$ the nucleus is opaque, i.e., we have a saturated $`\mathrm{\Gamma }(𝐛,𝐫)1`$ independent on the dipole size. The new large scale $`Q_A^2`$ must be compared to $`\overline{Q}^2`$ of eq. (3).

First, there is a trivial case of $`\overline{Q}^2Q_A^2`$. In this case $`r_S^2r_A^2`$, so that $`2\mathrm{\Gamma }(𝐛,𝐫)=\sigma (r_s,0)T(𝐛)`$, i.e., the impulse approximation is at work and the nuclear amplitude has precisely the same structure as the free nucleon one,

$$𝒜_{fi}^{(A)}(𝚫)=𝒜_{fi}^{(N)}(𝚫)d^2𝐛\mathrm{exp}(i𝐛𝚫)T(𝐛)=𝒜_{fi}^{(N)}(𝚫)AG_{em}(𝚫),$$
(29)
where $`G_{em}(𝚫)`$ is the charge form factor of a nucleus.

Much more interesting is the case of coherent diffractive DIS at $`Q^2Q_A^2`$, when the color dipoles in the photon have the dipole size $`r\frac{1}{Q}r_A^2`$, i.e., DIS proceeds in the regime of opacity and saturated color dipole cross section per unit area in the impact parameter space. Apart from this difference all the arguments of section 3 on the scanning radius will be fully applicable, only the scale for the scanning radius will change to $`a_S1`$. However, the, the power expansion in $`r_S`$ will change: instead of the common factor $`r_S^4`$ of section 3 for the free nucleon target one obtains $`r_S^2r_A^2`$ for diffraction off nuclei in the saturation scale, i.e., the substitution

$$\frac{1}{\overline{Q}^4}|_N\frac{1}{\overline{Q}^2Q_A^2}|_A$$
(30)
in the common prefactor of all helicity amplitudes. Otherwise there is no nuclear mass number dependent suppression of the relative strength of the helicity-flip and non-flip amplitudes compared to the free nucleon case. Finally,the $`𝚫`$ dependence of the effective nuclear form factor will be close to the $`𝚫`$ dependence of an amplitude of elastic scattering on a black disc,

$$G_{em}(𝚫)\frac{J_1(R_A\mathrm{\Delta })}{R_A\mathrm{\Delta }}$$
(31)
We emphasize that although the presence of those nuclear form factors limit the practical observation of coherent diffractive DIS to momentum transfers within the nuclear diffraction cone, $`𝚫^2\text{ }<R_A^2`$, and these small $`𝚫`$ cause the kinematical suppression of the helicity-flip within the coherent cone, there is no nuclear suppression of helicity flip even on a black nucleus.

## 5 Incoherent/quasielastic diffraction

In the incoherent (quasielastic, quasifree) diffractive vector meson production

$$\gamma ^{}AVA^{}$$

one sums over all excitations and breakup of the target nucleus without production of secondary particles. The process looks like a production off a quasifree nucleon of the target subject to certain intranuclear distortions of the propagating dipoles. The relevant multichannel formalism has been worked out in , the generalization to the color dipole formalism for $`z\frac{1}{2}`$ is found in . Here we notice that in the color dipole language, the calculation of the helicity amplitudes will be exactly the same as for the free nucleon target but with the extra attenuation factor of eq. (28) in all the integrands. i.e., the inoherent differential cross section equals

$$\frac{d\sigma (\gamma _i^{}AV_fA^{})}{d𝚫^2}=d^2𝐛T(𝐛)\frac{d\sigma _{qel}}{d𝚫^2},$$
(32)
where the amplitudes of the quasielastic production off a quasifree nucleon are given by

$$𝒜_{fi}^{(qel)}(x,𝚫)=i_0^1𝑑zd^2𝐫\sigma (𝐫,𝚫)\mathrm{exp}[\frac{1}{2}\sigma (r,0)T(𝐛)]\mathrm{exp}[\frac{i}{2}(12z)(𝐫𝚫)]I_{fi}(z,𝐫),$$
(33)
In the genuine hard regime of $`Q^2Q_A^2`$ the nuclear attenuation can be neglected and one recovers the free nucleon cross section times the number of nucleons. In the opposite regime of strong saturation, $`\overline{Q}^2Q_A^2`$, the $`r`$ dependence of the attenuation factor is stronger than that of the photon wave functions $`K_{0,1}(\epsilon r)`$. Then, repeating the derivation of the scanning radius in section 3, one will find

$$r_S^2\frac{3}{2}r_A^2.$$
(34)
The functional dependence of helicity amplitudes on the scanning radius $`r_S`$ will be the same as for the free nucleon target with one exception. Namely, the Bessel functions in the photon wave function shall enter with the argument

$$a_S=\epsilon r_S\frac{\sqrt{3}\overline{Q}}{Q_A}1.$$
(35)
In this limit

$$K_0(a_S)\mathrm{log}\left(\frac{Q_A}{\overline{Q}}\right)$$
(36)
which shows that some of the amplitudes will have a logarithmic enhancement, whereas

$$K_1(a_S)\frac{Q_A}{\overline{Q}}$$
(37)
is indicative of even stronger enhancement. However, the closer inspection of the helicity-flip amplitude $`𝒜_{LT}`$ (14) shows that $`K_1(a_S)`$ enters as a product $`a_SK_1(a_S)`$ which is a smooth function at $`a_S1`$. The helicity-flip amplitude $`𝒜_{TL}`$ of eq. (14) exhibits only the weak logarithmic enhancement (36). The case of the helicity-non-flip $`𝒜_{TT}`$ and double-flip $`𝒜_{TT^{}}`$ is a bit more subtle. Here one encounters

$$\epsilon K_1(\epsilon r_S)\psi _T^{}(z,r_S)\frac{1}{r_S}\psi _T^{}(z,r_S)$$
(38)
which for the soft short-distance wave function can be estimated as

$$\frac{1}{R_V^2}\psi _T(z,r_S)$$
(39)
whereas for the hard Coulomb wave function one finds an enhancement

$$\frac{Q_A}{R_C}\psi _T(z,r_S)$$
(40)
To summarize, strong nuclear absorption does not generate any special suppression of the helicity-flip amplitudes compared to the non-flip ones. Furthermore, the estimate (40) suggests even a possibility of an enhancement of the double-flip transitions depending on the hardenss of the short-distance wave function of the vector meson. Finally, this discussion shows that in the saturation regime the saturation scale $`Q_A^2`$ becomes the hard factorization scale for incoherent diffractive production. Namely, for $`\overline{Q^2}Q_A^2`$ this amounts to the substitution

$$\frac{1}{\overline{Q}^4}|_N\frac{1}{Q_A^4}|_A$$
(41)
in the common prefactor of all helicity amplitudes. Similarly, the diffraction slope for the vector meson production will be the same as that for the free nucleon target but taken for the hard scale $`Q_A`$.

Here we focused on the single incoherent scattering approximation. The higher order incoherent interactions can readily be treated following the technique of ref. , they wouldn’t change major conclusions on the interplay of the DIS hard scales $`\overline{Q}^2`$ and the saturation scale $`Q_A^2`$.

## 6 Conclusions

Our principal finding is a lack of nuclear suppression of the helicity-flip phenomena in hard diffractive production off strongly absorbing nuclei, which is in striking contrast to the familiar strong nuclear attenuation of the spin-orbit interaction effects as predicted by the Glauber theory. The QCD mechanism behind this finding is that absorption only affects the color dipole-nucleus scattering amplitude in which the $`s`$-channel helicity of the quark and antiquark is anyway conserved exactly. The helicity flip originates from the relativistic mismatch of the sum of helicities of the quark and antiquark and the helicity of the vector meson and photon. Within the color dipole approach we demonstrated how the expansion of helicity amplitudes in powers of the scanning radius and saturation radius changes from the free nucleon to coherent nuclear to incoherent (quasielastic) nuclear diffractive production.

It is proper to recall the early claims by Greenberg and Miller of the so-called vector color transparency - the vanishing spin-flip phenomena in hard processes on nuclear targets . These authors considered the polarization of ejected protons in quasielastic scattering $`A(e,e^{}p)`$ off nuclei in the version of spin-orbit interaction model and overlooked Zakharov’s helicity mismatch mechanism which, as we demonstrated above, persists in hard scattering. We discussed diffractive DIS, but all the arguments are equally applicable to the hadron-nucleus scattering.

The coherent and incoherent diffractive vector meson production off nuclei can be studied experimentally in the COMPASS experiment at CERN . Whereas we are confident in our predictions for perturbatively large saturation scale $`Q_A^2`$, the numerical estimates for the saturation scale give disappointingly moderate $`Q_A^21`$ GeV<sup>2</sup>. None the less, the qualitative pattern of predicted changes of the $`Q^2`$ dependence of diffractive vector meson production from the free nucleon to coherent nuclear and incoherent nuclear cases must persist even at moderately large $`Q_A^2`$.

This work was partly supported by the grant INTAS 00366.


text:
# A quasi-local measure of inter-scale transfer: An approach to understanding turbulence

## I Introduction

Turbulence is characterized by the transfer of an inviscid constant, such as kinetic energy or a passive scalar, between different length scales. One of the primary goals of turbulence research is to understand the mechanisms that drive this transfer process. Numerical simulations performed by Farge et al. suggest that in both two-dimensional (2D) and three-dimensional (3D) turbulence it is a relatively small number of “coherent structures” that dominate the turbulent dynamics Farge et al. (1999, 2001). As a result, it is important to understand how the existence and interaction of these coherent structures affect the transfer processes. Thus, the spatially-local transfer properties of the flow need to be measured and correlated with these structures.

In this manuscript a tool for obtaining local information about the scale-to-scale transfer of inviscid constants from experimental/numerical data is examined. The method, called the “filter approach” (FA), has traditionally been applied to large-eddy numerical simulations (LES) Lesieur (1987), but is developed here in the context of experimental data analysis. By applying a low-pass spatial filter to the equations of motion (the incompressible Navier-Stokes equation), separate equations for the filtered, or large-scale, fields and the remaining small-scale fields can be derived Germano (1992). Within the resulting equations are coupling terms that represent the interaction of the large- and small-scale fields with each other. In LES schemes the coupling terms in the large-scale equations are modeled thereby eliminating the necessity of directly computing the small scales Lesieur (1987). In our application of the filter approach, however, data from direct numerical simulations (DNS) or high resolution experiments are used to directly evaluate the coupling terms and obtain a quasi-local measure of the inter-scale transfer Tao et al. (2002); Rivera et al. (2003).

For readers familiar with LES and turbulence modeling it should be stressed that the philosophy driving our use of the filter approach is different from the a priori development of LES models discussed in Menevau and Katz (2000). There the objective was to determine empirically which of several LES modeling schemes for the large- to small-scale coupling terms most successfully emulates physical data. In that case, FA was used primarily as a benchmark, and only went as far as measuring the inter-scale coupling term. Rather than producing results for numerical benchmarking, FA can also be used as an analysis probe to determine where and when in a flow scale-to-scale transfer of inviscid constants takes place. In this way one can isolate important interaction events and form an understanding of turbulence inter-scale transfer mechanisms. Of course, this understanding could eventually be incorporated into LES models.

The filter approach is applied in this paper to both experimental and numerical data. The intent is not to investigate the underlying physics of the turbulence, which will be presented in later papers, but rather to determine the appropriate interpretation of the results and the limitations imposed by different filters, spatial boundaries, and finite measurement resolution.

## II Experimental and Numerical Systems

Experimental measurements were carried out in a flowing soap-film channel, a quasi-2D system in which decaying turbulence of low to moderate Reynolds number can be generated ($`10^2Re10^4`$). The channel was $`5`$ cm wide and was inclined at an angle of $`75^{}`$ with respect to vertical. The mean flow was $`120`$ cm$`/`$s and the film thickness was about $`10`$ $`\mu `$m. A more detailed description of the channel can be found in Rivera et al. (1998); Vorobieff et al. (1999). Using the empirical relationships measured in Vorobieff and Ecke (1999), the films’ kinematic viscosity was $`\nu 0.03`$ cm2/s. The turbulence generating grid consisted of rods of $`0.12`$ cm diameter with $`0.22`$ cm spacing between the rods. Thus, the blocking fraction is around $`0.3`$, which is typical for turbulence in 2D soap film flows Gharib and Derango (1989); Kellay et al. (1995); Rivera et al. (1998); Rutgers et al. (2001). The resulting Reynolds number, $`Re=UL/\nu `$, was $`880`$ based on the mean-flow velocity and an injection scale of $`L_{\text{inj}}=0.22`$ cm. The turbulent velocity, $`𝐮(𝐱)`$, and vorticity, $`\omega (𝐱)`$, fields generated by the grid were obtained by tracking $`35`$ $`\mu `$m polystyrene spheres (density approximately $`1.05`$ g/cc) within a $`1.8\times 1.8`$ cm2 region located $`6`$ cm downstream from the grid (20-30 eddy rotation times)Ishikawa et al. (2000); Ohmi and Li (2000). The particles were illuminated with a double pulsed Nd:Yag laser and their images captured by a 12-bit, $`2048\times 2048`$ pixel camera. Around $`3\times 10^4`$ particles were individually tracked for each image pair and their velocities and local shears were interpolated to a discrete $`135\times 135`$ grid. One-thousand velocity and vorticity fields were obtained in this way and were used to compute ensemble averages of the statistical measures described below. Typical velocity and vorticity fields are shown in Fig. 1.

To supplement the experimental data, a direct numerical simulation of the 2D Navier Stokes equation was performed. Computational details are presented in CHen et al. (2003). The equation

$$\frac{\omega }{t}+u_j\frac{\omega }{x_j}+\nu _m(1)^m\frac{^{2m}\omega }{x_j^{2m}}=F$$
(1)
was simulated in a square domain with side $`L=2\pi `$ and periodic boundary conditions. Here, $`u_i`$ is the $`i`$-th component of the velocity; $`\omega =ϵ_{ij}_ju_i`$ is the vorticity; and $`F`$ is a stirring force applied to wave numbers $`|𝐤|=47`$. The Einstein summation convention is used throughout. Two values of $`m`$ are considered: $`m=1`$ corresponding to Laplacian viscosity ($`\nu _1=0.01`$), and $`m=8`$ corresponding to hyper-viscosity ($`\nu _8=1.4\times 10^7`$), which has the effect of extending the inertial range.

The equation was solved using a fully de-aliased, parallel pseudo-spectral code with second-order Adam-Bashforth time-stepping. The resolution was $`2048^2`$. A statistically-stationary state was achieved after about 200 large-eddy turn-over times. Representative examples of vorticity fields generated using Laplacian and hyper-viscosity are shown in Fig. 2.

## III The Filter Approach

### III.1 General Application

The filter approach allows the direct measurement of the coupling between scales in non-linear systems. It was originally developed as a tool to truncate numerical simulations of turbulence by modeling small-scale behavior with the knowledge of large-scale behavior, i.e. LES Lesieur (1987). In this paper, the general application of FA to nonlinear systems will be considered first, followed by a specific example: the case of energy and enstrophy transfer between length scales in 2D turbulence.

General features of the transfer process are examined by considering a field, $`q(𝐱)`$, that is evolved by a non-linear evolution operator, $`F^{(\text{nl})}`$, as $`F^{(\text{nl})}q=0`$. Given a scale, $`l`$, one can separate the field, $`q`$, into large-scale, $`q_l`$, and small-scale, $`q_s`$, components. This is done by convolving a filter function, $`G_l`$, with $`q`$ and defining $`q_lG_lq`$ and $`q_sqq_l=(\delta G_l)q`$, where $`\delta `$ is the Dirac delta function. Applying $`G_l`$ to the evolution equation for $`q`$ yields $`G_l(F^{(\text{nl})}q)=0`$. Adding zero, written as $`F^{(\text{nl})}q_lF^{(\text{nl})}q_l`$, results in:

$$F^{(\text{nl})}q_l=(G_l(F^{(\text{nl})}q)F^{(\text{nl})}q_l).$$
(2)
Eq. (2) indicates that the large-scale field, $`q_l`$, evolves in exactly the same manner as the full field, $`q`$, up to a coupling term $`C(q_l,q)(G_l(F^{(\text{nl})}q)F^{(\text{nl})}q_l)`$. This coupling term, which can be considered as an external forcing or damping, arises from the small-scale field, $`q_s`$, interacting with $`q_l`$. Thus, $`C`$ provides a measure of the interaction between large and small scales. The measurement of $`C`$ requires no assumptions about the field, $`q`$, such as homogeneity or isotropy. Moreover, $`C`$ is a field quantity, not simply an average, and can reveal not only information about the coupling between scales, but also when and where such interactions are taking place and with what strength.

The field nature of the coupling term derived above, $`C`$, makes FA very valuable in the context of studying turbulence. As mentioned in the introduction, turbulence is characterized by an average scale-to-scale flux of quantities such as the kinetic energy. Many possible mechanisms underlying the transfer processes have been suggested, e.g., the stretching of vortex tubes into thin filaments may account for some fraction of the down-scale transfer of energy in 3D Frisch (1995). In 2D, vortex merger has been postulated as a way of transferring energy to larger scales Lesieur (1987). The difficulty is that neither of these pictures has been conclusively correlated with topological flow structures or the underlying transfer dynamics, though some attempts have been made Daniel and Rutgers (2002). By using FA to measure the inter-scale transfer and correlating this field with the position of flow structures (identified by other means), one can determine the veracity of these hypothesized transfer mechanisms.

There are a number of subtleties to consider when applying FA to experimental data. First, there is the choice of the filter function, $`G_l`$. The interpretation of the scale-to-scale transfer depends on the selection of the filter function. For example, if the convolution function is Gaussian defined by length scale $`l`$, the interpretation of the filtered functions is as given above. On the other hand, if the convolution is with the kernel $`H_l(\delta G_l)`$, where again $`G_l`$ is Gaussian and $`\delta `$ is a Dirac delta function, the resulting interpretation of the convolved fields, $`H_lf`$, as “large-scale” is incorrect. Indeed, the function $`H_l`$ produces “small-scale” fields. Other ramifications of changing the filter function will be explored shortly but the convention is adopted that the filter function is always low-pass.

Two additional considerations when applying FA to real data are finite measurement resolution and the existence of spatial boundaries. The extent to which the measurements are sensitive to either of these factors depends on the quality of the data and on the form of the filter used. Significant research in the context of LES has addressed similar concerns, but always for the purpose of approximating the physical system in a numerical simulation Lesieur (1987). These issues are explored in some depth in a later section.

### III.2 Application to two-dimensional turbulence: energy and enstrophy transport

In two-dimensional turbulence kinetic energy, $`Eu^2/2`$, and enstrophy, $`\mathrm{\Omega }\omega ^2/2`$, are conserved in the inviscid limit. Theory, numerics, and experiments all indicate that energy is transferred on average from small to large length scales (up-scale) and that enstrophy is transferred in the opposite direction (down-scale) Kraichnan and Montgomery (1980). FA, however, can be used to obtain much more detailed information about the energy and enstrophy transfer within the flow. We begin by examining the scale-to-scale coupling of energy in the 2D Euler equation (viscous terms are linear and, hence, cause no direct scale-to-scale transfer, eliminating the need to examine the full Navier-Stokes equation):

$$\frac{u_i}{t}+u_j\frac{u_i}{x_j}=\frac{p}{x_i},$$
(3)
where $`u_i`$ is the $`i`$-th component of the velocity field, $`p`$ is the density normalized pressure field, and summation over repeated indices is assumed. Contracting this evolution equation with a filter function, $`G_l`$, and extracting the coupling term as described above yields:

$`{\displaystyle \frac{(u_i)_l}{t}}`$ $`+`$ $`(u_j)_l{\displaystyle \frac{(u_i)_l}{x_j}}`$ (4)
$`=`$ $`{\displaystyle \frac{p_l}{x_i}}{\displaystyle \frac{}{x_j}}((u_iu_j)_l(u_i)_l(u_j)_l).`$
The notation $`(f)_l`$, or simply $`f_l`$, will be used to denote the large-scale field $`G_lf`$. Eq. (4) is almost the equivalent of Eq. (2). The one delicacy is that the term $`p_l`$ is not the large-scale pressure field, i.e., it is not the field obtained using the gradients of $`(u_i)_l`$. This is not an important issue here, however, because, as will be demonstrated later, $`p_l`$ does not contribute to the inter-scale transfer.

From Eq. (4) the scale-to-scale coupling term for the velocity is given by

$$C((u_i)_l,u_i)=\frac{\tau _{ij}^{(l)}}{x_j},$$
(5)
where $`\tau _{ij}^{(l)}=(u_iu_j)_l(u_i)_l(u_j)_l`$ is the subgrid-scale stress tensor. To obtain the equation for the evolution of large-scale energy, $`E^{(l)}=(u_i)_l^2/2`$ one multiplies Eq. (4) by $`(u_i)_l`$. (Note that the notation $`E^{(l)}`$ is used rather than $`E_l`$ since $`E_l=(u_i^2)_l/2(u_i)_l^2/2=E^{(l)}`$.) The resulting equation,

$`{\displaystyle \frac{E^{(l)}}{t}}+{\displaystyle \frac{}{x_j}}((u_j)_lE^{(l)}+(u_j)_lp_l)=(u_i)_l{\displaystyle \frac{\tau _{ij}^{(l)}}{x_j}}.`$ (6)
for the energy contained at scales larger than $`l`$, is identical to the full energy evolution equation up to the coupling term on the right hand side (again ignoring the pressure term).

This coupling term is not yet in the form to directly yield scale-to-scale energy transfer information. There are two ways in which the small-scale velocities can change the large-scale energy: by physically transporting it from point-to-point or by transferring it between scales. To separate the two, the Leibnitz rule is used to rewrite the right hand side of Eq. (6) as

$$=\frac{}{x_j}((u_i)_l\tau _{ij}^{(l)})+\tau _{ij}^{(l)}\frac{(u_i)_l}{x_j}.$$
(7)
Notice that the latter of these two terms is Galilean invariant, whereas the former is not. Boosts to the reference frame should not change the scale-to-scale transfer of energy but will change the point-to-point transport. Therefore, the former term is attributed to the point-to-point transport and the latter to the scale-to-scale transfer.

Another way to contrast the point-to-point coupling with scale-to-scale coupling is to consider the limit of a homogenous system. In this case taking an ensemble average, $`\mathrm{}`$, should eliminate all point-to-point transport terms, leaving only inter-scale transfer contributions. Since the action of ensemble averaging commutes with the derivative operation, and since the spatial derivative of an ensemble average is zero in the limit of homogeneity, the ensemble average of Eq. (6) is simply $`_tE^{(l)}=\tau _{ij}^{(l)}_j(u_i)_l.`$ Not only does this demonstrate that $`_j[\tau _{ij}^{(l)}(u_i)_l]`$ is a point-to-point term, but it also demonstrates that all of the terms on the left, other than the time derivative, are point-to-point as well. This fact allows us to ignore the delicacy with respect to the pressure term: it does not affect scale-to-scale transfer.

For simplicity Eq. 6 is rewritten as

$$\frac{E^{(l)}}{t}+\frac{J_j^{(l)}}{x_j}=\mathrm{\Pi }^{(l)},$$
(8)
where $`J_j^{(l)}`$ and $`\mathrm{\Pi }^{(l)}`$ are defined as

$`J_j^{(l)}`$ $``$ $`(u_j)_lE^{(l)}+(u_i)_l\tau _{ij}^{(l)}+(u_j)_lp_l,`$ (9)
$`\mathrm{\Pi }^{(l)}`$ $``$ $`\tau _{ij}^{(l)}{\displaystyle \frac{(u_i)_l}{x_j}}=\tau _{ij}^{(l)}(S_{ij})_l.`$ (10)
where $`(S_{ij})_l`$ is the large-scale strain tensor. The negative sign in the definition of $`\mathrm{\Pi }^{(l)}`$ is added so that down-scale transfer has a positive value whereas up-scale transfer is negative.

An almost identical method can be used to determine the scale-to-scale enstrophy transfer of the flow. The starting point, however, is the 2D Euler equation for vorticity,

$$\frac{\omega }{t}+u_j\frac{\omega }{x_j}=0.$$
(11)
As above, the equation is contracted with a filter function, $`G_l`$, and the coupling term is extracted,

$`{\displaystyle \frac{\omega _l}{t}}+(u_j)_l{\displaystyle \frac{\omega _l}{x_j}}={\displaystyle \frac{\sigma _j^{(l)}}{x_j}},`$ (12)
where $`\sigma _j^{(l)}=(u_j\omega )_l(u_j)_l\omega _l`$ is the subgrid scale vorticity transport vector. Notice that the term $`\sigma _j^{(l)}`$ defined in Eq. 12 has an almost identical form to $`\tau _{ij}^{(l)}`$ in the energy equations. This general form is typical of the filter approach. The coupling terms between large- and small-scale fields often take the form $`(ab)_la_lb_l`$ for quadratic nonlinearities.

To change the large-scale vorticity equation to an equation for large-scale enstrophy evolution, one must multiply by $`\omega _l`$. This yields

$`{\displaystyle \frac{\mathrm{\Omega }^{(l)}}{t}}`$ $`+`$ $`{\displaystyle \frac{}{x_j}}((u_j)_l\mathrm{\Omega }^{(l)})=\omega _l{\displaystyle \frac{\sigma _j^{(l)}}{x_j}}`$ (13)
$`=`$ $`{\displaystyle \frac{}{x_j}}(\omega _l\sigma _j^{(l)})+\sigma _j^{(l)}{\displaystyle \frac{\omega _l}{x_j}},`$
where again the Leibniz rule was used to separate point-to-point transport from inter-scale transfer. Grouping the appropriate terms, as was done for the energy equation, leads to the final form,

$$\frac{\mathrm{\Omega }^{(l)}}{t}+\frac{K_j^{(l)}}{x_j}=Z^{(l)},$$
(14)
where

$`K_j^{(l)}`$ $``$ $`(u_j)_l\mathrm{\Omega }^{(l)}+\sigma _j^{(l)}\omega _l`$ (15)
$`Z^{(l)}`$ $``$ $`\sigma _j^{(l)}{\displaystyle \frac{\omega _l}{x_j}}`$ (16)
where the negative sign in front of the scale-to-scale coupling term, $`Z^{(l)}`$, again makes down-scale transfer positive.

To make the above derivations more concrete, some examples are provided of filtered fields. For purposes of illustration, vorticity and enstrophy fields are presented and an analysis of velocity and energy fields is left to a later publication. Figure 3 displays steps in the calculation of the scale-to-scale enstrophy transfer, $`Z^{(l)}`$, for a typical vorticity field extracted from experimental data. For this calculation the filter function was a Gaussian with Fourier-space definition

$$G_l(𝐤)=𝐞^{\frac{|𝐤|^\mathrm{𝟐}}{𝐤_𝐥^\mathrm{𝟐}}},$$
(17)
where $`k_l2\pi /l`$. These figures illustrate that the application of FA is straightforward: (1) compute the secondary field, $`u_i\omega `$, from the measured velocity and vorticity fields; (2) perform a convolution of these fields to obtain $`\sigma _i^{(l)}`$; (3) take the scalar product of $`\sigma _i^{(l)}`$ with the appropriate gradient of the large-scale fields, namely $`_i\omega _l`$.

There is, however, a caveat for the general case. It may not always be possible to separate the point-to-point transport from the scale-to-scale transfer terms. For the energy equation, the terms were determined by using the Leibniz rule to separate out the Galilean invariant part of the energy flux. For the enstrophy, the separation was obtained by analogy with the energy equation rather than by a strict application of Galilean invariance. There is no a priori expectation that such a separation will be as simple, or even possible, for arbitrary nonlinear systems. In compressible flows, for example, it is possible to measure the coupling terms, but the point-to-point transport caused by small scales is tied to the scale-to-scale transfer in a non-trivial way. The interpretation of such results must, therefore, be done carefully.

## IV Systematic effects

This section will address the effect of varying the filter, $`G_l`$, on the interpretation of the results, the consequences of the data being limited in spatial extent (i.e. by boundaries), and variations in the results caused by experimental or numerical limitations on the resolution of the data.

### IV.1 Different Filters

The interpretation of the results of FA may depend on the form of the filter function, $`G_l`$. A Gaussian kernel was used in Eq. (17), but this particular choice was made only because it has a simple interpretation in both real- and Fourier-space, taking the same form in both. The FA technique imposes no such constraints in general. In some cases, one may want to preferentially constrain the filter in real-space to a well-defined length scale or, instead, may want to filter so as to select only a sharp band of wave numbers. Either possibility can be explored using this tool. It must be kept in mind, however, that sharpening the filter in real- or Fourier-space causes a corresponding broadening of the filter in the other. The impact and interpretation of varying the form of the filter on the resulting fields is considered here.

The form of the Fourier filters used is

$$G_l^{(n)}(𝐤)=𝐞^{\left(\frac{|𝐤|}{𝐤_𝐥}\right)^n},$$
(18)
where $`n`$ is the filter order. The case $`n=2`$ is the Gaussian filter considered earlier. As $`n`$ increases, the filter sharpens around the Fourier mode corresponding to filter length $`l`$. Similarly, the real-space filters are defined by:

$$G_l^{(n)}(𝐱)=\mathrm{𝐀𝐞}^{\left(\frac{|𝐱|}{𝐥/\pi }\right)^n}.$$
(19)
Note that for $`n=2`$, the real filter is equivalent to the Fourier filter. As the real filter is sharpened it approaches an area average over a box of diameter $`2l/\pi `$. Both of these sets of filters are shown in Fig. 4. Other types of filter are possible, but are not considered here.

From an experimental point of view, the sharper real-space filters are more attractive than the Fourier filters because the real-space envelope of the Fourier filters grows as the order of the filter is increased. Since experimental data is invariably windowed to the cross-section of the measurement apparatus, this means that sharper Fourier filters quickly grow to interact with boundaries. For sharp real filters this is not a problem: as the filter increases in sharpness it becomes more spatially compact. For the purpose of comparing with theory, however, sharper Fourier filters approach the ideal; see, for example, the discussion given in Frisch Frisch (1995) in which the fields are filtered by an infinitely sharp cutoff in Fourier-space.

In Fig. 5, the results of using a Gaussian, two low-order Fourier-space filters, and two low-order real-space filters in the calculation of the enstrophy transfer are shown (the associated vorticity field is also displayed). Superficially, the fields are fairly similar. The strength of the fluctuations, however, changes with filter (the grey scale limits of the fields are noted in the captions). The change in magnitude is stronger for Fourier filters than real-space filters, with the $`n=4`$ Fourier-space filter experiencing $`50\%`$ larger fluctuations than for $`n=2`$.

The similarity in the fields seems to indicate that the relative magnitudes of enstrophy flux remain more or less constant. In particular, the stronger values of flux associated with powerful vortices (middle top and middle left) have almost identical forms, though there may be a slight sharpening of the lobes for both higher order real-space and Fourier-space filters (which, we note, have a quadrupolar form). Higher-order Fourier-space filters increase the symmetry of many features, i.e. they appear less ovular (see for example the lobes indicated by the white arrows in Fig. 5), whereas the contrast between features and the background is enhanced by higher-order real-space filters. On the whole, however, the qualitative features are fairly insensitive to the form of filter used, maybe surprisingly so. One might expect that the beating of the sharper Fourier kernel in real space would show up more strongly in the fields. For the stronger enstrophy transfer signals this does not appear to be the case.

Weaker signals are more sensitive to the particular choice of filter. Consider the weak lobe in the lower middle of the field indicated by the black arrow in Fig. 5 (or the one in the middle right also indicated by a black arrow). For the Gaussian filter, the lobe is barely visible. For sharper filters in both real- and Fourier-space ($`n=4`$), however, the prominence of this lobe with respect to the stronger signals in the system is enhanced.

The most important observation is that the qualitative structure of the field is relatively unaffected by the choice of filter. Sharper filters increase the “contrast”, but neither eliminate nor create structural features in the flow. For the purposes of correlating scale-to-scale transfer to topological features, this is an important feature of FA.

The above comparison is only qualitative. A more accurate quantitative comparison is presented in Fig. 6, where the probability distribution function of the enstrophy transfer is presented for the same set of filters used in Fig. 5. The agreement between the different filters is quite good, but a little deceptive. Note, first, that the magnitude of the RMS fluctuations has been normalized out. Second, there may be (though it is, perhaps, below the noise level) a slight increase in asymmetry, in particular in the negative tail of the PDF for real-space filters (the open symbols). Given that in Fig. 5 the large values of enstrophy transfer, corresponding to the tails of the PDF, are relatively insensitive to changes in the filter, this collapse of the PDFs is reasonable.

On the other hand, in Fig. 5 the weaker values of inter-scale transfer have somewhat increased contrast. This change is emphasized in the lowest-order moments of the distribution, rather than in the tails. The average, $`Z^{(l)}`$, and fractional sign-probability comparison, $`\delta P(Z^{(l)})P(Z^{(l)}>0)P(Z^{(l)}<0)`$, are shown in Fig. 7. Here, there is a significant difference in the average enstrophy flux as a function of filter order. In particular, the average rises and falls more sharply for the higher-order Fourier-space filters than it does for the real-space filters. Also, the area fraction saturates at a smaller value, then falls more quickly. The signs of the average flux and the area difference, however, are quite robust, although the magnitudes seem to vary (in some places by a factor of two over this range of filter orders). The rise in the peak of average scale-to-scale transfer is reminiscent of ringing such as takes place in the Gibbs phenomenon. It is not clear whether or not this is the source of the change.

#### IV.1.1 Interpretation

At this point, one might ask: for which filter is the result closest to the “real” enstrophy flux? This question depends entirely on what one means by “enstrophy flux”. The scale-to-scale transfer between wave numbers is most closely approximated by higher-order Fourier filters. Because of the associated broadening of the filter in real-space, however, the resulting fields are not as good a measure of the spatially-local enstrophy transfer. On the other hand, for the flux to be localized in physical space for comparison with real-space structures, sharp real-space filters are preferable. In this limit, the enstrophy flux can no longer be defined as the movement of enstrophy from small Fourier modes to larger ones but is actually a measure of the flux out of some bands and into others (which do not necessarily have larger wave number). The strength of FA does not lie in its ability to measure the exact magnitude of the transfer, but rather that the sign of the transport and the qualitative form of the fields is robust to changes in the filter.

### IV.2 Boundaries

An important consideration when using FA, one that is also a major issue in LES, is the presence of boundaries (either physical or as limits of the viewing window). As a Fourier-space filter grows in order, and correspondingly in real-space extent, less and less of the data can be used as the boundaries begin to affect the computation of the convolution in regions farther and farther into the interior.

To investigate the influence of finite system size a random periodic streamfunction was generated on a $`512\times 512`$ grid. This streamfunction was then used to obtain velocity and vorticity fields from which the enstrophy flux field, $`Z_0^{(l)}`$, was computed. The left half of the fields was then set to zero and the enstrophy field, $`Z_b^{(l)}`$, recalculated. The normalized RMS difference,

$$Z_{\text{error}}^{(l)}=\frac{(Z_b^{(l)}Z_0^{(l)})^2^{1/2}}{(Z_0^{(l)})_{\text{RMS}}},$$
(20)
is shown in Fig. 8 as a function of $`x/l`$, the distance from the introduced boundary normalized by the filter size. The various plots are given for Fourier-space filters of different orders (see Eq. 18).

With the exception of the order $`3`$ filter there is a continuous increase in the boundary affects. Taking a nominal error rate of $`10^4`$ as acceptable, a boundary of $`bl(n1)`$ is appropriately sized for the data. This results in a loss of $`2l(n1)`$ in linear box size since one must apply the condition to left-right (top-bottom) boundaries. This effective boundary has been adhered to throughout this paper. For real-space filters the boundaries become less of an issue. Indeed, for the sharpest real-space filter (a simple average over a circle), the boundary is $`l/2`$.

### IV.3 Effects of Finite Resolution

In this section the effect of finite measurement resolution on the ability of FA to resolve the behavior of the enstrophy flux for a given filter function is explored. In the preceding sections experimental data was used in the analysis, but to estimate the effects of finite resolution it is necessary to consider numerical data where the range of scales and the uncertainty in the measured values is fairly well known. The numerical data also eliminates any concerns about the effect of boundaries since the boundary conditions are periodic.

Pre-filtering at the limit of the data resolution should result in no apparent difference in the transfer from that computed from the original field. As the pre-filter length, $`L`$, is increased above the resolution scale differences in the measured transfer will increase. The question this raises is: how far above $`L`$ are these finite-resolution effects felt? To answer this question a series of $`20^{\mathrm{th}}`$ order real-space filters with varying $`L`$ are applied to two data sets simulating the local averaging inherent in experimental measurements (such as image acquisition and particle tracking). From these pre-filtered fields the enstrophy transfer, $`Z_L^{(l)}`$, is calculated and compared with the full flux.

The results of these calculations are shown for both data sets in Fig. 9. For both, the viscous cut-off scale is between $`0.01`$ and $`0.02`$. For pre-filters with $`L`$ at or below this length there is little change in the magnitude of the average flux (slightly more prominent in the case of hyper-viscosity). Once the filter size exceeds the viscous scale, i.e., penetrates at all into the inertial range, the effects are quickly felt up to almost the injection scale, $`0.4`$. This may be a general feature of FA, or it may be a peculiarity of the 2D enstrophy transfer process being non-local. In either case, it demonstrates the importance of resolving the entire inertial range, including the viscous scale, for an accurate representation of the enstrophy transfer to be possible. Whether or not enstrophy transfer is local in Fourier-space is an important question in its own right, and one which is sufficiently complex to warrant exploration in a separate paper.

## V Conclusions

In the above discussion and use of FA, familiar turbulence assumptions, such as homogeneity, isotropy, and inertial range, were not requirements of the measurement. This is perhaps the greatest strength of FA: it can be applied to any type of turbulent flow, regardless of that particular flow’s properties. All that is needed is the evolution equation for the system from which the inter-scale transfer can be obtained. The flow doesn’t have to be turbulent, FA can be applied to laminar or periodic flow and reasonable results obtained. In other words, the interpretation of FA does not rely on a pre-existing theory. Rather, it is a tool that can directly test notions of how inter-scale transfer takes place in systems and can be used to build appropriate theories.

This “bottom up” approach to physics (rather than a theoretical trickle down) comes at a price: one must have highly resolved fields of data with a significant range of spatial scales. The former is a requirement imposed by the resolution issues discussed above, which can only be relaxed if one assumes a wavenumber local inter-scale transfer process. The latter is a limitation imposed by the data windowing and interaction of the filters with the measurement boundary. Of course, the standard techniques of measuring velocity fields in fluids, namely particle imaging velocimetry or particle tracking, more or less ensure that the measurements are not far from the viscous scale, and thus of high enough resolution to not assume local transfer. This is because these techniques rely on groups of particles moving coherently, which only holds when the scales being probed are small enough that local Taylor expansions describe the flow. The difficulty is in simultaneously measuring a significant range of scales above the viscous scale for meaningful information to be obtained. With standard PIV practices this is possible in 2D. Only with holographic PIV is this possible in 3D. There, however, the limitation is in the amount of data that can be obtained (order 10 fields is possible; currently, 1000 fields are not).

As discussed in the introduction, our ultimate intention for the FA technique is to probe the mechanisms driving the inter-scale transfer in 2D turbulence, with a particular eye towards the role of coherent structures in the transfer process. This will be done both via statistical analysis of fields in the laboratory frame of reference (Eularian frame), as well as following the motion of fluid parcels (Lagrangian frame). FA, however, is not limited to these measurements alone, but could find application to the inter-scale transfer and mixing of passive scalars, or other quantities not necessarily conserved in the inviscid limit (i.e. local topology). And as hinted at previously, FA can be applied (carefully!!!) to general non-linear systems where an equation of motion is known and inter-scale transfer is of interest.

## VI acknowledgments

We thank Greg Eyink, Phil Marcus, Misha Chertkov and Boris Shraiman for interesting discussions and suggestions.