text:
# The (-1) twist

## Appendix

The purpose of this appendix is to give an alternate and somewhat shorter proof of the main Lemma (Lemma 1.2) of \[CLR\], which goes as follows.

###### Lemma

Let $`Y()`$ be an irreducible subvariety spanning a $`^{p+1}.`$ Fix integers $`r,k`$ with $`0r1kd`$. Let $`_d^N`$ be a linear map of vector spaces such that for $`y_1,\mathrm{},y_{r1}`$ general points of $`Y`$, the restriction to $`_d(y_1\mathrm{}y_{r1})`$ induces a surjection

$$\mathrm{\Psi }:_d(y_1\mathrm{}y_{r1})^{k+1}.$$

Then for a general choice of elements $`h_{k+1},\mathrm{},h_d`$ and for general subsets $`Y_1,\mathrm{},Y_kY`$ each of cardinality $`p`$, with $`Y_iy_i`$, $`i=1,\mathrm{},r1`$, the restriction of $`\mathrm{\Psi }`$ to the subspace $`(Y_1)\mathrm{}(Y_k)h_{k+1}\mathrm{}h_d`$ surjects.

###### Demonstration proof

By a good chain in $``$ we mean a (connected) chain whose components are straight lines (i.e. pencils) of the form $`(S)`$ where $`S`$ is a general $`p`$-tuple in $`Y`$ and whose ’vertices’ (i.e. singular points) are general in $`()`$. Clearly two general elements of $``$ can be joined by a good chain. Let $`_d_d`$ be the set of monomials. By a good chain in $`_d`$ we mean a chain which is a union of subchains of the form $`𝒞_i^{}=h_1\mathrm{}𝒞_i\mathrm{}h_d`$ where $`𝒞_i`$ is a good chain in $``$. It is easy to see that two general monomials can be joined by a good chain.

Next, let us say that a monomial $`h_{de+1}\mathrm{}h_d_e`$ is rel $`\overline{y},\overline{y}=\{y_1,\mathrm{},y_{r1}\}`$ if $`h_i(y_i),i=de+1,\mathrm{},r1`$ (this condition is vacuous if $`r1<de+1`$); denote by $`_e(\overline{y})`$ the set of these. Again it is easy to see that two general elements of $`_e(\overline{y})`$ can be joined by a good chain within $`_e(\overline{y})`$.

Now to prove the Lemma it suffices to show by induction on $`q,0qk`$, that, with the above notations,

$$dim\mathrm{\Psi }((Y_1)\mathrm{}(Y_q)h_{q+1}\mathrm{}h_d)\mathrm{min}(q+1,k+1)$$

(for general choices rel $`\overline{y}`$). For $`q=0`$ this is clear as $`\mathrm{\Psi }`$ is nonzero on a general element of $`_d(\overline{y})`$, because these span $`_d(\overline{y}`$). Assume it is true for $`q`$ and false for $`q+1`$, and suppose first that $`q+1r1`$. Now because

$$Z:=\mathrm{\Psi }((Y_1)\mathrm{}(Y_q)h_{q+1}\mathrm{}h_d)=\mathrm{\Psi }((Y_1)\mathrm{}(Y_q)(Y_{q+1})h_{q+2}\mathrm{}h_d)$$

by assumption (i.e $`Z`$ doesn’t move as $`h_{q+1}`$ varies in $`(Y_{q+1})`$, and because two general elements of $`(y_{q+1})`$ can be joined by a good chain (whose components are of the form $`(Y_{q+1})`$), it follows that $`Z`$ is independent of $`h_{q+1}(y_{q+1})`$, fixing the other $`h`$’s. Since $`y_1,\mathrm{}y_{r1}`$ are all interchangeable, it follows that a similar statement holds for any permutation of them. In particular $`Z`$ contains all pencils of the form $`\mathrm{\Psi }(h_1\mathrm{}(Y_j)\mathrm{}h_d),j=1,\mathrm{},r1`$ and from connectivity of $`_{r1}(\overline{y})`$ by good chains it follows that in fact $`Z=\mathrm{\Psi }(_{r1}(\overline{y})h_r\mathrm{}h_d)`$. As for the remaining $`h_j`$’s, say $`j=r`$, we may pick $`y_rf^{}(h_r)_0`$ and apply similar reasoning to $`y_1,\mathrm{},y_{r1},y_r`$ in place of $`y_1,\mathrm{},y_{r1}`$ . The foregoing argument yields that $`Z`$ is independent of $`h_r(y_r)`$, fixing $`y_r`$ and the other $`h`$’s (indeed that $`Z=\mathrm{\Psi }(L_r(\overline{y}y_r)h_{r+1}\mathrm{}h_d)`$. Then another similar argument with good chains of $`h_j`$’s not fixing $`y_j`$ yields easily that $`Z`$ is actually independent of $`h_j`$. Now since we can connect two general elements of $`_d(\overline{y})`$ by a good chain, we conclude that $`Z=\mathrm{\Psi }(_d(\overline{y}))`$, which is a contradiction.

The case $`q+1r`$ is similar but simpler: we may conclude directly that $`Z`$ is independent of $`h_{q+1}(Y_{q+1})`$, hence of $`h_j(Y_j)`$ for all $`rjk`$ and use good chains as above to deduce a contradiction.∎


format: text

text:
# Velocity Structure of the ISM as Seen by the Spectral Correlation Function

## 1 Introduction

In the past few decades, the turbulent character of the velocity field in the interstellar medium has been increasingly recognized (Zuckerman & Evans 1974; Dickman 1985; Scalo 1987; for recent works, see the volume by Franco & Carramiñana 1999 and the review by Vázquez-Semadeni et al. 2000a). Due to the highly fluctuating character and sensitivity to initial conditions (initial conditions arbitrarily close to each other separate exponentially in time – see, e.g., Lesieur 1990) of the turbulent motions, statistical tools become important in trying to characterize and understand the nature of interstellar turbulence. These tools become even more important with the continuous growth of both observed and theoretical data cubes, as it becomes necessary to compress such large data sets into a more manageable parameter space. On the observational side, studies of the centroid velocity probability distribution function ($`v`$-PDF) (Miesch & Scalo 1995; Lis et al. 1996, Miesch, Scalo & Bally 1999) have compared the shape of the $`v`$-PDF for various observed regions and numerical simulations of incompressible or weakly compressible turbulence; through the density and velocity autocorrelation functions, Kleiner & Dickman (1984) have tried to infer the correlation length scale in Taurus Molecular Cloud; using structure tree statistics, Houlahan & Scalo (1992) tried to discriminate between hierarchically-nested and random collections of clouds; through Principal Component Analysis and autocorrelation functions, Heyer & Schloerb (1997) and Brunt & Heyer (2001) obtain the velocity dispersion-size relationship for both observations and pseudo-simulations, in a manner independent of how clouds are defined; applying $`\mathrm{\Delta }`$variance analysis, Stutzki et al. (1998) and Bensch et al. (1999) have measured the fractal dimension of the projected intensity towards the Polaris Flare. Finally, Lazarian & Pogosyan (2000) have given a method for recovering the power spectra of both the three-dimensional density and of the line-of-sight velocity from the emissivity power spectrum in velocity channels of spectroscopic data.

On the numerical side, recent simulations of interstellar turbulence have allowed statistical studies of the simulated ISM and clouds. In fact, again due to the sensitivity to initial conditions, modeling of individual turbulent objects is not feasible, and statistical modeling is again the appropriate way to proceed. It thus becomes appropriate to compare numerical simulations and observations using statistical methods. For example, Falgarone et al. (1994) studied optically thin line-profiles in 5123 simulations, concluding that the shapes of the spectra, their moments, and their spatial variation exhibit similarities with observational line-profiles of non-star forming regions; Vázquez-Semadeni, Ballesteros-Paredes & Rodríguez (1997) searched for Larson-type relationships in numerical simulations, finding in particular that the density-size relationship is not satisfied, but instead there exists a whole family of low-column density clouds that may however be missed by observational surveys limited by integration time (see also Scalo 1990); Padoan & Nordlund (1999) compare sub- and super-Alfvénic simulations with observations, concluding that the first case may be in conflict with observations; Pichardo et al. (2000) studied the projection of numerical MHD simulations, finding that the morphology of the channel maps resembles maps of the line of sight (LOS) of the velocity more closely than maps of the density field; Mac Low & Ossenkopf (1999) have recently used the $`\mathrm{\Delta }`$-variance to find that the models with high Mach numbers (more than M $``$ 4) are in better agreement with observations; Lazarian et al. (2001) have shown that the theoretical predictions of Lazarian & Pogosyan (2000) are verified for simulated spectroscopic data from numerical simulations of ISM turbulence.

In particular, Rosolowsky et al. (1999 = RGWW) introduced a new method, called the “Spectral Correlation function” (SCF), which estimates the spatially-averaged111In this paper, the word “spatial” refers to the plane of the sky, unless otherwise stated. correlation of spectra inside a box of size $`s`$ in spectroscopic data cubes. These authors used the SCF to measure the correlation between spectra in neighboring positions, as an indicator of the degree of “smoothness” of both observational and numerical data cubes. With this method, they were able to distinguish between, for example, simulations with and without self-gravity, by comparing the variation of the SCF upon randomization of the positions of the spectra on the plane of the sky (POS).

In the present work, we extend the work of RGWW, further exploring the capabilities of the SCF to discriminate among numerical simulations of various turbulent regimes (isothermal with large- and small-scale forcing, and ISM-like), and to quantitatively compare them to spectroscopic observational data of diffuse HI gas. In order to do this, we construct spectral maps from the numerical data cubes assuming optically thin lines, with and without thermal broadening of lines. For non-isothermal gas, the velocity structure of the colder gas may be “hidden” by the thermal broadening of the warmer gas.

The plan of the paper is as follows: in §2 we describe both the observed HI data (Hartmann & Burton 1997 ) and the simulated data cubes used here, and explain the construction of the spectral maps for the simulations. In §3 we briefly review the SCF. In §4 we apply the SCF to the data sets, discussing the information that can be extracted from this procedure towards distinguishing between the various simulations and quantitatively estimating their similarity with observational data. In §5 we discuss the nature of the spatial structure of the SCF, and discuss how the SCF can be used to guide theoreticians to match their simulations with observational data. Finally, in §6 we draw the main conclusions.

## 2 Data

### 2.1 Numerical Simulations

We consider three numerical simulations in three dimensions. The first one, called simply ISM, represents the ISM in a box of 300 pc on a side, at a resolution of 100 points per dimension, centered on the Galactic midplane at the solar Galactocentric distance. Full equations, details and a complete description of the model can be found in Passot et al. (1995 =PVP) and Pichardo et al. (2000). Here we just mention that the code solves the full self-gravitating MHD equations, including that of the internal energy conservation with additional model terms representing the background plus local stellar heating and the radiative cooling. The local heating mimics stellar ionization heating, and is turned on (“an O star is formed”) when the local density reaches a threshold value of 10 times the mean density, producing local bubbles of warm gas that expand and feed the global turbulence. “Stars” remain on for 2 Myr. No supernova-like energy injection is included because of numerical limitations. The cooling is parameterized as a piecewise power-law function, as fitted to the cooling functions of Dalgarno & McCray (1972) and of Raymond, Cox & Smith (1976) by Rosen, Bregman & Norman (1993) and Rosen & Bregman (1995) (see PVP and Pichardo et al. 2001). The “isothermal” Jeans length is twice the integration box size, although the presence of cooling gives an effective Jeans length of roughly 2/3 the box size (Ballesteros-Paredes, Vázquez-Semadeni, & Scalo 1999). The mean density is 3.3 cm$`3`$ and the velocity, temperature and magnetic field units are, respectively, 6.4 km s$`1`$, 3000 K and 5 $`\mu `$G.

Two limitations of the ISM simulation should be mentioned. First of all, the lack of SN energy input implies that the simulations are somewhat less energetic than the actual ISM. This fact will be quantified in §5.2. On the other hand, the heating and cooling functions used do not produce any thermally unstable temperature ranges, although we do not consider this to be a serious problem, as recent discussions by Vázquez-Semadeni et al. (2000b) and Gazol et al. (2001) suggest that the structural properties of the ISM are more a consequence of the dynamical processes than of the thermal ones. On the other hand, a true shortcoming is the absence of hot ($`T10^6`$ K) gas, although recent work by various groups (e.g., Gazol-Patiño & Passot 1999; de Avillez 2000) suggests that the filling factor of the this gas at the midplane is not larger than $`20`$%. In this case, the error commited by not including it is probably minor, as the main effect of this omission will be that in actual observations a fraction of the line of sight will contain “holes” in the cold/warm gas distribution, filled with hot ionized gas, which will be absent in our simulation data. Given the low resolution of our simulations, this may actually constitute an advantage, as better sampling of the gas contributing to the HI line spectra along each LOS will be accomplished.

The second simulation, called ISM-IT, is similar to run ISM, except that an isothermal regime is considered, with no cooling and heating terms, the isothermal Jeans length is 1.1 times the box size, and a lower density threshold is used for initiating star formation, so that the star formation rate is similar to that of run ISM in spite of the weaker gravity. In this run, which may be thought to represent gas within a molecular cloud at $`T`$ 10 K, the stellar energy input proceeds via winds (locally divergent source terms in the momentum equation – see Vázquez-Semadeni et al. 1996; Avila-Reese & Vázquez-Semadeni 2001), and the internal energy evolution equation is bypassed altogether, but self-gravity and the magnetic field are still included. Finally, the third run, also at a resolution of $`100^3`$ and called IT, is a purely hydrodynamic isothermal run without self-gravity and with random compressible forcing applied at scales of 1/4 the simulation box size. As a summary, we present in Table 1 the main features of each one of the runs.

All three runs are started with Gaussian fluctuations with random phases, uncorrelated among all variables. Note that for run IT, the forcing is active since the beginning of the simulation, while for runs ISM and ISM-IT, the forcing is turned on only after the conditions for star formation mentioned above are met.

#### 2.1.1 Observing the Simulations

In order to produce “spectra” from the simulations, we generate two types of line profiles from the simulations. In one, we simply calculate density-weighted $`z`$-velocity histograms (position-position-velocity cubes, or PPV). This is equivalent to assuming that the emission is optically thin, but also that the thermal broadening is negligible, so that what we see is directly the emission from each pixel at its physical velocity along the LOS. Such spectra are produced by considering all points $`(x_0,y_0,z)`$ over the $`z`$-coordinate, and then creating a velocity axis by adding the contribution from each $`z`$-position to a velocity bin corresponding to its LOS-velocity ($`u_z`$). We refer to these as not-thermally-broadened (NTB) spectra. These synthetic observations, although unsuitable for direct comparison with real observational data because the latter can not escape some level of thermal broadening, allow us to study the effects of velocity structure exclusively.

For the second type of spectra, we calculate the radiative transfer taking the density field in the simulation to correspond to the HI density, and assuming local thermodynamic equilibrium, optically thin lines, and Gaussian thermal broadening. In this case, every pixel “emits” at all velocities, but its emission is weighted by a Gaussian profile centered at the fluid velocity of the pixel, and of width determined by the temperature at the pixel. We will refer to this kind of spectra as thermally-broadened (TB) spectra. In both cases, “channel maps” are then produced as intensity images over the remaining two spatial coordinates ($`x`$ and $`y`$) at a given $`u_z`$ interval.

In Fig. 1 we show the corresponding TB and NTB spectra for nine pixels around position $`(x=3,y=3)`$ of run ISM. The dotted lines represent the TB spectra, while the bold solid lines represent the NTB spectra at low (16 velocity channels) resolution, and the light solid lines represent the NTB spectra at high velocity resolution (64 channels). It is seen that the NTB spectra are much more irregular and spiky, because there are only 100 events in each LOS, which are then redistributed in either 64 (at high velocity-resolution) or 16 (at low velocity-resolution) possible bins. Note also that the spectra change substantially between the low- and high-resolution cases. In contrast, the TB spectra are much smoother, and they do not change so much with position on the POS, implying that thermal broadening “hides” the real velocity structure of the simulations. We will return to this point in §4.1.2.

It is worth pointing out that previous works constructing velocity histograms out of 3D numerical simulations have frequently added up the contribution from several neighbouring LOSs. For example, Falgarone et al. (1994), using $`512^3`$ simulations, considered data from a square box of 32 $`\times `$ 32 neighbouring LOSs, thus including 524288 events for each individual spectrum. However, at the risk of being velocity under-sampled, in the present work we use only one LOS per spectrum, in order to avoid further deterioration of the already low on-the-sky spatial resolution imposed by the machine available to us (a Cray Y-MP 4/64).

### 2.2 Observational Data

As the source of observational data, we use the 21 cm HI data from the Atlas of Galactic Neutral Hydrogen (Hartmann & Burton 1997) obtained with the 25 m Leiden/Dwingeloo telescope. Details of the observations and analysis may be found in the atlas. Here we just mention that the velocity resolution is 1.03 km s-1, and that we use the velocity range from $``$20 to $`+`$20 km s-1 (the original data was taken from $``$450 to 400 km s-1.)

We have chosen the North Celestial Pole Loop as the comparison region. This is an extended HI region located at $`120l160`$, $`10b50`$ (see § 4.3.) There are several reasons for this choice, mostly related to the fact that the ISM simulation represents a well-defined region of space, while HI observations suffer from confusion along the LOS. First, it is a region that lies off the Galactic plane, thus avoiding excessive contribution to the emission from very distant, physically unrelated regions. Secondly, the presence of a well-defined physical structure (the NCP Loop) reassures us that the emmision is mostly due to a well-confined region is space. Finally, there are molecular clouds and stars that allow to have an estimation of the distance, at least to some of the molecular clouds in the region (MBM 30, $`d=110\pm 10`$ pc), as well as an upper limit ($`d=346`$ pc) to other members (MBM 29; see Penprase 1993). This implies that, if the emission lies at a distance of $``$ 230 pc (the mean between these two extremes), the dimensions of the analyzed region has a size $``$160 pc, comparable to the physical size of the ISM simulation box.

Finally, it is worth noting that the resolution of the HI data changes with galactic latitude: it is 0.5 for $`b`$ 35; 0.6 for 35$`b<`$ 45; and 0.7 for 45$`b<`$ 55. This has the effect of slightly overestimating the SCF at high latitudes, since we assume a uniform grid spacing of 0.5.

## 3 The Spectral Correlation Function

A detailed discussion of the development and advantages of the SCF is given by RGWW. For completeness, here we repeat the necessary definitions, give a physical interpretation of the SCF, and then a discussion of the strategies for its application, and the effects of noise.

### 3.1 Definitions and Physical Interpretation

Given two arbitrary spectra with antenna temperatures $`T_1(v)`$ and $`T_0(v)`$, the SCF is defined by RGWW as:

$$\mathrm{SCF}1\sqrt{\frac{[sT_1(vl)T_0(v)]^2𝑑v}{s^2T_1^2(vl)𝑑v+T_0^2(v)𝑑v}}$$
(1)
where $`l`$ and $`s`$ are free parameters chosen as to minimize the integral $`\mathrm{{\rm Y}}=[sT_1(vl)T_0(v)]^2𝑑v`$. Thus, for two identical spectra with infinite signal to noise (s/n) ratio, SCF$`=1`$, while for two completely uncorrelated spectra, SCF = 0.

Given the parameters $`l`$ and $`s`$, RGWW define four modes of the SCF: a) SCF0 $``$ SCF$`(s=1,l=0)`$, which compares the two spectra in their raw form. Values of SCF0 close to unity imply that the spectra resemble each other in shape, amplitude and velocity offset. b) SCFl$``$ SCF$`(s=1,l)`$, with $`l`$ chosen as to minimize $`\mathrm{{\rm Y}}`$. This form eliminates differences between the spectra arising from pure velocity offsets. Values of SCFlclose to unity imply that the spectra resemble each other in shape and amplitude, regardless of velocity offset. c) SCFs$``$ SCF$`(s,l=0)`$, with $`s`$ chosen as to minimize $`\mathrm{{\rm Y}}`$. This form eliminates differences between the spectra arising from pure intensity scalings. Values of SCFsclose to unity imply that the spectra resemble each other in shape and velocity offset, regardless of absolute amplitude. d) SCFl,s$``$ SCF$`(s,l)`$, with both $`s`$ and $`l`$ chosen as to minimize $`\mathrm{{\rm Y}}`$. This is then the general definition of the SCF (eq. ), and eliminates differences arising from both velocity offset and amplitude. In this case, values close to unity imply that the two spectra resemble each other in shape only.

### 3.2 Effects of Instrumental Noise

As has been shown by RGWW (see also Padoan, Rosolowsky and Goodman 2001 = PRG), the SCF can be strongly affected by noise and their effects has to be evaluated. PRG propose to “clean” the effects of the noise by dividing the actual value of the SCF by an “optimal value” of the SCF, i.e., by the maximum value that the SCF would have if the only difference between both spectra were the random noise. The optimal value is given by (see PGR)

$$\mathrm{SCF}(Q)=1\frac{1}{Q},$$
(2)
where $`Q`$ is the “spectrum quality”, defined as

$$Q_i=\frac{1}{\sigma }\sqrt{\frac{\mathrm{\Sigma }_vT_i^2\delta v}{N_{\mathrm{ch}}\delta v}}=\frac{1}{\sigma }\sqrt{T_i^2},$$
(3)
where in turn $`\sigma `$ is the mean rms noise, $`N_{\mathrm{ch}}`$ is the number of velocity channels, $`\delta v`$ is their width, and $`T_i^2`$ is the mean value of the emission in the velocity interval. As can be seen, for high values of the signal-to-noise ratio, the parameter $`Q`$ goes to infinity, and the SCF is renormalized by a number close to unity.

PGR show that adding noise to both the observed (molecular, with signal-to-noise values of the order of a few) and the synthetic (without intrinsic noise) spectra until the signal-to-noise ratio is constant throughout the data, and renormalizing the SCF as indicated in eq. (2), gives values of the SCF essentially independent of the noise level.

In the present work we deal with data that are almost not affected by noise: on one hand, the numerical simulations are intrinsically noiseless; on the other, for the HI data (Hartmann & Burton 1997, see §2), the characteristic value of $`\sigma `$ is 0.07 K and its characteristic antena temperature value is of several K, suggesting high values of $`Q`$. Indeed, we have verified that for the observed data set (HI) used here, the mean value of $`Q`$ is 98.8 (see §4.3). Therefore, in this paper we do not correct for the effects of noise.

### 3.3 Velocity Window

In order to compute the four kinds of the SCF, in particular, the modes SCFl,s and SCFl, which imply shifting one spectrum with respect to the other, RGWW use velocity channels within 3 FWHMs of the velocity centroid, with the FWHM and the velocity centroid computed by fitting a Gaussian to the line profiles. In turn, PRG use a 10 $`\mathrm{\Delta }_2v`$ velocity window, where $`\mathrm{\Delta }_2v`$ is the standard deviation of the velocity of a spectrum created by averaging over the whole map. However, in our case of negligible noise, we simply extend the original data sets (observational and synthetic) by adding zeros on both sides of the spectral lines to produce a spectral window three times wider than the original.

### 3.4 Maps and Histograms of the SCF

The SCF, like the regular autocorrelation function, is in general a function of spatial separation, and this aspect of it is dicussed by PRG. However, RGWW has used it in a “local” mode, which we follow here. In this form, we construct maps of the average value of the SCF in boxes of 3 pixels per side in a given spectral map. That is, to each point on the $`(x,y)`$ plane (the POS of the simulations) we assign the average value of the SCF evaluated using that point and each of its eight nearest neighbors. Note that with this approach, the SCF maps are useful for exhibiting small-scale spatial variations in the line spectra. In addition, in order to quantify the similarity between spectral maps, RGWW also considered the distribution (histogram) of the SCF in a given map and its moments, and we follow this approach here.

In practice, the nearest-neighbor variant of the SCF used here acts like an angle-average of the spatial (i.e., in the POS) derivative of the spectral profiles. Sites of strong profile variation appear in the SCF maps as low values (dark), and soft variations appear as large values (bright). That is, the SCF stresses small-scale changes (on the POS) in the line spectra. Moreover, the fact that the SCF considers the whole velocity structure in the line profiles at each position, implies that maps of the SCF give more information than the channel maps alone, because they make velocity variations (as well as density ones) evident to the eye, which might escape examination of pairs of channel maps. The SCF histogram moments then quantify such velocity and intensity differences.

## 4 Results and Interpretation

### 4.1 Understanding the SCF

In order to “calibrate” the SCF, and to determine its response upon changes in the spectral resolution and thermal broadening, we have simulated four different “observations” of run ISM, labeled 16-TB, 16-NTB, 64-TB and 64-NTB, where the number indicates the number of velocity channels, and the letters TB/NTB indicate whether or not thermal broadening was included. We also have produced another set (TB and NTB) of “observations” of this run, but with the velocity artifitially multiplied by a factor of 6, in order for the line spectra and the SCF histograms to better match those for the NCP Loop. This will be discussed in §5.2. Velocity channel maps, maps of the SCF, and SCF histograms are respectively shown in figs. 2, 3 and 4. Each case (16/64 velocity channels, and TB/NTB) is indicated on the top of each figure.

#### 4.1.1 Effects of the velocity resolution

Comparing the NTB channel maps (Figs. 2b and d), we note that at higher velocity resolution the small-scale structure is more prominent, a fact already noted by Pichardo et al. (2000), and interpreted by Henney, Vázquez-Semadeni & Pichardo (2001) in terms of the presence of “caustics” or “cusps” (loci of high intensity due to velocity crowding) in the position-position-velocity (PPV) data cube. We will return to these caustics in §5.1. However, this effect is almost non-existent in the more realistic TB case.

Additionally, it can be seen that the NTB map at low velocity (or “spectral”) resolution (fig. 3b) exhibits a clear “wormy” or “stringy” structure, probably related to the caustics mentioned above (§5.2. However, the wormy structure is almost absent from the high-resolution SCF maps (fig 3d). This is an artifact of the low spatial resolution of our simulations (100 grid points per dimension), which causes any given point in the simulated POS to have very few points along the LOS contributing to the intensity in any one velocity channel (less than two on average at a velocity resolution of 64 channels), in turn causing the line profiles to be highly irregular (cf. fig. 1). Thus, the little signal available is diluted over many velocity channels, and only the most promintent features are preserved, with a very “blurred” quality. This shows that the effective velocity resolution is in effect limited by the spatial resolution. The scattered character of the NTB line profiles at high velocity resolution is also evident in their corresponding SCF histograms (fig. 4d), which are clearly much wider and have lower mean values than the low-resolution ones (fig. 4b), indicating that the profiles vary strongly among neighboring LOS.

#### 4.1.2 Effects of thermal broadening

We now turn to the effects of thermal broadening on velocity channel and SCF maps, and on the SCF histograms. For the discussion, we consider only the low-velocity resolution observations, due to the poor statistics present in the high-resolution ones.

It is first important to note that run ISM contains large regions occupied by warm ($`T1.5\times 10^3`$ K gas, whose thermal velocity dispersion dominates the line profiles at any given position on the POS, as seen in fig. 1. This is because, even though the low-temperature regions (clouds) contain supersonic motions, these motions are subsonic compared to the sound speed of the warm medium. Moreover, the warm gas is barely transonic with respect to itself. Thus, the thermal broadening due to the warm gas along the LOS masks the bulk velocity structure seen in the profiles, which are nearly perfectly Gaussian in the TB case, in contrast with the extremely irregular profiles seen in the NTB case. This effect is expected to occur in actual observations of multi-temperature atomic gas, so that the velocity structure of the cold HI gas may only be seen in absorption (Gibson et al. 2000; Heiles 2001).

This effect is also seen in the channel maps. The structure in the NTB channel maps, regardless of resolution, can be described as having the appearance of superposed veils, with sharp, bright edges delimiting fainter regions. In contrast, the TB channel maps are much smoother, losing the veil-like appearance, and becoming more similar to maps of the total integrated intensity. Thus, only the gross features such as large-scale velocity gradients remain noticeable. Again, this effect is also expected to occur in actual observational data of warm gas.

In the SCF maps, the “wormy” nature of the NTB cases, mentioned in §4.1.1, is seen to lose its small-scale part in the TB maps. That is, the “worms” in the NTB case appear to “bundle up” into larger-scale filaments, in a sort of self-similar or hierarchical manner. However, in the TB SCF maps, only the bundles appear to be preserved, while their smaller constitutive filaments are not seen. As will be discussed in §5, the small-scale filaments, if originated by the caustics in the PPV cube, would be an artifact of the velocity segregation of the data upon the process of a spectroscopic observation. Thus, thermal broadening has the fortunate effect of counteracting this spurious small-scale generation by the velocity separation effected by the observation process.

Another important feature is that the SCF maps for the TB observations have a much smaller dynamic range than the NTB ones, indicating a much lower variability of the SCF in the former case although, on the other hand, the same overall structures are seen in both types of maps. The SCF histograms clearly reflect the above effects. While the histograms of the NTB observations have large variances (fig. 4b) and moderate mean values, those for the TB (fig. 4a) cases are much narrower and with means very close to unity, the maximum possible value of the correlation. The small variances and large means of the TB cases again indicate the “erasure” of the velocity structure by the thermal broadening, which causes all spectra to be much more similar to each other than in the NTB cases. This furthermore indicates that the structure seen in the TB SCF maps should be due mostly to variations in intensity rather than in velocity structure. This is indeed confirmed by the fact that, in the TB cases, the histograms of SCFs, which can compensate for intensity variations but not for velocity shifts, are much tighter (have smaller variances) than those of SCFl, which compensates in the opposite manner. Instead, in the NTB cases, the histograms of SCFs have comparable (in the low velocity resolution case) or even larger (in the high-resolution case) variances than the histograms of SCFl, indicating that in this case the variability in the maps is due mainly to the velocity structure.

### 4.2 Comparison between numerical simulations

We now compare these results with two other numerical simulations. As mentioned in §2.1, we consider an isothermal simulation including the magnetic field and self-gravity (called ISM-IT) on one hand, and a simple isothermal simulation (called IT) on the other. Given the strong masking effect of thermal broadening discussed in the previous section, in this section we discuss only NTB cases, so that the velocity structure is fully accesible to the SCF. However, note that this cannot be done in the case of real observational data.

Figures 5a, and b show the velocity channels for cases ISM-IT and IT, respectively. The veil-like structure described for run ISM is also present in runs IT and ISM-IT, although in the latter the veils are more roundish, somewhat resembling soap bubbles, or peanut shells. This is most likely a reflection of the intermediate-scale forcing applied to run IT (in scales of about 1/4 of the box), in contrast with the small-scale forcing acting on runs ISM and ISM-IT, due to the stellar sources. The presence of the magnetic field in runs ISM and ISM-IT is probably also a factor contributing to the more elongated structures present there222Note that in the ISM simulations, $`\beta `$, the magnetic to thermal pressure ratio, is less meaningful than in isothermal simulations, since both the magnetic and thermal pressure are not constant and, in fact, $`\beta `$ is a local quantity, spanning a factor of at least 3 orders of magnitude..

Another evident distinction between the runs is that very-large-scale structures (runing diagonally from one side of the box to the other) of roughly half of the box size), which are present in run ISM (see e.g., velocity channel at $`4`$ km s-1 in fig 2b), are absent from cases IT and ISM-IT. This is because of the smaller effective Jeans length in run ISM (§2.1). Thus, such large-scale structures in run ISM have a gravitational origin, but have no counterpart in runs ISM-IT and IT.

The SCF maps of runs ISM-IT and IT (figs. 6a and b) continue to show the “wormy” appearance described for run ISM, although the “worms” form rounder bundles in run IT, due to the also rounder shapes of the structures seen in its channel maps. The SCF histograms for all three runs (figs. 3b and 6a and b) have comparable means and variances. However, the histograms for runs in which star formation is active (runs ISM and ISM-IT) seem to have more extended low-correlation tails, and therefore more negative values of the skewness. This appears to be a consequence of the greater abundance of shells and shocks in cases with active star-formation, and thus suggests that the histograms of the SCF are capable of “sensing” the presence of these structures, representing them by large negative skewness values. Unfortunately, the correspondence is not unique: the phenomena that can cause a large low-SCF tails are not limited to shocks and shells. Any sharp changes in the gas structure, either in density or in LOS-velocity, or even in the PPV cubes (e.g., the caustics), can cause them. Thus, the method cannot be used to unambiguously identify shocks and shells. Instead, the existence of low-SCF tails in the histograms should be taken only as suggestive of their possible presence, which should be then verified by other means.

### 4.3 SCF for the North Polar Loop

We now apply the SCF to the HI maps of the North Celestial Pole (NCP) Loop, a region possibly formed by a nearby supernova explosion (see, e.g., Pound & Goodman 1997). In fig. 8a we show the channel maps for this region, which contains an arc-shaped structure appearing at positive velocities, and a clear large-scale velocity gradient: at large negative velocities most of the emission is concentrated at high longitude and low latitude, while at large positive velocities most of the emission comes from the arc and is concentrated in the upper half of the region.

As we mentioned before (§2.2), the reason for choosing this region is that it is a high-latitude region, avoiding the contribution to the emission from physically disconnected regions far away in the line of sight, and that the determination of the distance to some molecular clouds allows us to assume that the size of the region is about 160 pc, comparable to the size of the simulations (300 pc).

Figure 8b then shows the SCF maps for this region. Several points are apparent. First, the “wormy” structure seen in all the simulations is also seen in the NCP Loop, particularly in the arc region, indicating the presence of shocks or caustics in the field (§§4.1.1 and 5)333The “worms” in the observations’ SCF maps at first sight appear larger than those in the simulations’ SCF maps. However, this is only due to the larger pixel size of the NCP Loop maps. The filaments are in all cases only a few pixels wide.. In this case, it is natural to assume that it is primarily shocks that are causing the “worms”, the NCP Loop being an expanding shell. Second, we see that the SCF maps only contain structure around the arc, and not near the bottom, where rather featureless emission is seen in the channel maps, indicating that the local form of the SCF used in this paper eliminates roughly uniform regions, and highlights regions of strong variability. Another consequence of this behavior is that the local form of the SCF eliminates large-scale velocity gradients, as indicated by the fact that no sign of the gradient seen in the channel map set is apparent in the SCF maps. Although this is in a sense obvious, it makes the local SCF maps a powerful tool for studying the small-scale intensity and velocity structure. In particular, the SCF should make a good shell-identification tool, in the sense that we can quantify, using the skewness of the histograms, the relative importance of sharp edges (shells) of gas.

Also worth noting is the fact that the map of SCFl contains significantly lower values of the correlation compared to the map of SCFs, indicating that most of the variability is due to the intensity rather than to the velocity structure. This is also seen in the SCF histogram for SCFl (fig. 8c), which shows a more extended low-correlation tail, and consequently a larger negative skewness than the SCFs histogram. However, these features may instead be an indication that thermal broadening is not negligible for this region, thus partially “blinding” the SCF to the velocity structure. If this is so, then the thermal broadening should be of 3-5 km s-1, since the total width of the HI line profiles for this region is about 8-10 km s-1, such that the total width comes from the turbulent structure, but the thermal broadening does not allow to see in detail such turbulent velocity structure. Indeed, it is noteworthy that the histograms of SCFl for the simulations have larger negative skewnesses than those of SCFs in the thermally broadened cases, but the opposite is true for the NTB cases, supporting this view. This is also probably the reason why the featureless emission at the bottom of the channel maps does not contain “worms” associated to the small-scale velocity structure in the SCF maps.

## 5 Discussion

### 5.1 “Worms” as shocks or caustics in PPV space

The local form of the SCF as used in this paper is by definition sensitive only to small-scale features (either in intensity or in velocity). Specifically, structural features in the SCF maps correspond to steep gradients in the plane of the sky of the line profiles’ shape and/or amplitude. In particular, shocks and the caustics discussed by Henney et al. (2001) are especially well suited for detection by the local SCF.

For clarity, a brief summary of the results of Henney et al. (2001) is in order here. Essentially, the process of performing a spectroscopic observation amounts to performing a mapping from the real, three-dimensional physical space onto position-position-velocity (PPV) space (see also Lazarian & Pogosyan 2000; Pichardo et al. 2000). As is well known, “caustics” form in PPV space at places where the mapping is not single-valued. i.e., where finite-length line-of-sight segments contribute to infinitesimal velocity intervals, causing a diverging intensity at those sites in the PPV cube (a phenomenon known as “velocity crowding”). With finite-width velocity channels, the intensity does not diverge, but still reaches large values at the caustics. The latter are surfaces in the PPV cube, and in general intersect the velocity channels (constant-velocity planes in the PPV cube) at finite angles along lines which show up as elongated, “worm”-like features in the channel maps. This structure is spurious in the sense that it does not correspond to any actual density features in the original physical space. In the case of sharply-defined velocity channels, the spurious structure reaches arbitrarily small scales as the velocity resolution is increased. Upon thermal broadening, the minimum scale present is explicitly calculated by Henney et al. (2001).

This spurious small-scale structure may however be picked out by the local SCF, showing up as the “worms” we have mentioned repeatedly. However, if the thermal broadening is very strong (i.e., the velocities are strongly subsonic), the caustics become thicker and less intense in PPV space, so that the “worms” become wider and fainter. Thus, as already mentioned in §4.1.2, thermal broadening present in real spectroscopic observations has the fortunate effect of counteracting the spurious small-scale generated by the synthetic observation process, although at the cost of blurring the velocity information.

The other possible type of small-scale feature that is naturally picked out by the local SCF is shocked shells. Contrary to the caustics, shocks are real features in 3D space, and, since they imply supersonic velocity differences occurring over very small scales, it is possible that they “survive” thermal broadening. We conclude that, when thermal broadening is not negligible, “worms” in SCF maps are most likely to correspond to real velocity features – shocks. However, for extremely supersonic motions where thermal broadening is negligible, it is possible that the SCF may not be able to distinguish between them and the spurious caustics. Further work is necessary to quantify and distinguish caustics and shocks through the SCF.

### 5.2 “Diagnosing” the simulations through the SCF

The SCF is intended to provide a quantitative means of comparing numerical simulations and observations, which should then allow workers to adjust their simulations to make them match the observational data as closely as possible. In this section we show an example of how this can be done.

Comparing both the NTB and the TB 16-channel SCF histograms of the ISM simulation (figs. 4a and b) with those of the NCP Loop (fig. 8c), we note that neither set of the simulation histograms is a very good match of those of the NCP Loop, and instead seem to bracket the latter. The NTB histograms are too broad and have smaller means than the NCP Loop ones, while the opposite is true of the TB histograms. In particular, the narrowness of the TB histograms suggests that the thermal width dominates the line spectra of the simulation, causing them to be very similar, all being very close to Gaussians. Indeed, as seen in fig. 1, the thermally broadened profiles have lost essentially all the velocity information and are strongly dominated by the thermal component. Instead, the 16-channel NTB profiles are very spiky and do not look like realistic profiles.

As mentioned in §4.1.2, we attribute this problem to the fact that, even though the cold gas contains supersonic motions, these motions are not supersonic compared to the thermal speed of the warm gas, and thus the velocity information of the cold gas is swamped in the thermal width of the warm gas. This is a real effect, and we expect it to affect the NCP Loop data. However, it moreover occurs that the latter include the contribution from gas moving at significantly larger velocities, most likely due to large shell expansion velocities (see Pound & Goodman 1997). Instead, the simulation only contains motions in the interval $`6`$ km s$`{}_{}{}^{1}v6`$ km s-1, as it includes only HII region-like expanding shells, but not supernova remnant-like ones, which would have much larger velocities. Thus, in the HI data there exists quite a larger velocity range, which “escapes” the thermal smoothing. This suggests that, in order to make the simulated line profiles more realistic, we should “expand” the LOS-velocity axis by some factor, until the “right” ratio of non-thermal velocity range to thermal width is obtained. The expansion factor can be determined by matching the SCF histograms of the simulation with thermal broadening and velocity expansion to those of the HI data. We find that expanding the velocity axis by a factor of 6 gives the best match, as shown in figs. 9a and b, which respectively present the SCF histograms for the thermally-broadened, velocity-expanded simulation data and some selected line profiles. Note also that, since the line profiles of fig. 9b appear significantly more realistic, this figure also shows the extent to which the velocity structure is hidden by the thermal broadening, by comparing the solid lines (thermally broadenend spectra) to the dotted lines (actual velocity histograms along the LOS).

Note that the matching can still be improved, as other differences between the simulation and the NCP Loop data still remain at the next level of refinement: first, although we have matched the mean and standard deviation of the SCF histograms, the skewnesses and kurtoses of the simulation data are considerably smaller than those of the NCP Loop data. This is a reflection that in the NCP Loop channel SCF maps there exist regions of very low values of the SCF corresponding to the arc, which is a strong and large feature, while no such prominent structures are present in the simulation data. Second, the histograms of SCFl have larger (negative) skewnesses than the histograms of SCFs, while the opposite is true for the simulation histograms. As mentioned in §4.3, we have interpreted this as a consequence that in the NCP Loop SCF maps the variability has a larger contribution from intensity variations, while in the simulation data there seems to be a larger contribution from velocity variation along an LOS. Thus, it seems that better matches could be obtained by comparing the simulation data to a region of sky without such a prominent feature, although for brevity we do not pursue this further here, as we believe that the potential of the SCF to guide the simulations to match the observational data has been sufficiently exemplified.

Certainly, the velocity expansion we have done here is not a self-consistent fix to the simulation, as the correct thing to do would be to re-scale the whole simulation to Galactic scales, and to correctly take into account the disk rotation and stratification, as well as more powerful sources of energy (supernovae) capable of accelerating the gas to higher velocities. The conclusion is that our simulations represent a much less turbulent medium than the actual ISM, as implied by their smaller velocity dispersion compared to the velocity dispersion of the H I gas. This exercise has shown the potential of the SCF to pick out differences between the HI and simulation PPV data, and to point towards the physical ingredients that are lacking in the simulations for them to be better models of a given type of observation.

## 6 Summary and Conclusions

In this paper we have presented further discussion on the Spectral Correlation Function (SCF), in addition to that given by RGWW and PRG, and then applied it to data from numerical simulations of the ISM and from HI data (Hartmann & Burton 1997) of the North Celestial Pole (NCP) Loop. The application to the numerical data, for which all the information on the distribution of the physical variables is known, allows us to understand the response of the various SCF measurements used here on the velocity structure in the simulations. This understanding was then used to interpret the results of applying the SCF to the NCP Loop data, and to suggest modifications to the parameters of the ISM simulation to better match the observational data.

In this paper, as in RGWW, the SCF has been used in a “local” form, in which the SCF between one position in a map and its nearest neighbors is calculated and assigned as an “intensity” to that position, thus producing maps of the SCF. Thus, structural features in the SCF maps correspond to sharp gradients in the line profiles’ shape and/or amplitude on the plane of the sky. Also, histograms of the SCF values in the maps have been studied. This procedure has been applied to all four “modes” of the SCF, SCF0, SCFs, SCFl and SCFl,s. The behavior of the SCF with separation, similarly to the more familiar case of regular correlation functions in hydrodynamics (see, e.g., Lesieur 1990), is discussed by PRG.

We first reviewed basic definitions concerning the SCF, in particular the four “modes”, which allow the SCF to focus on spectral similarities in either total intensity, velocity structure, none, or both. We then considered the effect of noise, and showed that in the case of HI (and certainly in the case of numerical simulations) it is not an issue, allowing us to neglect its effects in the subsequent discussion.

In the absence of thermal broadening, the SCF maps of all simulations exhibited a strikingly “wormy” structure, with an interesting kind of self-similarity in which the smallest-scale “worms” bundle up to form larger-scale ones. The “worms” may be associated to caustics in the position-position-velocity (PPV) data cube (loci of high intensity due to velocity crowding), discussed by Henney et al. (2001). Such caustics are not a real feature of the physical fields, but are generated by the sepctroscopic observation process, and they are certainly present in the PPV cube, strongly influencing the structure seen in the channel maps. Their presence causes the SCF histograms to broaden, and to have lower mean values and larger negative skewnesses.

Application of the SCF to a multi-temperature, self-gravitating MHD simulation of the ISM at intermediate scales (3-300 pc), showed that thermal broadening may dominate the spectral lines when the warm gas is primarily subsonic, even if cold, locally supersonic regions exist along the line of sight (assuming optically thin lines), as these regions are still sub- or trans-sonic with respect to the warmer gas. This caused the small-scale structure (which also typically consists of smaller velocity amplitudes) to disappear from the thermally-broadened (TB) SCF maps compared to the non-thermally broadened (NTB) ones, and the TB SCF histograms to be much narrower and with higher mean values than their NTB counterparts. This implies that in the case of the multi-temperature HI data much of the small-scale velocity structure associated with the colder gas may be unobservable in emission, and is possibly only accessible through self-absorption features in the line profiles.

We then applied the SCF to two isothermal runs, one including stellar energy injection, self-gravity and the magnetic field (run ISM-IT), which may be representative of molecular gas, and one without all these ingredients and forced at large scales (run IT). We did not consider thermal broadening for these runs because nearly isothermal (molecular) gas in the ISM is generally very supersonic, at least above scales $`0.1`$ pc. The channel maps for these runs were seen to be markedly different from those for run ISM, lacking the largest-scale structures formed by large-scale gravitational instability in the latter run (run ISM-IT has a Jeans length larger than the simulation size, and run IT has no gravity at all). However, SCF maps of the isothermal runs were not found to be so different, because the “local” SCF used in this paper focuses on small-scale structure.

The SCF histograms of simulations ISM and ISM-IT, which have active star-forming regimes and thus contain stronger, more numerous shocks, tend to have more extended tails at low values of the SCF, and thus large negative skewnesses. This appears to be due to the strong spectral and intensity variations caused by the shocks, suggesting that large negative skewnesses in the SCF histograms is indicative of the presence of shocks and/or shells in the physical fields. Unfortunately, other structures in the gas, or even artifacts in the PPV cube such as the caustics, can also produce such low-SCF tails in the histograms, and thus the indication cannot be completely unambiguous.

Then, SCF maps and histograms of the NCP Loop region were analyzed. The “wormy” structure is seen in the maps for this case also, indicating that the caustics are a real concern. Moreover, the histograms also exhibit large negative skewnesses, due to the prominence of the shell-like structure in this region, in agreement with the suggestion from the simulation data.

We then proceeded to show how the SCF can be used to “fine tune” the simulations to better match the observations. Indeed, the means and standard deviations of the NTB and TB simulation histograms are seen to bracket the corresponding values for the NCP Loop histograms, and in fact the TB line profiles are seen to be completely thermally dominated, with virtually no velocity information left. Attributing this effect to the fact that the ISM simulation contains a rather limited velocity contrast, we were able to obtain a much better match to the NCP Loop SCF histograms, and much more realistic line profiles, by “expanding” the velocity range by a factor of 6, indicating that the simulation requires stronger sources of energy to generate larger velocity dispersion. In other words, we found that the simulations presented here are less energetic than the gas in the NCP region.

We conclude that the local mode of the SCF used here has the ability to quantify the presence of specific velocity and structural features in observational data such as shocks/shells, as well as to allow a quantitative comparison between observational and simulation data.

We are pleased to acknowledge C. Heiles for enlightening and enjoyable discussions. This work was supported in part by NASA Astrophysical Theory Program grant no. NAG5-10103 and CONACYT grants 88046-EUA to J. B.-P., and 27752-E to E.V.-S.


text:
# 1 Introduction

## 1 Introduction

A convenient approach which may be applied to the investigation of the infrared behavior of QCD is the formulation of the non-Abelian gauge theory on loop space that becomes possible by virtue of the close correspondence between gauge and chiral fields . One of the advantages of this method is that it allows one to consider both the perturbative and nonperturbative contributions. The basic object of study in such an approach is the gauge invariant vacuum average of the Wilson loop operator

$$W(C)=\frac{1}{N_c}\text{Tr}0|𝒫\mathrm{exp}\left(ig_C𝑑x_\mu \widehat{A}_\mu (x)\right)|0,$$
(1)
where the integration goes along the closed contour $`C`$ and the gauge field

$$\widehat{A}_\mu (x)=T^aA_\mu ^a(x),T^a=\frac{\lambda ^a}{2},$$
(2)
belongs to the Lie algebra of the gauge group $`SU(N_c)`$, while the Wilson loop operator $`𝒫\text{e}^{ig{\scriptscriptstyle 𝑑xA(x)}}`$ lies (for quark lines) in its fundamental representation.

In particular, the objects such as path-ordered contour integrals of the gauge field emerge naturally in the study of hard collisions of hadrons (Drell–Yan process, deep inelastic scattering, etc.) and the hadron form factors in the eikonal approximation . The use of the Wilson integral formalism in the heavy quarks effective theory has also been discussed (see, e.g., ). In the case of hard processes, the integration contours become infinite, since they present the classical trajectories of the partons participating in the collision. Considering the parton scattering processes at large total energies and small momentum transfers it is revealed that the role of the resumed soft gluons exchanges becomes significant. While in QED this problem is solved in closed form yielding the well-known Sudakov form factor (by virtue of exponentiation of the leading one loop perturbative term), in QCD the situation appears to be more complicated due to its non-Abelian nature and the influence of original nonperturbative effects. Nevertheless, it has been shown in perturbation theory that in QCD processes the exponentiation takes place with the exponent given by a series in the strong coupling constant. Additional soft contributions may appear due to various nonperturbative effects, a number of which are closely related to the nontrivial structure of the QCD vacuum. In the case of the instanton ensemble modeling the QCD vacuum, the exponentiation of the single instanton contribution can also be proven (see Section 4). Indeed, although the QCD vacuum plays an important role in the high-energy collisions, the direct investigation of these effects remains a difficult task. The aim of the present work is to evaluate them within the framework of the instanton liquid model of QCD vacuum for the simplest possible quantity – the quark form factor.

It is generally believed that the nonperturbative structure of QCD can be well understood within the instanton liquid model . Considering the QCD vacuum as an ensemble of instantons, one can describe a lot of the low-energy phenomena in strong interactions on the qualitative, as well as the quantitative level. The importance of the instanton induced effects of strong interactions is also supported by lattice simulations . The instanton picture is generally considered as a fruitful and perspective framework for hadronic physics. The role of the instantons in hard hadronic processes has been studied intensively, both theoretically and experimentally. The perspectives for an unambiguous experimental detection of instanton contributions are believed to be optimistic and very promising. The main problem arising in this context is the evaluation of instanton induced effects at a hard scale of a given process . Hence, we are going to study the instanton induced contributions to the Wilson integral over contours of a various geometry and apply the obtained results to the analysis of nonperturbative effects in the high-energy hadron collisions. Recently a similar approach has been developed in the study of high energy elastic and quasielastic parton-parton scattering . The vacuum averages of Wilson operators with finite-length contours were considered in the instanton model in , where the properties of the nonlocal quark and gluon condensates were discussed.

In the present paper, we start the investigation of the instanton induced effects in the high-energy QCD processes. We propose an approach which allows one to evaluate the instanton contributions to the Wilson integrals made of several (in)finite lines containing specific cusp and/or cross singularities. For this purpose, we start with one of the simplest configurations, i. e., the angle with infinite sides that corresponds to the integration path for the Wilson operator describing the soft part of the quark Sudakov form factor . Such configurations enter as constructive units in the Wilson contours for many high energy processes.

This paper is organized as follows. In Section 2, we briefly describe the general renormalization properties of the Wilson loops with cusp singularities and find the explicit formula in leading order, applicable both to perturbative and nonperturbative fields. Then we analyze the UV and IR divergencies and calculate the perturbative part in an arbitrary covariant gauge demonstrating explicitly the gauge invariance of the obtained (renormalized) result. Section 3 is devoted to the study of the instanton induced nonperturbative contributions. The expression for them is found in the closed form in the weak-field limit in terms of the instanton profile function. It is shown that the instantons actually give powerlike corrections to the perturbative result. The magnitude of the instanton induced effects is found to be comparable with the perturbative part at a certain low momentum scale. In Section 4, we prove the exponentiation in the small instanton density parameter of the all-order single instanton contribution to the Sudakov form factor. Afterwards, we discuss the possible applications of the results obtained above to investigations of the physically interesting processes and make some conclusive remarks.

## 2 Renormalization and perturbative contribution

Let us consider in short the renormalization properties of the Wilson loop. It is well known, that evaluation of $`W(C)`$ (1) encounters divergencies. The divergencies appearing for smooth contours without self-intersections can be removed by means of the convenient R operation :

$$W_R(C,g_R,\mu )=\underset{\epsilon 0}{lim}\stackrel{~}{W}(C,g_R,\mu ,\epsilon )\underset{\epsilon 0}{lim}RW(C,g,\epsilon ),$$
(3)
where $`\epsilon `$ is the dimensional regularization parameter, $`\mu `$ is the UV normalization scale. In this case, the R operation consists in the multiplicative renormalization of the coupling constant $`gg_R`$ and the field $`\widehat{A}_\mu \widehat{A}_\mu ^R`$.

In the more general case, the contour may contain a number of cusps, i. e., the points where the derivative is not smooth. In this situation, divergencies of another type arise which depend on the corresponding cusp angles $`\gamma _i`$. Although they were considered as a shortcoming of the loop formulation of a gauge theory about two decades ago, these singularities (and the anomalous dimensions related to them), are shown to play an important role in the partonic hard processes controlling the asymptotics of the scattering amplitudes . In order to remove these divergencies, the generalized $`K_\gamma `$ operation has been proposed . For a contour with one cusp $`C_\gamma `$ this operation is reduced to

$$W_R(C_\gamma ,g_R,\mu ,\overline{C}_\gamma )=\underset{\epsilon 0}{lim}K_\gamma \stackrel{~}{W}(C_\gamma ,g_R,\mu ,\epsilon ),$$
(4)
where the contour $`\overline{C}_\gamma `$ is treated as a subtraction point. It has been proved that the $`K_\gamma R`$ operation renormalizes multiplicatively any loop integral with finite number of cusps . The multiplicative renormalizability means that the $`K_\gamma `$ operation acts as

$$K_\gamma \stackrel{~}{W}(C_\gamma ,g_R,\mu ,\epsilon )=Z_{cusp}(g_R,\gamma ,\mu ,\epsilon )\stackrel{~}{W}(C_\gamma ,g_R,\mu ,\epsilon ),$$
(5)
where

$$Z_{cusp}=\stackrel{~}{W}^1(\overline{C}_\gamma ,g_R,\mu ,\epsilon )=1+\underset{n=1}{\overset{\mathrm{}}{}}\left(\frac{\alpha _S}{\pi }\right)^nZ_n(\gamma ,\mu ,\epsilon )$$
(6)
is an additional renormalization constant. The coefficients $`Z_n`$ are the sum of simple poles plus some finite terms

$$Z_n(\gamma ,\epsilon )=\underset{k=1}{\overset{n}{}}\left(\frac{1}{\epsilon ^k}A_{kn}(\gamma )\right)+(\text{finite    terms})_n.$$
(7)
A choice of the finite part of $`Z_n`$ to each order actually determines the $`K_\gamma `$ subtraction scheme. In order to remove pole terms we shall use the modified minimum subtraction $`\overline{MS}`$-scheme $`(K_\gamma ^{\overline{MS}})`$.

The first nontrivial cusp dependent term in the expansion of $`W(C)`$ (1) in powers of $`g^2`$ for the angle with two infinite straight line rays \[(Fig. 1a)\] contains the contributions from both perturbative \[(Fig. 1b)\] and nonperturbative \[(Fig. 1c)\] fields

$$W^{(1)}(\gamma )=W_P^{(1)}(\gamma )+W_{NP}^{(1)}(\gamma ).$$
(8)
Regardless of which part is considered – either perturbative, or nonperturbative – we will write the one-loop contribution in the following form:

$$W^{(1)}(\gamma )=\frac{g^2C_F}{2}_{C_\gamma }𝑑x_\mu _{C_\gamma }𝑑y_\nu D_{\mu \nu }(xy),$$
(9)
where $`C_F=\frac{N_c^21}{2N_c}`$. It is convenient to present the gluon propagator $`D_{\mu \nu }(z)`$ in the form

$$D_{\mu \nu }(z)=\delta _{\mu \nu }_z^2d_1(z^2)_\mu _\nu d_2(z^2).$$
(10)
Here and in what follows, we use the dimensional regularization with $`n=42\epsilon `$, $`\epsilon <0`$ in order to control the IR-divergent terms in the integrals. The remaining UV singularity (due to the infinitely small $`z^2`$ in the vicinity of the cusp) will be regularized by the corresponding UV cutoff.

The trajectories of the incoming and outgoing quarks \[(Fig. 1a)\] may be parametrized as

$$x=v_1s(0<s<\mathrm{}),y=v_2\tau (\mathrm{}<\tau <0).$$

The angle between the vectors $`v_1`$ and $`v_2`$ is given in the Minkowski space by

$$\mathrm{cosh}\chi =(v_1v_2)=\frac{(p_1p_2)}{m^2}=1+\frac{Q^2}{2m^2},Q^2=(p_2p_1)^2,v_{1,2}^2=1,$$
(11)
where the quark momenta are supposed to be on shell: $`p_1^2=p_2^2=m^2`$. The continuation to the Euclidean space is defined as :

$$\chi i\gamma .$$
(12)
By virtue of the interpretation of the graph Fig. 1 as an amplitude of the elastic scattering of an on-mass-shell one-dimensional fermion on a color singlet potential , we have to consider the quantity

$$W^{(1)}(\gamma )=\stackrel{~}{W}^{(1)}(\gamma )\stackrel{~}{W}^{(1)}(0),$$
(13)
where (for some technical details of the calculations, see the Appendix)

$$\stackrel{~}{W}^{(1)}(\gamma )=g^2C_F\left[(n2)d_1(0)\gamma \text{cot}\gamma +d_2(0)\right].$$
(14)
It follows from Eq. (13) that the integrals (9) in which both points $`x`$ and $`y`$ belong to the same side of the angle do not contribute to the quantity $`W^{(1)}(\gamma )`$. Hence we have within the one loop accuracy

$$W(\gamma )=14\pi \alpha _SC_F(n2)h(\gamma )d_1(0),$$
(15)
where

$$h(\gamma )=\gamma \text{cot}\gamma 1$$
(16)
is the universal cusp factor. We should emphasize here that the expression (15) holds for perturbative part as well as for the nonperturbative part depending on the value $`d_1(0)`$. For the perturbative field, Eq. (15) reflects the explicit gauge invariance in the set of covariant gauges, since the gauge fixing parameter $`\xi `$ enters only in the function $`d_2(z^2)`$.

Let us consider first the perturbative part $`W_P^{(1)}(\gamma )`$. By using the free propagator in the Euclidean space (see the Appendix) the IR-regularized value of $`d_1(0)`$ can be written in the form

$$d_1(0;\epsilon ,\lambda )=\frac{(\lambda ^2\pi )^\epsilon }{16\pi ^2}_0^{\mathrm{}}𝑑\alpha \alpha ^{(1+\epsilon )},$$
(17)
where $`\lambda ^2`$ is the IR regularization parameter. This integral diverges at the upper (UV) limit for $`\epsilon <0`$, and hence we must regularize it. To this end, we may introduce the UV cutoff $`\mu ^2`$, that corresponds to the following replacement in the denominator of the perturbative propagator: $`z^2z^2+\mu ^2`$, and finally we get

$$d_1(0;\epsilon ,\mu /\lambda )=\frac{1}{\epsilon }\frac{1}{16\pi ^2}\left(\frac{\lambda ^2\pi }{\mu ^2}\right)^\epsilon .$$
(18)
Thus one obtains, with one loop accuracy,

$$W_P(\gamma ,\alpha _S,\mu /\lambda ,\epsilon )=1+\frac{\alpha _S}{2\pi }C_Fh(\gamma )\frac{1}{\epsilon }(1\epsilon )\left(\frac{\lambda ^2\pi }{\mu ^2}\right)^\epsilon ,$$
(19)
and the cusp renormalization constant $`Z_{cusp}^{(1)}`$ within the $`\overline{MS}`$-scheme:

$$Z_{cusp}^{(1)}(\gamma ,\epsilon )=1\frac{\alpha _S}{2\pi }C_Fh(\gamma )\left(\frac{1}{\epsilon }+\text{ln}\pi 1\right),$$
(20)
which is in agreement with . Therefore, the perturbative part of the finite renormalized function reads

$$W_P(\gamma ,\mu )=\underset{\epsilon 0}{lim}Z_{cusp}^{(1)}(\gamma ,\mu ,\epsilon )\stackrel{~}{W}_P(\gamma ,\mu ,\epsilon )=1\frac{\alpha _S}{2\pi }C_Fh(\gamma )\text{ln}\frac{\mu ^2}{\lambda ^2}.$$
(21)
The one-loop cusp anomalous dimension which satisfies the RG equation reads

$$\left(\mu \frac{}{\mu }+\beta (g_R)\frac{}{g_R}\right)W_P(\gamma ,\mu ,\alpha _S(\mu ))=\mathrm{\Gamma }_{}^{cusp}{}_{P}{}^{(1)}(\gamma ,\alpha _S(\mu ))W_P(\gamma ,\mu ,\alpha _S(\mu )),$$
(22)
where $`\beta (g_R)`$ is the usual QCD $`\beta `$-function $`\beta (g_R)=\mu \frac{}{\mu }g_R=\text{const}g_R^3+O(g_R^5)`$, which can be found in the one-loop order

$$\mathrm{\Gamma }_{}^{cusp}{}_{P}{}^{(1)}(\gamma ,\alpha _S)=\frac{d}{d\text{ln}\mu }W_P^{(1)}(\gamma ,\mu )=\frac{\alpha _S}{\pi }C_Fh(\gamma ).$$
(23)
Hence, we reproduce the Wilson operator value for the infinite contour with the Euclidean cusp parameter $`\gamma `$.

## 3 Instanton contribution

Let us estimate the nonperturbative contribution to $`W(C)`$ in the instanton model. The instanton field is given by

$$\widehat{A}_\mu (x;\rho )=A_\mu ^a(x;\rho )\frac{\sigma ^a}{2}=\frac{1}{g}\text{R}^{ab}\sigma ^a\eta _{}^{\pm }{}_{\mu \nu }{}^{b}(xz_0)_\nu \phi (xz_0;\rho ),$$
(24)
where $`\text{R}^{ab}`$ is the color orientation matrix $`(a,b=1,2,3)`$, $`\sigma ^a`$’s are the Pauli matrices, and $`(\pm )`$ corresponds to the instanton or antiinstanton. The averaging of the Wilson operator over the nonperturbative vacuum is reduced to the integration over the coordinate of the instanton center $`z_0`$, the color orientation and the instanton size $`\rho `$. The measure for the averaging over the instanton ensemble reads $`dI=d\text{R}d^4z_0dn(\rho )`$, where $`d\text{R}`$ refers to the averaging over color orientation, and $`dn(\rho )`$ depends on the choice of the instanton size distribution. Taking into account (24), we write the Wilson integral (1) in the single instanton approximation in the form

$$w_I(C)=\frac{1}{N_c}0|\text{Tr}\mathrm{exp}\left(i\sigma ^a\varphi ^a\right)|0,$$
(25)
where the phase is

$$\varphi ^a=\text{R}^{ab}\eta _{}^{\pm }{}_{\mu \nu }{}^{b}_{C_\gamma }𝑑x_\mu (xz_0)_\nu \phi (xz_0;\rho ).$$
(26)
We omit the path ordering operator $`𝒫`$ in (25) because the instanton field (24) is a hedgehog in color space, and so it locks the color orientation by space coordinates. Performing the averaging over the color rotations, making use of the known property of ’t Hooft symbols,

$$\eta _{}^{\pm }{}_{\mu \nu }{}^{a}\eta _{}^{\pm }{}_{\rho \sigma }{}^{a}=\delta _{\mu \rho }\delta _{\nu \sigma }\delta _{\mu \sigma }\delta _{\nu \rho }\pm \epsilon _{\mu \nu \rho \sigma },$$
(27)
and performing the subtraction of self-energy part (13) we obtain the all-order single instanton contribution to the cusp-dependent part of Wilson loop (1):

$$w_I(\gamma )=d^4z_0𝑑n(\rho )\left[\text{cos}\varphi (\gamma ,z_0,\rho )\text{cos}\varphi (0,z_0,\rho )\right],$$
(28)
where the squared phase $`\varphi ^2=\varphi ^a\varphi ^a`$ may be written as

$$\varphi ^2(\gamma ,z_0,\rho )=\underset{i,j=1,2}{}[(v_iv_j)z_0^2(v_iz_0)(v_jz_0)]\times $$

$$_0^{\mathrm{}}𝑑\sigma \phi [((1)^{i+1}v_i\sigma z_0)^2;\rho ]_0^{\mathrm{}}𝑑\sigma ^{}\phi [((1)^{j+1}v_j\sigma ^{}z_0)^2;\rho ].$$
(29)
Although the expressions (28, 29) give the complete formula for the all-order single instanton contribution, in what follows we restrict ourselves to the investigation of the weak field limit. Phenomenologically, it is assumed that the instanton size distribution is sharply peaked at a certain finite value. As we will see soon, this finite instanton size $`\rho `$ provides an UV cutoff for the nonperturbative gluon propagator. Therefore, for the instanton field we use the same dimensional IR regularization with $`\epsilon <0`$ and the IR parameter $`\lambda `$ as in the perturbative case333 Another possible approach to divergencies in the integrals over the instanton fields is to use a constrained instanton solution which is characterized by the profile function exponentially decreasing at large distances .. So, in this limit the leading instanton induced term reads

$$w_I^{(1)}(\gamma )=\frac{g^2\lambda ^{n4}}{2}𝑑n(\rho )_{C_\gamma }𝑑x_\mu _{C_\gamma }𝑑y_\nu \frac{d^nk}{(2\pi )^n}\stackrel{~}{A}_\mu ^a(k;\rho )\stackrel{~}{A}_\nu ^a(k;\rho )\text{e}^{ik(xy)}.$$
(30)
By using the Fourier transform of the instanton field

$$\stackrel{~}{A}_\mu ^a(k;\rho )=\frac{2i}{g}\eta _{}^{\pm }{}_{\mu \sigma }{}^{a}k_\sigma \stackrel{~}{\phi }^{}(k^2;\rho ),$$
(31)
Eq. (30) can be written in the form of Eq. (10) with the instantonic analogue of the function $`d_1(z^2)`$:

$$d_1(z^2)d_1^I(z^2)=\frac{1}{g^2C_F}𝑑n(\rho )D_I(z^2;\rho ,\lambda ),$$
(32)
where

$$D_I(z^2;\rho ,\lambda )=\lambda ^{4n}\frac{d^nk}{(2\pi )^n}\text{e}^{ikz}\left[2\stackrel{~}{\phi }^{}(k^2;\rho )\right]^2.$$
(33)
Above, $`\stackrel{~}{\phi }(k^2;\rho )`$ is the Fourier transform of the instanton profile function $`\phi (z^2;\rho )`$ and $`\stackrel{~}{\phi }^{}(k^2;\rho )`$ is it’s derivative with respect to $`k^2`$. Now using the result (15) of the previous section, we get the instanton contribution in the form

$$w_I(\gamma ;\epsilon ,\lambda )=(n2)h(\gamma )𝑑n(\rho )D_I(0;\epsilon ,\lambda ,\rho ).$$
(34)
Consider now the renormalization of the nonperturbative part for the instanton field in the singular gauge,444 We consider the singular gauge, since in the regular gauges the profile function decreases insufficiently rapidly at infinity, and therefore the contribution of the infinitely distant part of the contour must be taken into account, in contrast to the singular gauge where this contribution vanishes. where the profile function is

$$\phi (u;\rho )=\frac{\rho ^2}{z^2(z^2+\rho ^2)}.$$
(35)
By using the Fourier transform of this function one gets the IR dimensionally regularized $`D_I(0;\epsilon ,\lambda ,\rho )`$ in the form

$$D_I(0;\epsilon ,\lambda ,\rho )=$$

$$=\frac{\rho ^4\pi ^{2+\epsilon }\lambda ^{2\epsilon }}{4}_0^{\mathrm{}}𝑑\alpha _1𝑑\alpha _2_0^{\mathrm{}}𝑑\beta _1𝑑\beta _2\frac{\text{e}^{\rho ^2(\alpha _1+\beta _1)}}{(\alpha _1+\alpha _2)^{1+\epsilon }(\beta _1+\beta _2)^{1+\epsilon }\left(\alpha _1+\alpha _2+\beta _1+\beta _2\right)^{2\epsilon }}.$$
(36)
It is clear from this expression that the instanton size plays a role of the UV cutoff, since it provides an exponential suppression of the integrand at large $`(\alpha ,\beta )`$ what corresponds to the large momentum. Computing this integral we have

$$D_I(0;\epsilon ,\lambda ,\rho )=$$

$$\frac{\rho ^4\pi ^2}{4}\left(\rho ^2\lambda ^2\pi \right)^\epsilon \frac{\mathrm{\Gamma }(1\epsilon )}{\epsilon }\left(1+\epsilon _0^1_0^1_0^1𝑑x𝑑y𝑑z\text{ln}\left[\frac{xz+y(1z)}{z(1z)}\right]+O(\epsilon ^2)\right).$$
(37)
Applying the renormalization procedure as described in the previous section, we find in the leading order the instanton contribution to the renormalization constant

$$Z_I(\gamma ,\epsilon )=12h(\gamma )\left(\frac{1}{\epsilon }𝑑n(\rho )\frac{\rho ^4\pi ^2}{4}+\mathrm{\Delta }_{NP}\right),$$
(38)
where the finite term is

$$\mathrm{\Delta }_{NP}=\gamma _E+\text{ln}\pi +\frac{\pi ^2}{6}\frac{3}{2},$$
(39)
and obtain the corresponding contribution to the Wilson loop

$$w_I^{(1)}(\gamma ,\lambda )=1+\pi ^2h(\gamma )𝑑n(\rho )\rho ^4\text{ln}\left(\rho \lambda \right).$$
(40)
In order to estimate the magnitude of the instanton induced effect we consider the distribution function which has been suggested in  (and discussed in  in the framework of constrained instanton model) in order to describe the lattice data , namely the following one

$$dn(\rho )=\frac{d\rho }{\rho ^5}C_{N_c}\left(\frac{2\pi }{\alpha _S(\rho )}\right)^{2N_c}\mathrm{exp}\left(\frac{2\pi }{\alpha _S(\rho )}\right)\mathrm{exp}\left(2\pi \sigma \rho ^2\right),$$
(41)
where the numerical constant $`C_{N_c}`$ is determined by the number of colors

$$C_{N_c}=\frac{0.466\text{e}^{1.679N_c}}{(N_c1)!(N_c2)!}0.0015,$$
(42)
and the string tension is accepted to be $`\sigma (0.44GeV)^2`$ . Then, using the one loop expression for the running coupling constant

$$\alpha _S(\rho )=\frac{2\pi }{\beta _0\text{ln}\rho \mathrm{\Lambda }_{QCD}},\beta _0=\frac{11N_c2n_f}{3},$$
(43)
we find the instanton contribution (40) in the form \[in the distribution (41), the slow varying logarithmic factor due to the power of the coupling $`\alpha _S`$ is assumed to be constant, and taken at the point of the mean instanton size $`\overline{\rho }`$\]

$$w_I^{(1)}(\gamma ,\lambda )=1+\pi ^2h(\gamma )\frac{C_{N_c}\mathrm{\Gamma }(\beta _0/2)}{4}\left(\frac{2\pi }{\alpha _S(\overline{\rho })}\right)^{2N_c}\left(\frac{\mathrm{\Lambda }_{QCD}}{\sqrt{2\pi \sigma }}\right)^{\beta _0}\text{ln}\frac{\lambda ^2}{2\pi \sigma },$$
(44)
where we omitted the constant term on the right hand side since it can be easily removed by the finite renormalization. The expression (44) shows explicitly that the instantons yield the powerlike corrections to the perturbative result, which is expected from a general consideration. Indeed, the powerlike behavior of the nonperturbative corrections to the high energy QCD processes may be obtained from the renormalon analysis (see, e.g., ).

It is instructive to express the result (44) in terms of the mean instanton size $`\overline{\rho }`$ and the instanton density $`\overline{n}`$ calculated directly from the distribution (41), which read respectively

$$\overline{\rho }=\frac{\mathrm{\Gamma }(\beta _0/23/2)}{\mathrm{\Gamma }(\beta _0/22)}\frac{1}{\sqrt{2\pi \sigma }},$$
(45)
$$\overline{n}=\frac{C_{N_c}\mathrm{\Gamma }(\beta _0/22)}{2}\left(\frac{2\pi }{\alpha _S(\overline{\rho })}\right)^{2N_c}\left(\frac{\mathrm{\Lambda }_{QCD}}{\sqrt{2\pi \sigma }}\right)^{\beta _0}(2\pi \sigma )^2.$$
(46)
Numerically, these quantities correspond to the usual mean size $`\rho _0`$ and density $`n_0`$ of the instanton liquid model, which parametrize the deltalike approximation to the distribution (41) (this approach, which was proposed about two decades ago , now gets a strong support in lattice data ):

$$dn(\rho )=n_0\delta (\rho \rho _0)d\rho ,$$
(47)
where the parameters were estimated in :

$$n_01fm^4,\rho _01/3fm.$$
(48)
To compare the instanton induced and perturbative parts, which are, respectively,

$$\{\begin{array}{c}w_I(\gamma )=1Kh(\gamma )\pi ^2\overline{n}\overline{\rho }^4\left[\text{ln}\frac{\mu ^2}{\lambda ^2}\text{ln}\overline{\rho }^2\mu ^2\right]\hfill \\ W_P(\gamma )=1\frac{\alpha _S(\mu )}{2\pi }C_Fh(\gamma )\text{ln}\frac{\mu ^2}{\lambda ^2},\hfill \end{array}$$
(49)
where

$$K=\frac{\mathrm{\Gamma }(\beta _0/2)[\mathrm{\Gamma }(\beta _0/22)]^3}{2[\mathrm{\Gamma }(\beta _0/23/2)]^4}0.74,$$
(50)
we assume that the factorization scale $`\mu `$ (which divides the soft and hard regions of momenta in the factorized quark form factor) is of order of the inverse instanton size $`\mu \overline{\rho }^10.6GeV`$. Then we write the total leading order contribution to the Wilson loop expectation value in the form

$$W(\gamma ,\rho _0\lambda )=1+\frac{\alpha _S(\overline{\rho }^1)}{2\pi }C_Fh(\gamma )\text{ln}\left(\overline{\rho }^2\lambda ^2\right)\left(1+K\frac{S_0\pi ^2\overline{n}\overline{\rho }^4}{C_F}\right),$$
(51)
where $`S_0=\frac{8\pi ^2}{g^2(\overline{\rho }^1)}10`$ is the “classical enhancement” factor with the renormalized coupling constant $`g(\mu )`$ at the energy scale $`\mu \overline{\rho }^1`$. The ratio of the instanton correction to the perturbative leading term is given by

$$\frac{\text{Instanton}}{\text{1-Loop Perturb.}}=K\frac{S_0\pi ^2\overline{n}\overline{\rho }^4}{C_F}0.5,$$
(52)
which is estimated using the conventional value for the packing fraction

$$\pi ^2\overline{n}\overline{\rho }^40.1.$$
(53)
One can see using the main formula (44) that the strong power suppression of the instanton part is partially compensated by the large factor $`S_0^{2N_c+1}`$. This means that at the energy scale of order of $`\overline{\rho }^1`$ the magnitude of the instanton induced effects is comparable to the leading perturbative part, and must be taken into account as well.

Let us estimate now the two-loop perturbative contribution which may be equally important at the chosen low scale, where the strong coupling constant is not small enough \[$`\alpha _S(0.6GeV)0.5`$\]. The two-loop result has been obtained in  and in the limits of small and large angles is given by

$$W^{(2)}(\gamma )\{\begin{array}{c}\left(\frac{\alpha _S}{\pi }\right)^2C_FN_c\frac{\gamma ^2}{3}\left[\frac{11}{48}\text{ln}^2\frac{\mu ^2}{\lambda ^2}\left(\frac{3}{8}\frac{\pi ^2}{12}+\frac{67}{72}\right)\text{ln}\frac{\mu ^2}{\lambda ^2}\right]\mathrm{for}\mathrm{small}\gamma ,\hfill \\ \left(\frac{\alpha _S}{\pi }\right)^2C_FN_c\gamma \left[\frac{11}{48}\text{ln}^2\frac{\mu ^2}{\lambda ^2}\left(\frac{67}{72}\frac{\pi ^2}{24}\right)\text{ln}\frac{\mu ^2}{\lambda ^2}\right]\mathrm{for}\mathrm{large}\gamma .\hfill \end{array}$$
(54)
From these expressions we get that the ratio of the instanton contribution to the perturbative two-loop term is approximately equal to 1.4 and varies slightly (only about $`10\%`$) with changing of $`\gamma `$ over whole allowed kinematical interval. Thus, the complete consideration of the quark form factor at the low momentum scale must include both the two-loop perturbative part and the leading order instanton one, which appear to be of the same order of magnitude.

## 4 Exponentiation of perturbative and instanton corrections

Expressions (51, 54) define the first terms of the Wilson loop expansion in gauge fields. On the basis of the exponentiation theorem  for the non-Abelian path-ordered exponentials it is well known that perturbative corrections to the Sudakov form factor are exponentiated to high orders in the QCD coupling constant. The theorem states that the contour average $`W_P(C)`$ can be expressed as

$$W_P(C)=\mathrm{exp}\left[\underset{n=1}{\overset{\mathrm{}}{}}\left(\frac{\alpha _S}{\pi }\right)^n\underset{WW(n)}{}C_n(W)F_n(W)\right],$$
(55)
where summation in the exponential is over all diagrams $`W`$ of the set $`W(n)`$ of the two-particle irreducible contour averages of $`n`$th order of the perturbative expansion. The coefficients $`C_n(W)C_FN_c^{n1}`$ are the “maximally non-Abelian” parts of the color factor corresponding to the contribution coming from a diagram $`W`$ to the total expression (55) in the contour gauge, and the factor $`F_n(W)`$ is the contour integral presented in the expression for $`W`$. This means that the essential diagrams are only those, which do not contain the lower-order contributions as subgraphs and, as a result, the higher-order terms are non-Abelian.

Let us now demonstrate how the single instanton contribution is exponentiated in the small instanton density parameter, treating the instanton vacuum as a dilute medium . The gauge field is taken to be the sum of individual instanton fields in the singular gauge (24, 35), with their centers at the points $`z_j`$’s. In this gauge, the instanton fields fall off rapidly at infinity, so the instantons may be considered individually in their effect on the loop. Moreover, the contribution of infinitely distant parts of the contour may be neglected and only those instantons will influence the loop integral, which occupy regions of space-time intersecting with the quark trajectories. Since the parametrization of the loop integral along rays of the angle plays the role of the proper time, a time-ordered series of instantons arises and has an effect on the Wilson loop. Thus, the contribution of $`n`$ instantons to the loop integral $`W_I(\gamma )`$ can be written in the dilute approximation as

$$W_I^{(n)}(\gamma )=\text{Tr}\left(U^1U^2\mathrm{}U^nU^n\mathrm{}U^2U^1\right),$$
(56)
where the ordered line integrals $`U_i`$’s

$`U^j(\gamma )`$ $`=`$ $`T\left\{\mathrm{exp}\left(ig{\displaystyle _0^{\mathrm{}}}𝑑\sigma v_1^\mu A_\mu (v_1\sigma z_j)\right)\mathrm{exp}\left(ig{\displaystyle _{\mathrm{}}^0}𝑑\tau v_2^\mu A_\mu (v_2\tau z_j)\right)\right\}`$
are associated with individual instantons with the positions $`z_j`$’s. Because of the wide separation of the instantons in the dilute phase and rapid falloff of fields in the singular gauge, the upper and lower limits of the line integrals are extended to infinity. The line integrals $`U^i`$’s take into account the part of the contour that goes at infinity from $`+\mathrm{}`$ back to $`\mathrm{}`$ and in the singular gauge $`U^i=1`$. For $`U^j(U^j)`$, the integral is taken over the increasing (decreasing) time piece of the loop.

Then, the expression is simplified when averaging over the gauge orientations of instantons. The averaging is reduced to substitution of $`U^j`$ by $`g_jU^jg_j^1`$, where $`g_j`$ is an element of colour group, and independent integration of each $`g_j`$ over the properly normalized group measure is performed. Under this averaging one gets

$$U^nU^n\frac{1}{N_c}\text{Tr}\left(U^nU^n\right),$$
(57)
which is just the single instanton contribution $`w_I^{(n)}(\gamma )`$ as it is given by Eqs. (25, 28). But then if the averaging is done in the inverse order, from $`n`$ down to $`1`$, the entire loop integral collapses to a product of traces

$$W_I^{(n)}(\gamma )\underset{n\mathrm{}}{lim}\underset{j=1}{\overset{n}{}}w_I^{(j)}(\gamma ).$$
(58)
Since the individual instantons are considered to be decoupled in the dilute medium, the total multiple instanton contribution to the vacuum average of the Wilson operator simply exponentiates the all-order single instanton term $`w_I(\gamma )`$ in (28), and one has

$`W_I(\gamma )`$ $`=`$ $`\underset{n\mathrm{}}{lim}\left\{1+{\displaystyle \frac{1}{n}}w_I(\gamma )\right\}^n=\mathrm{exp}[w_I(\gamma )].`$ (59)
Thus, we have proved that in the dilute regime, the full instanton contribution to the quark form factor is given by the exponent of the all-order single instanton result. The exponentiation arises due to taking into account the many-instanton configurations effect. As it is well known, in QED there occurs the exponentiation of the one-loop result due to Abelian character of the theory. In the instanton case, the analogous result takes place since instantons belong to the $`SU(2)`$ subgroup of the $`SU(3)`$ color group and the path-ordered exponents coincide with the ordinary ones. Exponentiation applied to the lowest order terms (51, 54) leads to vanishing the Wilson loop as a power of an IR cutoff $`\lambda `$. Then, the additional instanton contributions suggest that the Wilson loop decreases in the instanton medium faster with $`\lambda `$ than in purely perturbative theory.

## 5 Conclusion

We have calculated the instanton contribution to the soft part of the quark form factor, described in terms of the vacuum expectation value of the Wilson loop for the contour of a special form (28, 29). Further analysis of high energy behavior of this form factor must take into account the contributions of the hard and collinear parts as well. In order to do this, we would need to “free” the scale $`\mu `$ and consider it as a factorization parameter lying between the hard scale $`Q`$ and the IR cutoff $`\lambda `$ ($`\lambda <\mu <Q`$). Then we would find that the instantonic part, being UV finite due to the finite instanton size, does not contain an explicit dependence on the UV cutoff $`\mu `$, and therefore does not contribute to the cusp anomalous dimension. However, this issue is out of scope of the present paper, and we did not address this question here.

The situation may be changed essentially in some other hadronic hard processes. For example, the Wilson path for (in)elastic high energy quark-quark scattering consists of two infinite straight lines lying in two different planes, separated by the finite transverse distance . The calculation of the Wilson operator for this contour reduces formally to the evaluation of the path integrals analogous to the ones considered above. The only difference is that the sides of these contours would belong to the different planes. The transverse distance between them provides then an UV cutoff for the integrals, both for perturbative and instanton fields, and the parameter of the dimensional regularization $`\lambda `$ becomes the only (IR) cutoff for them. Therefore, the corresponding anomalous dimension would contain the instanton induced part, and the instantons would play a nontrivial role in the solution of the corresponding renormalization group equations. This topic will be a subject of the forthcoming study.

Moreover, in the instanton calculations of the high energy parton scattering processes , the integrals over the instanton field can be computed in all orders by virtue of the infinite length of lines that form the integration path specific for this process. The latter fact leads to the trivial factorization of the $`\gamma `$-angle dependence from the corresponding Wilson integral, and extract the instanton contribution into the numerically calculable form factors. In our case, one deals with the cusp, rather than cross singularities. The all-order calculation given formally by Eqs.(28, 29) is quite complicated for the contour with cusp, so we postpone the solution of this task.

To summarize, within the instanton vacuum model we have developed an approach which allows one to calculate the nonperturbative contributions to the Wilson integrals over the infinite contour with a cusp that represent, e. g., the classical trajectories of partons participating in hard collisions. We have proved that in the dilute regime, the full instanton contribution to the quark form factor is given by the exponentiated all-order single instanton result, see (59). In the weak-field limit, the instanton contribution to the soft part of the color singlet quark form factor is found explicitly in terms of the instanton profile function in the singular gauge. It is shown that the instanton induced effects are of a power type (44), but nevertheless they are comparable in magnitude to the perturbative ones at the scale of order of the inverse average size of the instanton in the instanton vacuum, see (51).

## 6  Acknowledgments

We are grateful to D. Antonov, N.I. Kochelev, S.V. Mikhailov, and S. Tafat. This work is partially supported by RFBR (Grant Nos. 02-02-16194, 01-02-16431) and INTAS (Grant No. 00-00-366). The work of I. Ch. was also supported by RFBR Grant No. 00-15-96577. He is grateful to the Abdus Salam ICTP in Trieste for the kind invitation and hospitality, where a part of this work was done.

## 7  Appendix

Below we write down some useful formulas, that have been used in the calculations. The dimensionally regularized free propagator in the Euclidean space reads $`(n=42\epsilon )`$:

$$D_{\mu \nu }(z;\xi )=\lambda ^{4n}\frac{d^nk}{(2\pi )^n}\text{e}^{ikz}\left(\frac{\delta _{\mu \nu }}{k^2}\xi \frac{k_\mu k_\nu }{k^4}\right)=$$

$$=\frac{(\lambda ^2\pi )^{(2n/2)}}{4\pi ^2}\mathrm{\Gamma }(n/21)\left(\delta _{\mu \nu }\frac{(1\xi /2)}{z^{2(n/21)}}\xi z_\mu z_\nu \frac{(n/21)}{z^{2(n/2)}}\right).$$
$`(1a)`$
The form factors $`d_{1,2}(u)`$ defined in Eq. (10) are given by ($`u=z^2`$ )

$$d_1(u)=\frac{(\lambda ^2\pi )^{(2n/2)}}{16\pi ^2}\frac{\mathrm{\Gamma }(n/22)}{u^{(n/22)}},d_2(u)=\xi d_1(u).$$
$`(1b)`$
The partial derivatives in terms of the derivative with respect to the interval $`u`$ in $`n`$-dimensional space-time

$$_z^2=2n_u+4u_u^2,$$
$`(2a)`$
$$_\mu _\nu =2\delta _{\mu \nu }_u+4z_\mu z_\nu _u^2.$$
$`(2b)`$
Scalar products in the Euclidean space

$$v_\mu ^1v_\nu ^2\delta _{\mu \nu }=\text{cos}\gamma ,v_\mu ^1v_\nu ^2z_\mu z_\nu =z^2\text{cos}\gamma +st\text{sin}^2\gamma .$$
$`(3)`$
For calculations of the perturbative part we used the following integrals:

$$_0^{\mathrm{}}𝑑s𝑑t\text{e}^{\alpha (s^2+t^2+2st\text{cos}\gamma )}=\frac{1}{2\alpha }\frac{\gamma }{\text{sin}\gamma }$$
$`(4)`$
and

$$_0^{\mathrm{}}𝑑s𝑑tst\text{e}^{\alpha (s^2+t^2+2st\text{cos}\gamma )}=\frac{1\gamma \text{cot}\gamma }{4\alpha ^2\text{sin}^2\gamma }.$$
$`(5)`$
Therefore, one has for a function $`d(u)`$:

$$d_i^{\{k\}}(u)=()^k_0^{\mathrm{}}𝑑\alpha \alpha ^k\text{e}^{\alpha u}\overline{d}_i(\alpha )$$
$`(6)`$
and gets

$$_0^{\mathrm{}}𝑑s𝑑td^{}(u)=\frac{\gamma }{2\text{sin}\gamma }d(0),$$
$`(7a)`$
$$_0^{\mathrm{}}𝑑s𝑑tud^{\prime \prime }(u)=\frac{\gamma }{2\text{sin}\gamma }d(0),$$
$`(7b)`$
$$_0^{\mathrm{}}𝑑s𝑑tstd^{\prime \prime }(u)=\frac{1\gamma \text{cot}\gamma }{4\text{sin}^2\gamma }d(0).$$
$`(7c)`$