text:
# Deep VLT search for globular clusters in NGC 5128: color-magnitude diagrams and globular cluster luminosity function Based on observations collected at the European Southern Observatory, Paranal, Chile, within the Observing Programme 63.N-0229, and at La Silla Observatory, Chile.

## 1 Introduction

Globular Clusters are among the oldest objects in the Universe. Moreover, they are intrinsically luminous and found around almost all galaxies in which searches have been made. There are growing evidences for globular cluster formation during major star formation episodes (Schweizer schweizer (1997); Larsen & Richtler lr99 (1999)). This makes them useful tracers of the formation and the evolution of galaxies (Harris harris91 (1991, 2000); Kissler-Patig et al. kissler3 (1998); van den Bergh vandenbergh (2000)).

The globular cluster luminosity function (GCLF) of a galaxy is the relative number of its globular clusters per luminosity (or magnitude) interval. It is usually unimodal and nearly symmetric with a peak at a characteristic (turn–over) magnitude. Empirically it was found that turn–over magnitudes for galaxies at the same distance (e.g. in the same galaxy cluster) differ by $`<0.1`$, which makes the GCLF a good standard candle (Harris harris91 (1991); Jacoby et al. jacoby92 (1992)). However, the full calibration of the GCLF is incomplete. The essential problem is that the galaxies at large distances, for which the GCLF method is mostly used, are giant ellipticals in which globular clusters are found in largest numbers and in which there are usually fewer problems with extinction, whilst the primary calibration of the standard candle $`M_V^0`$ rests on the MW and M31, which are spiral galaxies. In fact, evidence that GCLFs depend on Hubble type and environment has been claimed (Fleming et al. fleming (1995)), but the former is probably due to differences in mean metallicity of the globular clusters and the latter was never reliably demonstrated (Kissler-Patig kissler00 (2000)).

There is strong observational evidence that most luminous giant ellipticals have a bimodal globular cluster metallicity distribution (Gebhardt & Kissler-Patig gebh99 (1999); Kundu kundu (1999)). The globular cluster specific frequencies (the number of clusters normalized to the galaxy luminosity; Harris & van den Bergh hvdenb81 (1981)) are typically higher in ellipticals than in spirals (Harris harris91 (1991, 2000)). To explain the high number of globular clusters and the presence of several populations indicated by the bimodal metallicity distributions, different globular cluster formation mechanisms have been proposed (e.g. Ashman & Zepf ashmanzepf92 (1992); Côté et al. cmw98 (1998); Forbes et al. forbes97 (1997)). The number, chemical abundance and spatial distribution of a globular cluster system (GCS) put strong constraints on the formation of its parent galaxy.

At the distance of 3.6 Mpc (Soria et al. soria (1996); Harris et al. hetal99 (1999)) NGC 5128 is the closest giant elliptical galaxy. The first globular cluster in NGC 5128 was discovered by Graham & Phillips (GP80 (1980)). Its brightness was high with $`V=17.2\pm 0.1`$ mag, which corresponds to $`M_V=10.6`$ mag. Subsequent investigation of photographic plates (van den Bergh et al. vBHH (1981), Hesser et al. HHvBH (1984)hetal86 (1986)) revealed a rich GCS in this galaxy. The estimated number of globular clusters in NGC 5128 is N$``$1700 (Harris et al. hetal84 (1984)) compared with N$``$150 for the MW.

On the basis of the near IR colors for a magnitude–limited sample of 12 globular clusters, Frogel (Frogel80 (1980, 1984)) found one third of the sample to be super metal rich (above solar metallicity) and luminous. However, the Washington photometry of 62 spectroscopically confirmed clusters (Harris et al.  hetal92 (1992)) as well as a spectrophotometry for 5 clusters (Jablonka et al.  jablonka (1996)) showed less extreme metallicity. Harris et al.  (hetal92 (1992)) found a large spread in metallicities with a mean of \[Fe/H\]$`{}_{CT_1}{}^{}=0.8\pm 0.2`$ dex. They also find evidence for a weak spatial metallicity gradient, but state that the same could be produced by unknown reddening effects. Zepf & Ashman (zepfashman93 (1993)), using the same photometric data, suggest that the metallicity distribution of globular clusters in NGC 5128 is bimodal, with the higher metallicity peak at \[Fe/H\]=+0.25 dex, and proposed that this was due to the formation of a second population of globular clusters in a previous merger which created NGC 5128.

Several photometric studies in the inner regions of NGC 5128 have been conducted, with the aim of constraining the dependence of globular cluster metallicity on galactocentric radius and of finding super-metal-rich clusters (Minniti et al. minniti96 (1996); Holland et al. holland (1999)). These studies suggest the existence of radial gradients and bimodal distribution in globular cluster metallicity. However, their conclusions are limited due to the restricted samples, possible reddening within NGC 5128 that is especially severe in the central parts around the prominent dust lane, and possible contamination by stars within NGC 5128 as well as background galaxies.

In this paper, I present a sample of globular clusters from two fields in the halo of NGC 5128. Deep images combined with high resolving power of the VLT and the proximity of NGC 5128 allow the identification of globular clusters primarily on the basis of their non-stellar PSF. Constructing the GCLF for a statistically significant number of clusters that span the whole range of magnitudes, I pursue the question of the uniqueness of GCLF in elliptical galaxies.

This paper is organized as follows. The data, the photometric calibrations and completeness tests are described in Sect. 2. In Sect. 3 the color–magnitude diagram for the NGC 5128 globular clusters is compared to the ones of the MW and M31. The luminosity function is presented in Sect. 4 and the color distribution in Sect. 5. In the latter, the mean abundances are determined for the metal-poor and metal-rich populations. The results are summarized in Sect. 6.

## 2 The Data

### 2.1 Observations and Data Reduction

I used the VLT Antu (UT1)+FORS1 (FOcal Reducer/low dispersion Spectrograph) to obtain the Bessel U- and V-band photometry of globular clusters and stars in two different fields in NGC 5128. The field of view was $`6\stackrel{}{.}8\times 6\stackrel{}{.}8`$ and the pixel scale $`0\stackrel{}{.}2/`$pixel. The FORS1 detector is a 2048$`\times `$2048 CCD with 24$`\mu `$m pixels.

Field 1 was centered on the prominent N-E shell $``$14 away from the center of the galaxy, while Field 2 was chosen to overlap with Soria et al. (soria (1996)) HST observations at a distance of $``$9 from the center of the galaxy. Observations were carried out on the 11th and 12th of July 1999 in service mode. Pairs of images of 15 min exposure time were taken per field in each filter. The journal of observations is given in Tab. 1.

In addition to VLT data, I obtained a 45 min long exposure of Field 2 in Ks using SOFI at the New Technology Telescope (NTT) at the ESO/La Silla Observatory. The field of view of SOFI is $`4\stackrel{}{.}94\times 4\stackrel{}{.}94`$ and the pixel scale $`0\stackrel{}{.}292/`$pixel. This field is smaller than the one of FORS1 and so only about 40% of objects in one field have Ks-band magnitudes.

Both fields were also observed with the Wide Field Imager (WFI) at the 2.2m telescope on La Silla. The data from WFI observations are used here to assess the V-band magnitude of some of the brightest cluster candidates that were saturated on V-band FORS1 images. I calibrated the WFI photometry by comparing magnitudes of globular clusters well exposed (but not saturated) on VLT and WFI images in order to determine the zero point of WFI observations.

For optical data (VLT and 2.2m images), the standard image processing, including the overscan and bias subtraction and the flat-field correction, was performed within the IRAF111IRAF is distributed by the National Optical Astronomy Observatories, which is operated by the Association of Universities for Research in Astronomy, Inc., under contract with the National Science Foundation environment. The pairs of exposures in each filter were registered using the imalign task and averaged to obtain the final images. Cosmic rays were rejected in this process using the crreject algorithm within the imcombine task in IRAF. Bad pixels were masked out.

The standard procedure in reducing infrared (IR) data consists of sky subtraction, flat-field correction, registering and combining the images. Good sky subtraction in a crowded field like that of a galactic halo is particularly important. For that step I used the DIMSUM package (Stanford et al. dimsum (1995)) within IRAF. In DIMSUM the sky subtraction is made in two passes. In the first one a median sky is computed for each image from 6 temporarily closest frames. The shifts between the sky subtracted frames are then computed and all the images stacked together using a rejection algorithm to remove cosmic rays. An object mask is computed for the coadded image and then shifted back in order to create object masks for the individual frames. In the second pass, the sky subtraction is made using the object masks to avoid overestimation of the sky level. These masks are also used to check that the bright object cores were not removed as cosmic rays in the previous pass. After the mask-pass sky subtraction, the frames are registered with imalign and combined with imcombine task in IRAF.

### 2.2 Cluster identification and photometry

NGC 5128 is close enough that its globular clusters can be distinguished from stars on the basis of their slightly non-stellar appearance on high resolution ground-based images taken in excellent seeing conditions. At the distance of NGC 5128 (3.6 Mpc) a typical MW globular cluster with a mean King core radius of $`r_c=2.3`$ pc had, on the best seeing VLT images, a FWHM of 2.6 pix, which is $`0.4`$ pix bigger than the stellar FWHM. The excellent seeing, especially on the V-band images of Field 2, allowed me to identify globular cluster candidates with $`r_c`$ as small as 1.7 pc, corresponding to observed FWHM of 2.4 pix (see f2.GC-6 in Tab. 4), purely on the basis of their non-stellar appearance through the following procedure within IRAF: (i) first the daofind algorithm was run in order to detect all the stars and clusters in the image; (ii) photometry through a 3 pixels aperture was performed on all the objects; (iii) relatively bright, isolated, stellar objects were used to create the PSF, which was then fitted to all the objects, subtracted and the residual images created. DAOPHOT II (Stetson stetson (1987)) permits creation of the variable PSF, which was necessary for the U-band images that showed fixed pattern residuals upon subtraction of the non-variable PSF. In order to create a variable PSF I used 35 stars per image.

Because of their slightly larger FWHM and non-stellar PSF, globular clusters were easily detected on residual images created by subtracting the objects fitted with a stellar PSF. This subtraction was performed using the allstar task. Globular cluster residuals were oversubtracted in the center and undersubtracted in the wings (Fig. 1). Images were visually inspected and the residuals in the U- and V-band images compared. This was necessary to discard obvious galaxies from the candidate list.

After visual inspection of the residual images, 50 and 81 objects in Field 1 and 2, respectively, were retained as candidates for globular clusters. Good globular cluster candidates, the ones that have round ring-like residuals, are listed in Tab. 3. Other extended objects for which the visual distinction was not very clear, some of them being blended stars or round galaxies, but some of which might be globular clusters as well, are listed in Tab. 4. The Sextractor programme (Bertin & Arnouts sextractor (1996)) was used to measure precise coordinates, FWHM ($`f_U`$ and $`f_V`$ columns in Tab. 3 and 4) and ellipticities ($`ϵ_V`$ and $`ϵ_U`$ columns) for these objects.

I performed aperture photometry within IRAF for all the objects. Since the fields were quite crowded, stellar objects had to be removed in order to obtain more precise photometry of globular clusters. To do so, I first subtracted all the stars (this time excluding the globular cluster candidates from the aperture photometry list) from the original images using the allstar task. Aperture photometry of globular cluster candidates was then measured through circular apertures of 16, 20 and 25 pixel radii. In some cases a bright saturated star or another globular cluster candidate was found near the cluster, or the latter was located near the edge of the image and thus only a 3 pix aperture magnitude was used. The 3 pix magnitudes were later extrapolated to 25 pix magnitude values using an aperture correction calculated as a median correction for more isolated globular clusters. The error in magnitude for these clusters is larger and hard to quantify, since the aperture correction depends on the concentration and compactness of the cluster. The positions and magnitudes of globular clusters are presented in Tab. 3 and Tab. 4. The clusters with magnitudes measured through a smaller aperture and later corrected are flagged with “b” (near the bright star or another globular cluster) or “e” (near the edge of the image).

### 2.3 Photometric calibration

For the photometric calibration of the VLT images, a total of 14 stars in 4 different Landolt (landolt (1992)) fields, spanning the color range $`1.321<(UV)<4.162`$, were observed during the two nights. The following transformations were calculated:

$`u_{inst}`$ $`=`$ $`U24.262(\pm 0.089)+0.379(\pm 0.068)X`$ (1)
$`0.042(\pm 0.007)(UV)`$
$$v_{inst}=V27.348(\pm 0.042)+0.213(\pm 0.034)X$$
(2)
where X is the mean airmass of the observations, $`u_{inst}`$ and $`v_{inst}`$ are the instrumental magnitudes and U and V magnitudes from the Landolt (landolt (1992)) catalogue. The one sigma scatter around the mean was 0.031 mag for the U band and 0.023 mag for the V band (Fig. 2 left panel). Adding the color term (U$``$V) in the transformations slightly reduces the scatter in the V-band to 0.017 mag (Fig. 2 right panel). However, the calibration equation without the color term for the V-band was prefered, since the data in that filter go much deeper and U magnitudes for some objects could not be measured accurately enough.

SOFI images were calibrated using 6 standards from the Persson et al. (persson (1998)) list observed during the night and spanning an airmass range from 1.1 to 3.4. Each standard star was observed at 5 different positions on the IR-array. In this way total of 30 independent measurements was obtained. The least square fit (Fig. 3) yielded the following calibration equation:

$$k_{s,inst}=K_s22.34(\pm 0.01)+0.058(\pm 0.007)X$$
(3)
### 2.4 Completeness

I made completeness tests for detection of stars on the U- and V-band images. Using the IRAF task addstar within the DAOPHOT package I added 400 stars each time to each of the U- and V-band images, all with the same magnitude and uniformly distributed over the whole field, and re-computed the photometry. The completeness and the accuracy of photometry were measured at 0.5 mag bins. The 90% completeness for stellar sources is reached at magnitude 25 in the V-band and 23.5 in U-band. The completeness drops to 50% at V$`26`$ and U$`24.5`$. More detailed results of completeness tests will be presented in the forthcoming paper where the stellar photometry of NGC 5128 is analyzed (Rejkuba et al. 2001b ).

Globular clusters have larger PSFs than stars, their selection is based mainly on the recognition of a “globular-like” residual on a subtracted image and thus their identification as such is more subjective. In order to check the completeness level of globular cluster detections, I made similar simulations as for the stars. The fainter globular clusters, that are of interest in simulation, are not resolved and thus I just treat them as stars. To add a correct, typical globular cluster PSF, the best cluster candidates’ PSFs were used to create simulated images. In order to investigate the dependence of detection on core radius, simulations with different typical PSFs were made. Two sets of simulations were done using the PSFs of clusters with small core radii, f1.GC-17 and f2.GC-4, to add objects, while in the other two sets much less concentrated clusters f1.GC-10 and f2.GC-36 were used. At each pass 100 clusters were added, distributed uniformly all with the same magnitude. The completeness of the globular cluster detection was measured for each half magnitude.

The detection depends quite significantly on core radius in the sense that more small $`r_c`$ clusters are lost at faint end. The incompleteness magnitude is defined as magnitude where more than 50% of the added clusters are not detected as such (they are actually all detected, but their residuals are not recognized as globular cluster-like, because they are lost in noise). It sets in at V$``$22.7 mag for clusters like f1.GC-17 and f2.GC-4, but most of the larger clusters, like f1.GC-10 and f2.GC-36 can still be detected at 90% level at magnitude V=23. The incompleteness for clusters with large $`r_c`$ sets in at magnitude V=23.5. In the simulations clusters could be detected with magnitudes as faint as V=24. However, at those magnitude limits it is almost impossible to distinguish globular cluster residuals from galaxies.

From the simulations I conclude that clusters fainter than $`V23`$ mag would probably not be detected. MW clusters corresponding to these magnitudes are Pal 1, Pal 13, Terzan 1, E 3 and AM 4. However, at magnitude $`V=22.0`$ ($`2`$ magnitudes past the GCLF peak) more than 90% of the clusters are detected.

### 2.5 Comparison with published data

The most recent studies of the GCS in NGC 5128 (Holland et al. holland (1999), Alonso & Minniti am97 (1997)) have investigated the central regions of the galaxy away from the fields studied here. The only photometry of globular clusters found in common with the data presented here is the photographic plate measurements of five clusters in Hesser et al. (hetal86 (1986), hetal84 (1984)) and the Washington photometry of the three of these 5 clusters in Harris et al. (hetal92 (1992)). The clusters in common are C36, C37, C38, C12 and C15 (numbering from Hesser et al. hetal84 (1984), hetal86 (1986)). While in the Hesser et al. papers only the $`B`$ magnitude is given, Harris et al. (hetal92 (1992)) report also photographic V-band magnitudes as well as Washington photometry for C12, C36 and C37. In order to compare the Washington photometry with Bessell V magnitudes, I used the transformation equation (Harris et al. hetal92 (1992)): $`V=T_1+0.66(MT_1)`$.

The clusters in common are among the brightest ones and in this study their magnitudes were measured on 2.2m+WFI images. The comparison between the data from the literature and my measurements is given in Tab. 2. In columns 1 and 2 the identification of the objects in common is given, $`B_{pg}`$ magnitudes are from Hesser et al. (hetal84 (1984), hetal86 (1986)), $`V_{pg}`$ and $`V_{Wash}`$ are photographic magnitudes and the V-band brightness derived from Washington photometry. The last two columns present my measurements. There is a good agreement between the V-band magnitude derived from the Washington photometry and my measurements. The largest difference is for cluster C12 (f2.GC-81), which I measured 0.2 magnitudes fainter, probably due to the vicinity of a bad column on the WFI chip.

### 2.6 Selection of globular clusters

Typical Galactic globular clusters have mean King core radii $`r_c=2.3\pm 3.6`$pc, mean tidal radii $`r_t=45.1\pm 32.4`$pc and mean ellipticities $`ϵ=0.07\pm 0.01`$ (June 1999 update of Harris harris96 (1996) catalogue; White & Shawl rwhite (1987)).

Holland et al. (holland (1999)) used HST WFPC2 data to find globular clusters in the central parts of NGC 5128 and measured their structural parameters. They did not find any significant difference between the MW and NGC 5128 globular clusters’ King core radii ($`r_c`$), tidal radii ($`r_t`$) and half-mass radii ($`r_h`$). They also found no trend in $`r_c`$, $`r_t`$ or $`r_h`$ with galactocentric distance, contrary to what one may expect due to stronger tidal forces closer to the galactic center. However, they measured higher mean ellipticities for their globular cluster candidates with respect to the MW globular clusters. Their samples are restricted only to the brighter clusters in NGC 5128, corresponding to the absolute magnitude $`M_V=6.4`$. For comparison, the present sample reaches almost 2 magnitudes deeper.

I measured the FWHM and ellipticity for all the objects in Tab. 3 and 4 using the Sextractor programme (Bertin & Arnouts sextractor (1996)). U-band images had variable PSF across the frame and the objects were, in some places, elongated, thus showing higher ellipticity than in the V-band images.

The FWHM vs. ellipticity for all the objects in Fields 1 and 2 is plotted in Fig. 4. The upper panel shows the results obtained from the V-band, the lower panel the ones obtained from the U-band. The vertical line indicates the mean FWHM for stars. The objects from Tab. 4, some of which might as well be blended stars or galaxies are plotted as open symbols, while the filled dots were used for objects from Tab. 3 (i.e. good candidates). Objects that were slightly saturated on the V-band VLT images are shown only in the U-band diagram (lower panel). Due to the poorer seeing conditions during the observations with the 2.2m telescope, their profiles are more similar to the stellar profiles on the WFI images.

In order to classify an object as a globular cluster all of the following criteria were imposed:

1. non-stellar PSF that leaves ring-like residuals on subtracted images (Fig. 1);
2. FWHM$`>`$mean stellar FWHM and ellipticity(V)$`<`$0.4 (Fig. 4)
3. $`0<(UV)_0<2.5`$ (Fig. 6)

The selection of objects with ellipticity measured on the V-band images smaller than 0.4 excludes galaxies from the sample. The color range in the last criterion corresponds to the observed range of the (U$``$V)0 colors for the MW and M31 globular clusters. This range does not include the reddest M31 clusters (Fig. 6) because they were corrected only for the foreground average extinction towards M31 and thus might still be reddened. The cutoff at (U$``$V)$`{}_{0}{}^{}=2.5`$ takes into account the bulk of M31 clusters around $`V_021`$. The objects that did not pass all of the above criteria were flagged with “n” in Tab. 3.

## 3 Color-Magnitude and Color-Color Diagrams

The UV color-magnitude diagram (CMD) for all the objects is presented in Fig. 5. As before, the objects retained as good candidates after visual inspection of the images are plotted as filled symbols and the rest as open symbols. The objects from Field 1 are plotted as diamonds, the ones from Field 2 as triangles. Vertical lines indicate the color range within which the MW and majority of M31 globular clusters are located. Note that several open symbols are found in the same area of the CMD as globular clusters. Spectroscopic observations are necessary to assess the real nature of these objects.

The object at (U$``$V)$`{}_{0}{}^{}=0.99`$ and $`V_0=19.96`$ (f1.GC-8) is located very close to the area where a large number of bright young blue stars were found (Rejkuba et al. 2001a ). It might be a young object, precursor of a globular cluster, similar to the ones found in other merging galaxies (e.g. Whitmore et al. whitmore (1999)), as well as to some clusters found in the central regions of NGC 5128 (Holland et al. holland (1999)).

There is a number of objects brighter than $`V<17.5`$ that were all identified in Field 1. Probably these are bright stars blended with some faint companions. They all appear saturated in the V-band VLT images and were found in the lower right corner (south-west) of Field 1 on U-band images. Their FWHM, that was larger than stellar, is probably an artifact due to the distortion of the PSF in that part of the chip.

Figure 6 shows the CMD for all NGC 5128 globular clusters (filled circles) that satisfy all the selection criteria. This CMD is compared with the ones of the MW (stars) and M31 (squares) globular clusters. Globular clusters populate a similar region of the CMD in the three galaxies. Data for the MW clusters were taken from the web page of W. Harris (harris96 (1996)). M31 clusters are from the catalogue of Barmby et al. (barmby (2000)). I did not take into account globular clusters with colors bluer than B$``$V$`=0.55`$, as they might be younger objects (Barmby et al. barmby (2000)). The magnitudes were corrected for the mean reddening towards M31 (E(B$``$V$`)=0.08`$). For the globular cluster candidates in NGC 5128 a mean foreground reddening of E(B$``$V$`)=0.1`$ (Burstein & Heiles BH (1982); Schlegel et al. schlegel (1998)) was used. The reddening vector plotted in the upper right corner corresponds to E(B$``$V$`)=0.2`$.

By constraining the colors to be within the range of the reddest MW and bluest M31 globular cluster colors (cf. the third selection criterion) I exclude distant compact galaxies (Barrientos barrientos (2000)). The total number of bona fide globular cluster candidates is 71, with 23 in Field 1 and 48 in Field 2. Only those are taken into account in the further analysis.

Ks band magnitudes were measured on SOFI images for 23 of the bona fide globular cluster candidates and 10 of the extended objects, some of which might also be globular clusters. The U$``$V color is strongly dependent on both metallicity and age (Buzzoni buzzoni (1995); Kurth et al. kurth (1999); Maraston maraston98 (1998)). Adding the K-band magnitude information it is possible to disentangle the young from the old objects. The color-color plot is shown in Fig. 7. Filled dots are used to denote globular cluster candidates that satisfied all of the classification criteria. The lines are isochrones of evolutionary synthesis models of simple stellar populations (SSP) calculated for a Salpeter IMF and a range of metallicities from Z=0.0001 to 0.05 (see caption for the exact metallicities plotted for each model). The four different models displayed on panels (a), (b), (c) and (d) are from Kurth et al. (kurth (1999)), Bruzual & Charlot (bc99 (2001)), Maraston (maraston00 (2001)) and Worthey (wo94 (1994)), respectively. The ages are indicated with numbers denoting Gyr next to the lines. Unfortunately the models do not allow one to discriminate between ages in the range of 5-20 Gyr. However, it is comforting to see that almost all the objects are concentrated towards the isochrones in excess of 1 Gyr. This is additional evidence that they are indeed globular clusters. A spread of metallicities of more than 1.5 dex is also evident.

In Fig. 7 most of the clusters lie slightly to the right and below the model lines. There are two possible explanations. First, a difference with the SSP models, particularly significant in U-V and B-V colors, has already been pointed out by Barmby & Huchra (bh2000 (2000)). They suggested U fluxes from the stellar libraries, systematic problems with modeling HB color and observational errors as possible causes of this offset. Second, in the data presented here additional offset may arise from the difference between the Bessel and Johnson U-band transmission curves (Bessel bessel (1995)).

## 4 Luminosity Function

The V-band GCLF for NGC 5128 is compared to the ones of the MW and M31 in Fig. 8. The actual shape of the luminosity function depends on the size of the bins and on the bin-centers. Statistically it is possible to compare them independently of binning using the Kolmogorov-Smirnov (KS) test. In Fig. 9 the cumulative distribution functions for the three GCLFs are shown. The thick line is used for NGC 5128, thin line for the MW and dotted line for M31. Absolute magnitudes for NGC 5128 globular clusters were calculated using the distance modulus value of (m$``$M$`)_0=27.8`$ (Soria et al. soria (1996)) and reddening value of E(B$``$V$`)=0.1`$.

The KS statistic measures the maximum value of the absolute difference between two cumulative distribution functions from which a probability that the two distributions are drawn from a common distribution can be calculated. The probability that the GCLF of NGC 5128 matches the one of the MW is 89%, comparable to the probability that M31 GCLF matches the MW one (85%). The difference is slightly bigger in the U-band distributions (Fig. 9 right panel), but it is also probable that the three luminosity functions are essentially similar for the three galaxies. I conclude that the NGC 5128 GCLF is as similar (or as different) to the MW GCLF as the M31 one.

Usually the GCLF is fitted with a gaussian or $`t_5`$ distribution function in order to determine its peak magnitude. It should be noted, however, that the reasonably complete luminosity functions are known only for the two spiral galaxies, the MW and M31. Both of these GCLFs appear to depart from the gaussian distribution (Ashman et al. acz (1995)), while there is no complete GCLF for any elliptical galaxy. The dependence of the turn–over magnitude on the detection of the faint end of GCLF is discussed in Ashman et al. (acz (1995)). The determination of the turn–over magnitude for NGC 5128 is left for a subsequent paper where the complete analysis of WFI data is going to be presented.

The detection of the faint end of the GCLF in an elliptical galaxy is also important from the dynamical point of view. The faint globular clusters are typically also the less massive ones and their presence in the deep potential of a giant elliptical puts strong constraints on the potential itself as well as on the effects of tidal forces in the halo. Unfortunately the faintest globular clusters, corresponding to the MW Pal 1, AM 4, Terzan 1, Pal 13 are probably not detectable on FORS1 images (see Sect. 2.4).

The processes that destroy globular clusters can be divided into internal, like two-body relaxation and mass loss during stellar evolution, and external. The external processes depend greatly on environment and therefore on the position of the clusters in the galaxy. While tidal shocks preferentially destroy less dense clusters, dynamical friction brings more massive clusters very close to galactic centre. Both of these external processes are much stronger in the inner part of the galaxy. According to theoretical simulations (Vesperini vesperini (2000), Ostriker & Gnedin ostriker\_gnedin (1997)) in a galaxy like NGC 5128, with the total mass of $`4\times 10^{11}`$ M (Israel israel (1998)) and R$`{}_{e}{}^{}=5.24`$ kpc (Dufour et al. dufour (1979)), a strong disruption and mass loss are efficient only within a galactocentric distance of $`12`$ Re, while at larger radii most of the mass loss is due to effects of stellar evolution. Field 1 is located at $`23.5`$ Re and Field 2 at $`1.22.3`$ Re. While the dynamical effects may still be important for Field 2 clusters, they are probably completely negligible for clusters in the outer field. It is not possible to quantify these effects by simple comparison of the two subsamples, because of the small number statistics.

## 5 Color and Metallicity Distributions

The (U$``$V)0 color histogram for 71 globular clusters in NGC 5128 indicates a bimodal distribution (Fig. 10). Note that the (U$``$V)0 color distribution of the MW globular clusters does not show bimodality while the one of M31 does. The (U$``$V)0 color distribution of the MW globular clusters is in fact biased due to lack of U-magnitude and higher reddening values (only clusters with E(B$``$V$`)<0.5`$ are plotted) of bulge (more metal-rich and redder) cluster. The blue peaks of the MW, M31 and NGC 5128 globular cluster color distributions coincide. The blue cutoff in the three galaxies is quite sharp and is driven principally by the colors of blue horizontal branch stars. The red side of the color distribution is much broader. In M31 it can be partially driven by differential reddening inside the host spiral galaxy, because the colors for M31 clusters presented here were corrected only for the mean foreground reddening of E(B$``$V)=0.08. Examination of the stellar (U$``$V) vs. (V$``$K) diagram in the two fields in NGC 5128 indicates that the amount of reddening does not vary much in the halo of this giant elliptical (Rejkuba et al. 2001b ). Harris et al. (hetal92 (1992)) and Holland et al. (holland (1999)) also found more red (metal-rich) clusters in their NGC 5128 data.

The color distribution of the MW globular clusters in Fig. 10, being dominated by halo clusters, is representative of the GCS of our Galaxy as if observed from outside. The NGC 5128 clusters in Harris et al. (hetal92 (1992)) and in this work belong as well to the halo population. The differences in the color (metallicity) distribution of clusters and stars in the halo of a spiral, like MW and a gE like NGC 5128, are reflected in figure 10 with a much larger metallicity spread in gE galaxy.

The line overplotted on the color distribution of NGC 5128 is the best nonparametric kernel estimate (Silverman silverman (1986)) of the color distribution used to measure accurately the blue and the red peaks of the distribution. The blue and red peaks are measured at (U$``$V)$`{}_{0}{}^{}=0.73\pm 0.01`$ and $`1.28\pm 0.01`$, respectively.

From the color-color diagram (Fig. 7) it is obvious that the globular clusters in NGC 5128 are older than 1 Gyr. Unfortunately, the models do not allow one to better constrain the ages (Fig. 7; see also Barmby & Huchra bh2000 (2000)). Assuming that the globular clusters in NGC 5128 have similar ages as the MW ones, one can use the Galactic globular clusters to calibrate the metallicity scale. It is known that for the very low and very high metallicities the color-metallicity relation cannot be expressed by a simple linear function, but extrapolation from lower metallicities towards the higher ones by using the higher polynomial fit might lead to unphysical results. Thus I prefered to make a linear fit with gaussian errors in $`\sigma (Fe/H)=0.1`$ in order to calibrate the metallicity as a function of intrinsic (U$``$V)0 color. Only MW globular clusters from W. Harris (harris96 (1996)) web list with E(B$``$V)$`<0.5`$ (equivalent to E(U$``$V)$`<0.82`$; Rieke & Lebofsky rl85 (1985)) were used in the fit. The following color-metallicity relation is based on the Zinn & West (ZW (1984)) scale and is valid only for the range of metallicities observed in the MW globular clusters (Fig. 11):

$$[Fe/H]=2.93(\pm 0.07)+1.79(\pm 0.07)\times (UV)_0$$
(4)
In the Fig. 11 linear (U$``$V)0 vs. \[Fe/H\] color-metallicity relation (Eq. 4; full line) is plotted along with the evolutionary synthesis models of simple stellar populations (short dashed lines) from Kurth et al. (kurth (1999)). Using the empirical linear calibration and the SSP models I calculate the metallicity of the two peaks in the color-distribution. The blue and the red peak correspond to \[Fe/H\]$`=1.7`$ dex and \[Fe/H\]$`=0.6`$ dex using the linear relation (Eq. 4) and to \[Fe/H\]$`=1.7`$ dex and \[Fe/H\]$`=0.5`$ dex using the models.

Although it is not obvious from the (U$``$V)0 color distribution, but rather from the spectroscopically measured \[Fe/H\] distribution, the metallicity distribution of MW globular clusters is bimodal with the metal-poor and metal-rich peaks at -1.6 and -0.6 dex respectively (Harris harris2000 (2000)). These values are practically the same as the ones derived here for globular clusters in NGC 5128 from the (U$``$V)0 colors. On the other hand, on the basis of Washington photometry of 62 globular clusters in NGC 5128, Harris et al. (hetal92 (1992)) conclude that “in both the range and distribution of abundance, the NGC 5128 GCS is clearly different from the Galactic GCS and quite similar to that of the large elliptical NGC 1399”. They measured a mean metallicity of $`[\mathrm{Fe}/\mathrm{H}]_{CT_1}=0.8\pm 0.2`$ dex, substantially more metal-rich than the mean value for the MW GCS. They, however, mention that the values exceeding \[Fe/H\]$`0.25`$ are extrapolated, a fact that is also true in the metallicity calibration from (U$``$V)0 color distribution. The extrapolation leads to large errors and a possible overestimation of abundances of the metal-rich clusters. Because of this I do not report here the metallicities of individual clusters.

Note also that I restricted the (U$``$V)0 color to be within the range observed in MW and bulk of the M31 globular clusters. In this way the sample was selected to be as clean as possible from the very compact and red galaxies. It is not surprising that the clusters with supra-solar abundances are not present, reducing the mean abundance of the sample. However, relaxing the red color cut and adding the 6 “good” globular cluster candidates redder than (U$``$V)$`{}_{0}{}^{}=2.5`$ (filled symbols in Fig. 5) the mean of the color distribution shifts from $`(\mathrm{U}\mathrm{V})_0=1.05\pm 0.45`$ to $`(\mathrm{U}\mathrm{V})_0=1.22\pm 0.75`$, corresponding to mean metallicities of $`[\mathrm{Fe}/\mathrm{H}]=1.1\pm 0.1`$ dex and $`[\mathrm{Fe}/\mathrm{H}]=0.7\pm 0.2`$ dex, respectively. The latter value is very similar to the mean metallicity derived from Washington photometry by Harris et al. (hetal92 (1992)). In order to assess the real nature of the objects redder than (U$``$V)$`{}_{0}{}^{}=2.5`$ spectroscopic observations are necessary.

## 6 Conclusions

The main results and their implications presented in this work are the following:

* On the basis of high resolution ground-based images, taken with VLT and FORS1, I identified 71 bona fide globular clusters in the halo of the nearest giant elliptical galaxy NGC 5128. Only 5 of these clusters were previously known.
* In a UV color-magnitude diagram for the complete sample of 71 clusters and a UVK color-color diagram for a subsample of 23 clusters, the objects span a similar magnitude and color range as the globular clusters in the MW and M31. Note, however, that there is no cluster corresponding to $`\omega `$ Cen, the brightest MW cluster, in these two fields. The position of 23 objects in the (U$``$V)0 vs. (V$``$K)0 color-color diagram is consistent with their classification as old globular clusters.
* The GCLFs spanning $`10.1<M_V<4.9`$ and $`9.3<M_U<3.3`$ have been constructed. These are the deepest GCLFs of an elliptical galaxy made so far. Kolmogorov-Smirnov statistics show that the difference between the GCLFs of NGC 5128 and MW is not larger than the difference between the GCLFs of M31 and MW. Similarity of the GCLFs of an elliptical with respect to the spiral galaxy had never been tested at the faint end before.
* The presence of faint globular clusters in the halo of NGC 5128 puts constraints on the effectiveness of the tidal forces in the deep elliptical galaxy potential. The dynamical effects may be important for clusters that are found within $`2`$ Re from the galactic center. Unfortunately, the selection of globular clusters from VLT images is not sensitive to the most compact and the faintest clusters, similar to ones like Pal 1, Pal 13, AM 4 and Terzan 1 in our Galaxy, while the less dense ones start to be confused with background galaxies at faint magnitudes.
* The (U$``$V)0 color histogram of 71 clusters indicates a bimodal distribution, supporting the Zepf & Ashman (zepfashman93 (1993)) suggestion. Assuming that the clusters in NGC 5128 and in the MW span a similar age interval, and adopting the linear fit between the (U$``$V)0 color and metallicity, I derived the \[Fe/H\] of the red and blue peaks of the bimodal distribution to be $`1.7`$ dex and $`0.6`$ dex, respectively. Using the SSP models from Kurth et al. (kurth (1999)), instead of the linear fit, the results do not change significantly, giving values of $`1.7`$ dex and $`0.5`$ dex for the metal-poor and metal-rich peaks, respectively. This is different from the Harris et al. (hetal92 (1992)) and Zepf & Ashman (zepfashman93 (1993)) result, but is partially due to the adopted color cut of (U$``$V)$`{}_{0}{}^{}<2.5`$ in selecting globular clusters and to small number statistics, since both samples have $`<10`$% of the total cluster population. Relaxing the red color cut the mean color of the distribution corresponds to mean metallicity of \[Fe/H$`]=0.7`$ dex, very similar to Harris et al. (hetal92 (1992)) result. The true nature of the redder objects has to be assessed through spectroscopy.

###### Acknowledgements.

I am grateful to Dante Minniti for his guidance and help with this project. I thank Andres Meza for providing the code that calculates the kernel estimator. I would like to acknowledge the referee for valuable comments which led to improvements in the paper and useful discussions with Markus Kissler-Patig, Thomas Puzia and Steve Zepf. Thanks also go to Tim Bedding, Dave Silva and Elena Pancino. This research was partially supported by FONDECYT grant No. 1990440.


format: text

text:
# Introduction

## Introduction

One of the key issues in deeply inelastic scattering and deeply inelastic diffraction is whether parton saturation effects are present in the HERA energy regime. The answer to this question has important implications for our understanding of the very early stages of relativistic heavy ion collisions where parton saturation would result in the production of a high density and high field strength, $`F_{\mu \nu }1/\sqrt{\alpha _s},`$ state of QCD. Such a state would be a new and exciting regime of nonperturbative but small coupling QCD.

Parton saturation can be viewed as a manifestation of unitarity limits being reached. However, unitarity limits are difficult to see directly in deep inelastic scattering and diffraction. Thus, most studies of parton saturation at HERA have been done in terms of models which explicitly impose unitarity at moderate $`Q^2`$ and small $`x`$. The most widely used such model is that of Golec-Biernat and Wüsthoff  where the $`S`$-matrix for the scattering of a dipole of separation $`r`$ on a proton is given by $`S=e^{r^2Q_s^2/4}`$ with the saturation momentum, $`Q_s^2,`$ depending on the energy of the scattering. The fact that the Golec-Biernat Wüsthoff model and its generalization work so well at HERA is perhaps some evidence that deep inelastic scattering, and especially diffraction, have reached unitarity limits.

It would be nice to see the attainment of unitarity limits directly without having to use a model. This is not an easy task. Recall that in proton-proton scattering the energy dependence of total and elastic cross sections changes little between ISR and Tevatron energies. The near saturation of unitarity limits was found in proton-proton scattering only when Amaldi and Schubert  evaluated the proton-proton elastic scattering amplitude in impact parameter space and found near blackness for small impact parameters.

Vector meson production seems to be the best process from which to extract the dipole proton elastic scattering amplitude. For example, it is easy to check that our procedure (see Eqs. (4)-(12) below) would not work for the total diffractive cross section because there is not a single definite state produced, a condition which seems necessary in the way we proceed.

In this paper we copy the Amaldi-Schubert analysis in order to determine the dipole-proton scattering amplitude in impact parameter space in terms of diffractive $`\rho `$production data at HERA. Recall that diffractive electroproduction of a $`\rho `$meson can be viewed in terms of the following sequence of transitions: (i) The virtual photon goes into a quark-antiquark pair (dipole). (ii) The dipole scatters elastically on the proton. (iii) The dipole then goes into a $`\rho `$meson. Thus $`\rho `$electroproduction naturally carries information on the dipole-quark elastic scattering amplitude. The main uncertainty is in the wavefunction of the $`\rho `$meson which appears to be reasonably well constrained by previous phenomenology.

The dipole-proton $`S`$-matrix which we extract as a function of impact parameter is averaged over the dipole sizes which naturally appear in the $`\gamma ^{}`$ and $`\rho `$meson wavefunctions. This averaging does not affect the fact that for weak interactions $`S`$ is near 1 and for strong interactions $`S`$ is near zero. Thus the smallness of $`S`$ is a measure of how close one is to the unitarity limit, $`S=0.`$ Perhaps the best way to measure how close the dipole-proton scattering is to the unitarity limit is to note that $`1S^21S^2`$ gives the probability that a dipole passing the proton will induce an inelastic reaction at the impact parameter in question ($`\mathrm{}`$ denotes the average over dipole sizes as explained in Eq. (11)). From our analysis this probability is likely significantly greater than $`1/2`$ for $`Q^22\text{GeV}^2`$ and for $`b`$ near zero (see Fig. 6), although better large momentum transfer data would be needed to definitively determine $`S`$ at $`b0.`$ Thus, it is likely that saturation effects are important in this low $`Q^2`$ regime. Our estimate of $`Q_s^2`$ is about $`11.5\text{GeV}^2`$ at $`b0.3`$ fm, although this estimate has large uncertainties because it depends on knowing the average dipole size as a function of $`Q^2.`$ Finally, we determine the total dipole-proton cross sections for $`x10^4`$ to be 14.9, 10.6 and 7.5 mb at $`Q^2=0.45,3.5`$ and $`7\text{GeV}^2`$ respectively.

The outline of the paper is as follows: in section 1 we establish the relation between the $`S`$-matrix for dipole-proton scattering averaged over the distribution of dipole sizes and the diffractive differential cross section for electroproduction of vector mesons. Due to the particular properties of this distribution, this formula translates into a method of determining the $`S`$-matrix profile in impact parameter space for a dipole of given size $`r`$. The only theoretical input needed is the vector meson wave function. In section 2 we discuss the parametrizations for the wave functions that we used, and we study in detail the model independence of our procedure. In section 3 we apply this method to the analysis of the available HERA data and present our main results on the $`S`$-matrix profile in impact parameter space. We discuss the uncertainties due to the experimental errors on the measured cross section and due to the lack of data in the high $`t`$ region. We also discuss the possible implications for studies of saturation effects at high energies. Finally, in section 4 we state the summary of our results. We append some more technical points.

## 1 Dipole-proton $`S`$-matrix and diffractive meson production cross section

We consider the process of diffractive production of a vector meson, $`\gamma ^{}pVp`$, shown in Fig. 1. We adopt here the dipole formulation  in which the virtual photon $`\gamma ^{}`$ fluctuates into a $`q\overline{q}`$ pair (dipole) which subsequently interacts elastically with a proton, and eventually forms a meson bound state. This picture is justified at high center-of-mass energy of the $`\gamma ^{}p`$ system, since in this regime it is guaranteed that these successive processes are clearly separated in time. We assume that the energy is sufficiently high so that $`s`$-channel helicity conservation holds to a good accuracy, and we limit ourselves to the case of longitudinal vector meson production. The kinematics of this process is shown in Fig. 1: we define the 2-momentum transfer $`𝚫`$ related to the Mandelstam variable $`t`$ by $`t=𝚫^2`$. $`Q^2=q^2`$ and $`x=Q^2/(2Pq)`$ are as usual the virtuality of the photon and the Bjorken variable, respectively.

In this framework, the amplitude $`𝒜_{\text{el}}(x,𝚫,Q)`$ for the scattering can be written in the following factorized form:

$$𝒜_{\text{el}}(x,𝚫,Q)=\underset{h,\overline{h}}{}d^2𝐫𝑑z\psi _\gamma ^{}^{h,\overline{h}}(z,𝐫;Q)A_{\text{el}}^{q\overline{q}p}(x,𝐫,𝚫)\psi _V^{h,\overline{h}}(z,𝐫),$$
(1)
where $`A_{\text{el}}^{q\overline{q}p}(x,𝐫,𝚫)`$ is the elementary amplitude for the elastic scattering of a dipole of size $`𝐫`$ on the proton. $`\psi _\gamma ^{}^{h,\overline{h}}`$ is the photon light-cone wave function projected onto a state made of a quark-antiquark pair of charges $`e_q`$, masses $`m_q`$ and respective helicities $`h`$ and $`\overline{h}`$. This function is computed in first order light-cone perturbation theory of quantum electrodynamics, and for the longitudinally polarized virtual photon, it reads

$$\psi _\gamma ^{}^{h,\overline{h}}(z,𝐫;Q)=\delta _{h,\overline{h}}\frac{e_q}{\sqrt{4\pi }}\frac{\sqrt{N_c}}{2\pi }2z(1z)QK_0(\epsilon r),$$
(2)
where

$$\epsilon ^2=Q^2z(1z)+m_q^2.$$
(3)
Similarly, $`\psi _V^{h,\overline{h}}`$ is the wave function of the vector meson produced in the final state. We discuss it in more detail in the next section. The variable $`z`$ in Eq. (1) is the fraction of longitudinal momentum of the virtual photon carried by the quark.

The amplitude is normalized in such a way that the differential cross section for the full process is

$$\frac{d\sigma }{dt}=\frac{1}{16\pi }|𝒜_{\text{el}}(x,𝚫,Q)|^2.$$
(4)
Furthermore, the elementary amplitude $`A_{\text{el}}^{q\overline{q}p}(x,𝐫,𝚫)`$ can be related to the $`S`$-matrix element $`S(x,𝐫,𝐛)`$ for the scattering of a dipole of size $`𝐫`$ at impact parameter $`𝐛`$

$$A_{\text{el}}^{q\overline{q}p}(x,𝐫,𝚫)=d^2𝐛\stackrel{~}{A}_{\text{el}}^{}{}_{}{}^{q\overline{q}p}(x,𝐫,𝐛)e^{i𝐛𝚫}=2d^2𝐛[1S(x,𝐫,𝐛)]e^{i𝐛𝚫}.$$
(5)
We can check briefly the consistency of this formula, which defines $`S(x,𝐫,𝐛)`$ (see Ref. ). In the standard normalizations adopted here for the amplitudes, the optical theorem takes the form $`\sigma _{\text{tot}}^{q\overline{q}p}(x,𝐫)=miA_{\text{el}}^{q\overline{q}p}(x,𝐫,𝚫=0)`$. Using formula (5), this relation translates into the following expression for the total dipole-proton cross section:

$$\sigma _{\mathrm{tot}}^{q\overline{q}p}(x,𝐫)=2d^2𝐛(1eS(x,𝐫,𝐛)).$$
(6)
On the other hand, taking the amplitude (5) squared, dividing it by the flux factor and integrating it over phase space, it is easy to see that the elastic cross section is given by

$$\sigma _{\mathrm{el}}^{q\overline{q}p}(x,𝐫)=d^2𝐛\left|1S(x,𝐫,𝐛)\right|^2.$$
(7)
When $`S=0`$, the unitarity limit is reached which is equivalent to the scattering on a black body. One sees from formulae (6) and (7) that the elastic cross section is half the total one, as it should be. This justifies the consistency of the normalization adopted for $`S(x,𝐫,𝐛)`$ in Eq. (5). We will assume in the following that $`iA_{\text{el}}^{q\overline{q}p}(x,𝐫,𝚫)`$ is purely imaginary, i.e. that $`S(x,𝐫,𝐛)`$ is real.

The main goal of this paper is to extract $`S(x,𝐫,𝐛)`$ taking advantage of the present experimental knowledge of the cross section for diffractive vector meson production, see Eq. (1). To this aim we express the amplitude which appears in Eqs. (1,4) by means of $`S(x,𝐫,𝐛)`$ using relation (5). Then we take the inverse Fourier transform of the square root of Eq.(4), which gives

$$\frac{d^2𝚫}{(2\pi )^2}\sqrt{\frac{d\sigma }{dt}}e^{i𝚫𝐛}=\frac{1}{\sqrt{16\pi }}\underset{h,\overline{h}}{}d^2𝐫𝑑z\psi _\gamma ^{}^{h,\overline{h}}(z,𝐫,Q)\mathrm{\hspace{0.17em}2}[1S(x,𝐫,𝐛)]\psi _V^{h,\overline{h}}(z,𝐫).$$
(8)
In the following we suppress the angular dependence of the $`S`$-matrix since in this process one is only sensitive to the quantities which are angular averaged. We single out the (logarithmic) distribution of dipole sizes at the photon vertex, which is the overlap between the photon and the meson wave functions

$$p(r,Q)2\pi r^2\underset{h,\overline{h}}{}_0^1𝑑z\psi _\gamma ^{}^{h,\overline{h}}(z,r;Q)\psi _V^{h,\overline{h}}(z,r),$$
(9)
and we denote $`N(Q)`$ its integral

$$N(Q)_0^{\mathrm{}}\frac{dr}{r}p(r,Q).$$
(10)
We then define the mean value of a given function $`f(r)`$ with respect to the probability measure $`p(r,Q)/N(Q)`$ by

$$f(r)_p_0^{\mathrm{}}\frac{dr}{r}\frac{p(r,Q)}{N(Q)}f(r).$$
(11)
We now rewrite Eq. (8) using these definitions (9), (10), (11), and we obtain the following formula:

$$S(x,r,b)_p=1\frac{1}{2N(Q)\pi ^{3/2}}d^2𝚫e^{i𝚫𝐛}\sqrt{\frac{d\sigma }{dt}},$$
(12)
which shows that the measurement of the differential cross section $`d\sigma /dt`$ enables us to determine the $`S`$-matrix at fixed impact parameter $`b`$, averaged over the dipole size distribution $`p(r,Q)`$ defined by Eq. (9).

One should stress that $`N(Q)`$ which appears in Eq. (12) is the only source of theoretical uncertainty in this formula for the average $`S`$-matrix. This quantity depends on the parametrization of the wave function of the vector meson. However, as we shall discuss in detail in Sec. 2, $`N(Q)`$ is very well constrained and is hardly model-dependent. We now have to investigate the meaning of this average over dipole sizes, by studying $`p(r,Q)`$ in more detail.

The distribution $`p(r,Q)`$ has some interesting properties which are quite general and rather independent of the particular model for the meson wave function: it is sharply peaked at a specific value of $`r`$, which depends on $`Q^2`$ like $`r_Q=A/\sqrt{Q^2+m_V^2}`$ and its width, roughly independent of $`Q^2`$, is of order unity on a logarithmic scale. To illustrate these properties, we choose a particular model for $`\psi _V(z,r)`$  (details will be given in the next section), and we plot $`p(r,Q)`$ for various values of $`Q^2`$ in Fig. 2. One can easily obtain a rough estimate of $`r_Q`$ considering the fact that on one hand, the photon wave function behaves like $`\mathrm{log}r`$ at small $`r`$, and on the other hand, it can be approximated by $`\mathrm{exp}(\epsilon r)`$ in the asymptotic region $`r\mathrm{}`$, see Eq. (2). Assuming furthermore that the vector meson wave function is smooth in $`r`$, the integrand $`p(r,Q)/rrK_0(\epsilon r)`$ is then peaked around the maximum of the function $`r\mathrm{exp}(\epsilon r)`$, i.e. around $`r_Q1/\epsilon `$. In the case of longitudinal vector meson production, the dominant contribution comes from symmetric configurations for which $`z`$ is close to $`1/2`$. This leads to the above-quoted formula for $`r_Q`$, with $`A2`$.

When $`S(x,r,b)_p`$ is significantly below $`1`$, the quantity $`r_Q`$ can be interpreted as the mean value of the sizes $`r`$ of the dipoles which participate in the interaction. In this case, the dipole-proton amplitude $`\stackrel{~}{A}_{\text{el}}^{}{}_{}{}^{q\overline{q}p}(x,r,b)`$ is large for any dipole of size $`r`$ present in the initial wave function and thus does not filter these dipoles: the dominant dipole sizes involved in the interaction are then effectively distributed according to $`p(r,Q)`$. However, in the opposite case when $`S(x,r,b)_p1`$, the amplitude $`\stackrel{~}{A}_{\text{el}}^{}{}_{}{}^{q\overline{q}p}(x,r,b)`$ is small for the relevant values of $`r`$. In this situation, we are sensitive to the behaviour of the dipole cross section at small $`r`$, i.e. $`\stackrel{~}{A}_{\text{el}}^{}{}_{}{}^{q\overline{q}p}(x,r,b)r^2`$. This is known as the colour-transparency property. Then the interacting dipoles are effectively distributed according to $`r^2p(r,Q)`$. This results in the dominant dipole sizes being shifted towards larger values222The quantity $`r_s`$ in this case is the so called “scanning radius” introduced in . of $`rr_s`$ for which an estimate can be provided using the same technique as previously, and yields $`r_s6/\sqrt{Q^2+m_V^2}`$.

According to this discussion and considering the fact that one is a priori more interested in the case when $`S(x,r,b)<1`$, one can make the following approximation:

$$S(x,r,b)_pS(x,r_Q,b)\text{with}r_Q\frac{A}{\sqrt{Q^2+m_V^2}},A2.$$
(13)
Thus for fixed $`Q^2`$ and energy Eqs. (12,13) enable us to determine $`S(x,r_Q,b)`$, which is the $`S`$-matrix element for the scattering of a dipole of size $`r_Q`$ on a proton, at energy $`Q^2/x`$ and impact parameter $`b`$. A novel feature of our analysis is the fact that we are able to extract the profile of the $`S`$-matrix (or equivalently of the dipole-proton amplitude) in impact parameter space. The previous studies of Refs.  provided an estimate of the dipole cross section integrated over the impact parameter.

## 2 From mesons to dipoles

As already mentioned earlier in Sec.1 the only theoretical uncertainty in the extraction of $`S`$-matrix profile in impact parameter space is the form of the vector meson wave function which occurs in the overlap function $`p(r,Q)`$ (Eq. (9)) and in its integral $`N(Q)`$ (Eq. (12)). In this section, we present the different models that we choose for the vector meson wave function, then we compare them phenomenologically, and finally we evaluate the uncertainty on $`S(x,r_Q,b)`$ induced by the freedom of choice of the vector meson wave function.

### 2.1 Models for the vector meson wave function

Several different models for the vector meson wave function exist in the literature . All these approaches use the information from spectroscopic models of long distance physics. In these models one assumes that the meson is composed of a constituent quark and antiquark which move in an harmonic oscillator potential. This results in a wave function which has a gaussian dependence on the spatial separation between the quarks. Additionally these models are supplemented by the short-distance physics driven by the QCD exchange of hard gluons between the valence quarks of the vector meson. Finally, a relativization procedure has to be applied.

Of course, there are many uncertainties in the above procedure of obtaining the wave functions for the vector mesons, however they are constrained by model independent features. First of all, $`\psi _V^{h,\overline{h}}`$ has to satisfy the following normalization condition

$$1=\underset{h,\overline{h}}{}d^2𝐫𝑑z|\psi _V^{h,\overline{h}}(z,𝐫)|^2.$$
(14)
Second, the value of the wave function at the origin is related to the leptonic decay width $`\mathrm{\Gamma }(Ve^+e^{})`$ of the vector meson given by the following formula

$$_0^1𝑑z\psi _V(z,r=0)=\sqrt{\frac{\pi }{N_c}}\frac{f_V}{2\widehat{e}_V},\text{where}0|J_{\text{em}}^\mu (0)|Ve_qf_Vm_V\epsilon ^\mu ,$$
(15)
where $`f_V`$ is the coupling of the meson to the electromagnetic current and $`\widehat{e}_V`$ the isospin factor, which is the effective charge of the quarks in units of the elementary charge $`e`$: for the $`\rho `$ meson it is the charge of the combination $`(u\overline{u}d\overline{d})/\sqrt{2}`$, i.e. $`\widehat{e}_V=1/\sqrt{2}`$. Finally, one requires that the mean radius be consistent with the electromagnetic radius of the vector meson.

In our calculation we consider the models proposed in Ref.  (which we refer to as DGKP) and  (which we denote NNZ(1994) and NNPZ(1997) respectively).

#### Model DGKP

The following wave function is adopted:

$$\psi _V^{h,\overline{h}}(z,r)=\delta _{h,\overline{h}}z(1z)\frac{1}{\sqrt{4\pi }}\frac{\pi f_V}{\sqrt{N_c}\widehat{e}_V}f(z)\mathrm{exp}(\omega ^2r^2/2).$$
(16)
The parameter $`\omega `$ and the overall normalization are fixed by the condition (14) and by the value of the leptonic width. (The exact values of these parameters as well as the form of the function $`f(z)`$ are given in appendix A). Using the above form of the wave function $`\psi _V`$ one can write $`p(r,Q)`$ by replacing it in Eq. (9). In the case of model DGKP, it reads:

$$p(r,Q)=2\pi r^2\frac{dz}{4\pi }f_Ve_qz(1z)f(z)e^{\omega ^2r^2/2}\mathrm{\hspace{0.17em}2}z(1z)QK_0(\epsilon r).$$
(17)
#### Models NNZ(1994) and NNPZ(1997)

The wave function of the vector meson is given by

$$\psi _V^{h,\overline{h}}(z,𝐫)=\delta _{h,\overline{h}}\sqrt{\frac{N_c}{4\pi }}\frac{1}{m_Vz(1z)}[m_V^2z(1z)_{r}^{}{}_{}{}^{2}+m_q^2]\varphi (r,z).$$
(18)
The form of the radial wave function $`\varphi (r,z)`$ is given in appendix A. The expression for the overlap $`p(r,Q)`$ then is

$$p(r,Q)=2\pi r^2\frac{N_c\widehat{e}_V}{(2\pi )^2}\frac{2e_qQ}{m_V}𝑑z\left\{[m_V^2z(1z)+m_q^2]K_0(\epsilon r)ϵK_1(\epsilon r)\frac{}{r}\right\}\varphi (r,z).$$
(19)
The models NNZ(1994) and NNPZ(1997) differ only by the tuning of the parameters.

Expressions (17) and (19) differ on several points. We study the influence of these differences on the distribution $`p(r,Q)`$ and on its integral $`N(Q)`$ in the next sections.

### 2.2 Comparison of the theoretical predictions with the HERA data at $`t=0`$

Before applying the models (16,17) and (18,19) for the vector meson wave functions to our procedure of extracting the $`S`$-matrix element, we additionally checked how they compared with the data on forward diffractive vector meson electroproduction at HERA.

To this aim we compute $`d\sigma /dt|_{t=0}`$ using formulae (1,4) and substitute the dipole-proton amplitude $`A_{\text{el}}^{q\overline{q}p}`$ by the one proposed in Golec-Biernat and Wüsthoff’s saturation model333 A similar study was done in detail recently in  where different forms of the wave functions for $`\rho _0`$ and $`J/\mathrm{\Psi }`$ were considered and used in the calculation of the forward longitudinal and transverse total cross sections. The meson wave functions were chosen according to general principles rather than just relying on existing models. ($`\sigma (x,r)`$ in the notations of Ref. ).

The theoretical curves corresponding to all three different models DGKP, NNZ(1994), NNPZ(1997) are plotted in Fig. 3. Let us stress that there is no additional tuning of any of the parameters of the wave functions. One of the models, NNZ(1994), is able to describe the data quite well, whereas the other ones underestimate the normalization by as much as a factor of 2. However, the energy dependence of the cross section is well described by all models. Qualitatively similar conclusions were reached in the study of Ref. .

In order to understand this discrepancy between the different models, we plot in Fig. 4 the integrand $`p(r,Q)/r`$ which appears in Eq. (10) together with the dipole cross section from the Golec-Biernat and Wüsthoff model. One can see that different models for the wave functions result in different shapes in $`r`$ for $`p(r,Q)/r`$. It is clear from this plot that the exact shape of the function $`p(r,Q)/r`$ for large values of $`r`$ is very important in this case since it is weighted by the dipole cross section $`\sigma (x,r)`$ which is large in this regime.

### 2.3 Overlap integral $`N(Q)`$, mean radius $`r_Q`$ and theoretical uncertainty on $`S(x,r_Q,b)`$

The results on the forward production cross section presented in Fig. 3 suggest that it is essential to study in detail the model dependence of the $`S`$-matrix given by formulae (12) and (13).

Let us first study the properties of the quantity $`N(Q)`$. In Fig. 5 we plot it for the three different models for the vector meson wave function. In the region of interest, $`0.05<Q^2<10.0\mathrm{GeV}^2`$ they result in comparable $`N(Q)`$ within $`1025\%`$. This means that the average $`S`$-matrix obtained using Eq.(12) is quite robust to a given choice of the model for the vector meson wave function. In fact one can give some very general bounds on $`N(Q)`$, nearly independent on the model for $`\psi _V`$. From these upper bounds on $`N(Q)`$ we could obtain an absolute upper bound for the $`S`$-matrix which would be model-independent. We refer to appendix B for technical details. We represent these bounds in Fig. 5. One can see that the models happen to be quite close to these limits (within say $`25\%`$).

Second, we estimate the variation of $`r_Q`$ with respect to the different models for the wave functions. We compute the mean value $`r_Q`$ as:

$$r_Q=e^{\mathrm{log}r_p}.$$
(20)
The results for the different models at different values of $`Q^2`$ are given in Tab. 1. The values of $`r_Q`$ are consistent within $`25\%`$ for the various models. This confirms the validity of formula (13).

The three models appear to be more or less equivalent for our purpose. For simplicity, we then only consider model DGKP  in the following.

## 3 Impact parameter analysis of HERA data

In this section we extract the $`b`$-dependent $`S`$-matrix from the HERA data. To this aim, we apply formulae (12,13) to the data, together with expressions (17,19) obtained from the discussion of the vector meson wave function in the previous section. We also try to deduce from these results the value of the saturation scale.

### 3.1 Profile $`S(b)`$

The experimental data for diffractive production of vector mesons are usually parametrized by the forward diffractive cross section $`d\sigma /dt|_{t=0}`$ and the logarithmic slope in momentum transfer $`B(t)`$ as follows:

$$\frac{d\sigma }{dt}=\frac{d\sigma }{dt}_{|_{t=0}}e^{B(t)t}.$$
(21)
We take the available data  for the electroproduction of longitudinally polarized $`\rho _0`$ vector mesons. These data are given for momentum transfer $`t`$ below $`t_{\mathrm{max}}=0.6\text{GeV}^2`$, which enables us to determine reliably the $`S`$-matrix only for impact parameter values larger than $`b=1/\sqrt{t_{\mathrm{max}}}0.3\text{fm}`$. In order to compute the Fourier transform appearing in Eq. (12), we have to assume an extrapolation of $`d\sigma /dt`$ to larger values of $`t`$. The formula we use is

$$\begin{array}{c}S(x,r_Q,b)=1\frac{1}{2\sqrt{\pi }N(Q)}\sqrt{\frac{d\sigma }{dt}}_{|_{t=0}}(_0^{t_{\mathrm{max}}}dtJ_0(b\sqrt{t})e^{B(t)t/2}\hfill \\ \hfill +_{t_{\mathrm{max}}}^+\mathrm{}dtJ_0(b\sqrt{t})\sqrt{(t)}).\end{array}$$
(22)
Assuming that $`B(t)`$ is non-increasing with $`t`$, as various sets of data seem to indicate, we choose different functional forms for the extrapolating functions $`(t)`$:

* an exponential form $`\mathrm{exp}(dt)`$, with the constant $`d`$ being of the order of $`B(t_{\mathrm{max}})`$. This choice provides an upper bound on the $`S`$-matrix for $`b`$ close to 0.
* a power law form $`t^\alpha `$. This choice is motivated by the fact that the data for photoproduction at high $`t`$  indicate a $`t`$-dependence governed by such a power law, with $`\alpha =3`$. Such less steep $`t`$-dependence was also obtained in a theoretical calculation . We consider here two values for $`\alpha `$: $`\alpha =3`$ and $`\alpha =6`$.

The parameters of these functions $`(t)`$ are fixed from the fit to the experimental points corresponding to the highest values for $`t`$, in the region $`0.4\text{GeV}^2<t<0.6\text{GeV}^2`$. These specific assumptions give an idea of the uncertainty on the determination of $`S(x,r_Q,b)`$ due to our lack of knowledge of the differential cross section for $`t>t_{\mathrm{max}}`$. Another source of uncertainty is the experimental errors on the measurements of $`B(t)`$ and $`d\sigma /dt|_{t=0}`$. A complete error analysis would lie beyond the scope of this paper. Here we just estimate the influence of these uncertainties on $`S(x,r_Q,b)`$ by varying the measured quantities inside the 1-$`\sigma `$ error bars. By this method, we believe that we obtain a strict overestimate of the errors. The theoretical input $`N(Q)`$ neeeded in formula (22) is computed within the model DGKP. The value $`r_Q`$ is estimated using the procedure described in Sec. (1) and (2). The other models NNPZ(1997) and NNZ(1994) have also been tested and lead to very similiar results.

We take three values of $`Q^2`$: $`Q_1^2=0.45\text{GeV}^2`$, $`Q_2^2=3.5\text{GeV}^2`$ and $`Q_3^2=7\text{GeV}^2`$. These values correspond to bins of the ZEUS analysis , which we consider in the following. We note that the H1 data  lead to similar results. We recall that these values of $`Q^2`$ correspond to the respective values of the dipole size $`r_Q`$ (see Tab. (1)): $`r_{Q}^{}{}_{1}{}^{}=0.35\text{fm}`$, $`r_{Q}^{}{}_{2}{}^{}=0.21\text{fm}`$, $`r_{Q}^{}{}_{3}{}^{}=0.16\text{fm}`$. For these values of $`Q^2`$, we take similar low values of $`x`$, namely $`x_1=4.710^4`$, $`x_2=4.310^4`$ and $`x_3=5.810^4`$ respectively. The experimental slope $`B(t)`$ is parametrized in the region $`t<0.6\text{GeV}^2`$ by $`B(t)=B_0ct`$, where the values of $`B_0`$ and $`c`$ are $`B_0=9.5\mathrm{GeV}^2`$, $`c=4.0\mathrm{GeV}^4(\mathrm{for}Q^2=0.45\mathrm{GeV}^2)`$ and $`c=6.1\mathrm{GeV}^4(\mathrm{for}Q^2=3.5\mathrm{GeV}^2\mathrm{and}Q^2=7\mathrm{GeV}^2`$.

The results for $`S(x,r_Q,b)`$ are shown in Fig. (6), and are obtained using Eq. (22) applied to the ZEUS data. For each $`Q^2`$, the curves corresponding to three choices for the extrapolation function $`(t)`$ are drawn. The uncertainty corresponding to the errors on the measurements themselves are computed for the intermediate value of $`Q^2`$, and are represented in the lower part of the plot by a hashed band.

First, we observe that for $`b>0.3\text{fm}`$, the curves corresponding to different values of $`Q^2`$ are clearly separated and ordered according to:

$$S(x_1,r_{Q_1},b)<S(x_2,r_{Q_2},b)<S(x_3,r_{Q_3},b).$$
(23)
As $`x_1x_2x_3`$, this means that the proton is less transparent to dipoles of larger size $`r`$. We also note that in this region of $`b>0.3\text{fm}`$, the results are not very dependent on the choice of extrapolation of the $`t`$-dependence of the data. However, for $`b<0.3\text{fm}`$ (shaded region in Fig. (6)), the results become very sensitive to the asymptotic behaviour of the cross section $`d\sigma /dt`$, and data at larger $`t`$ are needed to be able to make firm statements about the value of $`S(x,r_Q,b)`$ near $`b=0`$. We also note that the errors on the measurements for $`t<0.6\text{GeV}^2`$ lead at most to an uncertainty of about $`25\%`$ for the $`S`$-matrix at impact parameter $`b=0.3\mathrm{fm}`$.

Finally, one can compute the total cross section for dipole-proton scattering using the following formula (which is nothing more than Eq. (6) averaged over $`r`$, for real $`S`$-matrix):

$$\sigma _{\text{tot}}^{q\overline{q}p}(x,r)_p=2d^2𝐛(1S(x,r,b)_p).$$
(24)
Using formula (12) to replace $`S(x,r,b)_p`$ and performing the integral over $`𝐛`$, one sees that this cross section is determined by $`N(Q)`$ and by the forward differential cross section $`d\sigma /dt_{|t=0}`$. For photon virtualities $`Q_1^2`$, $`Q_2^2`$ and $`Q_3^2`$, one obtains the values $`\sigma _{\text{tot}}^{q\overline{q}p}_p=14.9\mathrm{mb}`$, $`10.6\mathrm{mb}`$ and $`7.5\mathrm{mb}`$ respectively. These results are consistent with those of Ref. , however, the relationship between $`Q^2`$ and the radius $`r_Q`$ quoted in these papers remains a large source of uncertainty in this comparison.

### 3.2 Towards an estimate of the saturation scale

One can try to estimate the quark saturation scale $`Q_s^2`$ using the results presented in Fig. 6. From this analysis it is in principle possible to extract the dependence of this scale on the impact parameter $`b`$ of the collision. We try the following phenomenological formula  for the $`S`$-matrix:

$$S(x,r_Q,b)=\mathrm{exp}(r_Q^2Q_s^2(x,b)/4).$$
(25)
This exponential form comes from a Glauber-like summation of multiple independent scatterings of the quark-antiquark pair on the target nucleon . It is also supported by the Golec-Biernat and Wüsthoff model  where the radius $`R_0(x)`$ introduced there can be related to the saturation scale by $`R_0(x)=1/Q_s(x,b)`$. However, in that work, it was assumed that the $`b`$-profile of the $`S`$-matrix has the form of a sharp cutoff at a distance of order 1 fm.

Using the approximation $`r_Q2/\sqrt{Q^2+m_V^2}`$ one obtains $`Q_s^211.5\mathrm{GeV}^2`$ and $`Q_s^20.2\mathrm{GeV}^2`$ for $`b=0.3\mathrm{fm}`$ and $`b=1.0\mathrm{fm}`$ respectively. These values are consistent in order of magnitude when computed at different $`Q^2`$, which confirms the relevance of formula (25). These results would suggest that at small impact parameter values the saturation scale is in the semi-hard regime and support the onset of the shadowing effects at HERA. One should however stress that these estimates are rather rough and should be taken with care since they strongly depend on the precise value of $`r_Q`$. Indeed, $`r_Q`$ comes squared in the formula for $`Q_s^2`$.

## 4 Summary

In this paper we have shown that using the HERA data on diffractive electroproduction of vector mesons it is possible to extract the $`S`$-matrix element $`S(x,r,b)`$ for dipole-proton scattering in impact parameter space. By considering the full $`t`$ dependence of this process we have shown how to obtain the $`S`$-matrix element averaged over dipole sizes $`r`$. Due to the particular properties of the overlap function between photon and meson wave functions this average can be replaced by the $`S`$-matrix element evaluated at $`r=r_Q2/\sqrt{Q^2+m_V^2}`$. We have then shown that our results on the $`S`$-matrix are only weakly dependent on the choice for the model of the meson wave function. Since the data are available only for low momentum transfer, $`t<0.6\mathrm{GeV}^2`$, we have shown that our results are reliable for impact parameter values larger than $`b0.3\mathrm{fm}`$. It would be very interesting to extend the measurement of diffractive electroproduction of vector mesons to larger values of $`t`$. This would enable us to explore the very interesting regime of central impact parameter collisions at HERA. We have also estimated the value of the dipole cross section integrated over the impact parameter and have found it to be consistent with other analyses. Finally we have discussed how to use this result to estimate the saturation scale at HERA collider. We have suggested that the onset of the shadowing effects can be dependent on the impact parameter of the collision and that for central collisions it could occur in the semihard regime. However further theoretical and experimental studies are necessary.

## Acknowledgements

We thank Marcello Ciafaloni for his comments, and Allen C. Caldwell and Mara S. Soares for correspondance. S.M. thanks Columbia University for welcome at a preliminary stage of this work, and the Service de Physique Théorique de Saclay for support at that time. He also thanks Robi Peschanski and Bernard Pire for their suggestions. A.M.S. thanks Krzysztof Golec-Biernat and Jan Kwieciński for interesting discussions. A.H.M. wishes to thank Marcello Ciafaloni for his hospitality at the University of Florence where this collaboration began, and he wishes to thank Dominique Schiff for her hospitality in Orsay.

## Appendix A Models for the meson wave function

In this appendix, we detail the phenomenological forms we adopt for the meson wave function.

### A.1 Model by Dosch, Gousset, Kulzinger, Pirner (DGKP, Ref. )

The form of the vector meson wave function is given by:

$$\psi _V^{h,\overline{h}}(z,r)=\delta _{h,\overline{h}}\frac{𝒩}{\sqrt{4\pi }}4\omega (z(1z))^{3/2}\mathrm{exp}\left(\frac{1}{2}\frac{m_V^2}{\omega ^2}(z1/2)^2\right)\mathrm{exp}\left(\frac{1}{2}\omega ^2r^2\right).$$
(26)
The parameters chosen are the following:

$$\omega =0.330\text{GeV}\text{and}𝒩=4.48.$$
(27)
An additional feature of this model is that the mass $`m_q`$ of the quarks that appears in the photon wave function is running with $`Q^2`$. This gives important effects only for very low $`Q^2`$ and for photoproduction. In this regime, it is argued that the quarks should have constituent mass. The formula for the running mass reads:

$$m_q(Q^2)=\{\begin{array}{cc}\hfill 0.220\text{ GeV}(1Q^2/Q_0^2)& \text{for}Q^2<Q_0^2=1.05\text{ GeV}^2\hfill \\ \hfill 0& \text{for}Q^2>Q_0^2\hfill \end{array}$$
(28)
This model is referred to as DGKP.

### A.2 Models by Nemchik, Nikolaev, Predazzi, Zakharov (NNZ(1994), Ref.  and NNPZ (1997), Ref. )

In this model, the radial wave function $`\varphi (r,z)`$ appearing in eq.(19) satisfies the following normalization condition:

$$1=\frac{N_c}{2\pi }_0^1\frac{dz}{z^2(1z)^2}d^2𝐫\left\{m_q^2\varphi ^2(z,\mathrm{r})+(z^2+(1z)^2)(_r\varphi (z,\mathrm{r}))^2\right\}.$$
(29)
This equation is nothing else but the normalization condition (14), applied to the transversely polarized meson wave function. The obtained normalization factor $`\mathrm{\Psi }_0(1S)`$ is assumed to be the same for a longitudinally polarized meson. The function $`\varphi (r,z)`$ is defined by:

$$\begin{array}{c}\varphi (z,r)=\mathrm{\Psi }_0(1S)\hfill \\ \hfill \times \{4z(1z)\sqrt{2\pi R^2}\mathrm{exp}(\frac{m_q^2R^2}{8z(1z)})\mathrm{exp}(\frac{2z(1z)r^2}{R^2})\mathrm{exp}\left(\frac{m_q^2R^2}{2}\right)\\ \hfill +C^4\frac{16a^3(r)}{A(z,r)B(z,r)^3}rK_1(rA(z,r)/B(z,r))\}.\end{array}$$
(30)
The functions $`A`$ and $`B`$ are given by

$$\begin{array}{cc}\hfill A^2(z,r)& =1+\frac{C^2a^2(r)m_q^2}{z(1z)}4C^2a^2(r)m_q^2\hfill \\ \hfill B^2(z,r)& =\frac{C^2a^2(r)}{z(1z)},\hfill \end{array}$$
(31)
and $`a(r)3/(8m_q\alpha _s(r))`$. The following prescription is taken for the running coupling $`\alpha _s(r)`$:

$$\alpha _s(r)=\alpha _0\text{for}r>r_s\text{and}\alpha _s(r)=\frac{4\pi }{\beta _0\mathrm{log}(1/\mathrm{\Lambda }^2r^2)}\text{for}r<r_s,$$
(32)
where the parameter values are

$$r_s=0.42\text{fm},\alpha _0=0.8,\mathrm{\Lambda }=200\text{MeV}.$$
(33)
The masses of the quarks are taken to be $`0.15\text{GeV}`$.

The other parameters $`\mathrm{\Psi }_0(1S)`$, $`C`$ and $`R^2`$ are chosen taking into account several constraints: the normalization condition (29) for the wave function has to be satisfied, and the value of the leptonic decay width must agree with the experimental measurement. Additionally, the mean radius of the meson has to be of the order of a hadronic scale. Two different sets of parameters are considered:

$$\begin{array}{cc}\hfill C& =0.25,R^2=0.76\text{fm}^2\text{(reference }\text{[9]}\text{)}\hfill \\ \hfill \text{and}C& =0.36,R^2=1.37\text{fm}^2\text{(reference }\text{[10]}\text{)}.\hfill \end{array}$$
(34)
These two choices are refered to as NNZ(1994) and NNPZ(1997)respectively.

## Appendix B Model-independent bounds on the overlap integral $`N(Q)`$

In this appendix we explore more the properties of $`N(Q)`$ and show that upper bounds can be given on this quantity, regardless the model adopted for $`\psi _V`$. These provide upper bounds on $`S(x,r_Q,b)`$.

$`N(Q)`$ can be seen as a scalar product of the two wave functions, $`N(Q)=\psi _\gamma ^{}|\psi _V`$, see Eqs. (9) and (10). The Cauchy-Schwartz inequality then applies. It leads to an upper estimate for the integral: the product of the two wave functions is smaller than the square root of the product of the integrals of the squared wave functions. Using the normalization condition (14) for the meson wave function444 Strictly speaking, this bound is only valid for models for which the condition (14) is enforced for the longitudinally polarized vector meson.  and computing the full integral of the squared photon wave function, one eventually obtains:

$$N(Q)2\widehat{e}_V\frac{\alpha _{\text{em}}N_c}{6\pi }\left(16\frac{m_q^2}{Q^2}+24\frac{m_q^4}{Q^4}\frac{1}{\sqrt{1+\frac{4m_q^2}{Q^2}}}\mathrm{tanh}^1\frac{1}{\sqrt{1+\frac{4m_q^2}{Q^2}}}\right).$$
(35)
The r.h.s. of eq.(35) grows with $`Q^2`$ and thus this inequality is only interesting in the small-$`Q^2`$ region. Although independent of the precise form of the vector meson wave function, this bound nevertheless depends on the masses of the quarks.

Second, with the additional assumption that $`\psi _V(z,𝐫)`$ is maximum for $`r=0`$, we can write:

$$N(Q)2\widehat{e}_Vd^2𝐫𝑑z\psi _\gamma ^{}(z=1/2,𝐫;Q)\psi _V(z,𝐫=0).$$
(36)
The integrals over $`z`$ and $`𝐫`$ in the r.h.s. are now factorized. The integration over $`𝐫`$ is performed analytically, while the one over $`z`$, involving $`\psi _V`$, can be expressed as a function of the coupling $`f_V`$ of the meson to the electromagnetic current, using the relation (15). This finally leads to:

$$N(Q)2\sqrt{\pi \alpha _{\text{em}}}f_V\frac{Q}{Q^2+4m_q^2}.$$
(37)
The r.h.s. vanishes like $`1/Q`$ at large Q, which makes this inequality useful for large $`Q`$.

The region in the ($`Q^2`$,$`N(Q)`$) plane which is forbidden by these bounds is depicted in Fig. (5). One sees that all the models are very close to the upper bounds.