text:
# Self-trapping of strong electromagnetic beams in relativistic plasmas

## I Introduction

The problem of electromagnetic (EM) wave propagation and related phenomena in relativistic plasmas has attracted considerable attention in the recent past. From the nonthermal emission of the high-energy radiation coming from a variety of compact astrophysical objects it has become possible to deduce the presence of a population of relativistic electrons in the plasma created in the dense radiation fields of those sources . The principal components of these plasmas could be either relativistic electrons and nonrelativistic ions (protons), or relativistic electron-positron (e-p) pairs.

Relativistic e-p dominated plasmas may be created in a variety of astrophysical situations. Electron-positron pair production cascades are believed to occur in pulsar magnetospheres . The e-p plasmas are also likely to be found in the bipolar outflows (jets) in Active Galactic Nuclei (AGN) , and at the center of our own Galaxy . In AGNs, the observations of superluminal motions are commonly attributed to the expansion of relativistic e-p beams in a pervading subrelativistic medium. This model implies copious pair production via $`\gamma \gamma `$ interactions creating an e-p atmosphere around the source. The actual production of e-p pairs due to photon-photon interactions occurs in the coronas of AGN accretion disks, which upscatter the soft photons emitted by the accretion disks by inverse Compton scattering. The presence of e-p plasma is also argued in the MeV epoch of the early Universe. In the standard cosmological model, temperatures in the MeV range ($`T10^{10}K1MeV`$) prevail up to times $`t=1s`$ after the Big Bang . In this epoch, the main constituent of the Universe is the relativistic e-p plasma in equilibrium with photons, neutrinos, antineutrinos, and a minority population of heavier ions.

Contemporary progress in the development of super strong laser pulses with intensities $`I10^{2123}W/cm^2`$ has also made it possible to create relativistic plasmas in the laboratory by a host of experimental techniques . At the focus of an ultrastrong laser pulses, the electrons can acquire velocities close to the speed of light opening the possibility of simulating in the laboratory the conditions and phenomena that, generally, belong in the astrophysical realm .

Elucidation of the electromagnetic wave dynamics in a relativistic plasmas will, perhaps, be an essential tool for understanding the radiation properties of astrophysical objects as well as of the media exposed to the field of superstrong laser radiation. Although the study of wave propagation in relativistic plasmas has been in vogue for some time, it is only in the recent years that the nonlinear dynamics of EM radiation in e-p dominated plasmas  has come into focus. The enhanced interest stems from two facts : 1) e-p plasmas seem to be essential constituents of the universe, and 2) under certain conditions, even an ultrarelativistic electron-proton plasma can behave akin to an e-p plasma .

Wave self-modulation and soliton-formation is, perhaps, one of the more interesting and significant signatures of the overall plasma dynamics. The existence of stable localized envelop solitons of EM radiation has been suggested as a potential mechanism for the production of micropulses in AGN and pulsars . In the early universe localized solitons are strong candidates to explain the observed inhomogeneities of the visible universe .

In e-p plasmas, there does not exist a general satisfactory theory for the generation of soliton like structures by ultrastrong high-frequency EM fields of arbitrary spatio-temporal shape. Relative complexity of the fully relativistic equations (hydrodynamic or kinetic) has limited their solutions essentially to one-dimensional systems . A summary of the salient results is : in unmagnetized e-p plasmas, high-frequency pressure of the EM pulse pushes e-p pairs in the direction of its propagation thus creating a density hump in the region of the field localization. The effective refractive index of medium decreases in this region and as a result localized soliton-formation is not supported by the medium. (The increase in the refractive index due to the relativistic nonlinearity related with the particles’ high-frequency motion is not enough to cope with the decrease caused by the above mentioned effect). In Refs. it has been argued that localized solitons can be formed if the interaction between the EM field and acoustic phonons is taken into account- the envelope solitons propagating with subsonic velocities may, then, emerge. In magnetized e-p plasma, larger classes of soliton solutions are possible-typical examples may be found in Refs.. However, it is conceivable that soliton solutions obtained in a one dimensional formulation will turn out to be unstable in higher dimensions.

It is, therefore, a matter of utmost priority that we explore the possibility of finding stable multidimensional soliton solutions in e-p plasmas. Dynamics of $`3D`$ envelope solitons of arbitrary strong EM fields in e-p plasma with a small fraction of heavy ions has been analyzed in Ref.. It was shown that, in a transparent e-p plasma, EM pulses with $`L_{||}<<L_{}`$ (where $`L_{||}`$ and $`L_{}`$ respectively,the characteristic longitudinal and transverse spatial dimensions of the field), may propagate as stable, nondiffracting and nondispersing objects (light bullets) with large density bunching. It was further shown in Ref.  that these bullets are exceptionally robust: they can emerge from a large variety of initial field distributions, and are remarkably stable. In these studies, the EM field is pulse-like with longitudinal localization much stronger than the transverse; the localization is brought about by the charge separation electric field (usually absent in a pure e-p plasma) created by the presence of a small population of ions.

In the present paper we explore another mechanism for localization - we will show that the in the pure e-p plasma, the thermal pressure can provide the confining ”glue” just as the charge separation electric field did for an e-p plasma contaminated with ions. We will also deal with a complimentary manifestation of the radiation field- the ”beam” ($`L_{}<<L_{||}`$) rather than the pulse ($`L_{}>>L_{||}`$). Assuming the plasma to be transparent to the beam, and applying a fully relativistic hydrodynamical model, we demonstrate the possibility of beam self-trapping leading to the formation of stable $`2D`$ solitonic structures. The high-frequency pressure force of the EM field ( tending to completely expel the pairs radially from the region of localization) is overwhelmed by the thermal pressure force which opposes the radial expansion of the plasma creating conditions for the formation of the stationary self-guiding regime of beam propagation.

## II Basic Equations

In this paper the word ”relativistic” connotes two distinct regimes: the plasma becomes relativistic when either the directed fluid velocity approaches the speed of light or the thermal energy per particle is of the order of, or larger than the rest mass energy. Since both these paths to relativity are encountered in the astrophysical as well as laboratory plasmas( produced and accelerated by intense laser pulses), we will investigate a fully relativistic hydrodynamical model.

If the velocity distribution of the particles of species $`\alpha `$ ($`=e,p,i,\mathrm{}`$ where $`e`$,$`p`$, and $`i`$ denote respectively electrons, positrons and heavy ions (protons)) is taken to be a local relativistic Maxwellian, the hydrodynamics of such fluids is described by :

$$\frac{}{x_k}(U_\alpha ^iU_{\alpha k}W_\alpha )\frac{}{x_i}P_\alpha =\frac{1}{c}F^{ik}J_{\alpha k}$$
(1)
where $`U_\alpha ^i=[\gamma _\alpha ,\gamma _\alpha 𝐮_\alpha /c]`$ is the hydrodynamic four-velocity with $`𝐮_\alpha `$ as the three-velocity , $`\gamma _\alpha =(1u_\alpha ^2/c^2)^{1/2}`$ is the relativistic factor, $`J_{\alpha k}`$ is the four-current, $`F^{ik}`$ is the electromagnetic field tensor, and $`W_\alpha `$ is the enthalpy per unit volume

$$W_\alpha =\frac{n_\alpha }{\gamma _\alpha }m_{0\alpha }c^2G_\alpha \left(\frac{T_\alpha }{m_{0\alpha }c^2}\right)$$
(2)
Here $`m_{0\alpha }`$ and $`T_\alpha `$ are the particle invariant rest mass and temperature, respectively, $`n_\alpha `$ is the density in the laboratory frame of the fluid of species $`\alpha `$. The pressure $`P_\alpha =n_\alpha T_\alpha /\gamma _\alpha `$, and $`G_\alpha (z_\alpha )=K_3(1/z_\alpha )/K_2(1/z_\alpha )`$, where $`K_2`$ and $`K_3`$ are, respectively, modified Bessel functions of the second and third order and $`z_\alpha =T_\alpha /m_{0\alpha }c^2`$. The factor $`G_\alpha (z_\alpha )`$ has the following asymptotes, $`G_\alpha 1+5z_\alpha /2`$ for $`z_\alpha <<1`$ ( non-relativistic) and $`G_\alpha 4z_\alpha `$ for $`z_\alpha >>1`$ ( highly relativistic).

The set of equations (1)-(2) may be written in the standard form :

$$\frac{d_\alpha }{dt}(m_{0\alpha }c^2G_\alpha \gamma _\alpha )\frac{1}{n_\alpha }\frac{P_\alpha }{t}=e_\alpha 𝐮_\alpha 𝐄$$
(3)
$$\frac{d_\alpha }{dt}(G_\alpha 𝐩_\alpha )+\frac{1}{n_\alpha }P_\alpha =e_\alpha 𝐄+\frac{e_\alpha }{c}(𝐮_\alpha \times 𝐁)$$
(4)
where $`𝐩_\alpha =\gamma _\alpha m_{0\alpha }𝐮_\alpha `$ is the hydrodynamical momentum, $`𝐄`$ and $`𝐁`$ are the electric and magnetic fields, and $`d_\alpha /dt=/t+𝐮_\alpha `$ is the comoving derivative. The hydrodynamical velocity $`𝐮_\alpha `$ and the relativistic $`\gamma _\alpha `$ are related to the momentum by the standard relations: $`𝐮_\alpha =𝐩_\alpha /m_{0\alpha }\gamma _\alpha `$ and $`\gamma _\alpha =(1+𝐩_{\alpha }^{}{}_{}{}^{2}/m_{0\alpha }^2c^2)^{1/2}`$. It is interesting to note from Eqs. (3)-(4) that the fluid inertia is modified by the temperature; the expression $`M_{eff}^\alpha =m_{0\alpha }G_\alpha (z_\alpha )`$ denotes the effective mass of the particle. For ultrarelativistic temperatures ($`T_\alpha >>m_{0\alpha }c^2`$), the effective mass turns out to be $`M_{eff}^\alpha =4T_\alpha /c^2>>m_{0\alpha }`$. Thus the particles ”forget” their rest mass and the plasma turns into a kind of ”photon” gas. If an ultrarelativistic plasma is in thermodynamical equilibrium with the high-frequency photon gas ($`\mathrm{}\omega T`$), one should also take into account the radiation pressure $`P_R=\sigma T^4(\sigma =\pi /45h^3c^3)`$ . In this paper this effect will be neglected. We must also bear in mind that for extremely relativistic regimes, the model Eqs.(3)-(4) fail to adequately describe the plasma dynamics since heavy particle production has been neglected. This shortcoming will impose an upper limit on the temperature for the validity of the model.. Note that in the context of early universe, the epoch in which the e-p plasma is dominant, has a characteristic temperature $`T_\alpha 1MeV`$ and $`M_{eff}4m_{0e}`$. Since the particle masses are just a few times larger than their rest mass at these temperatures, the e-p plasma can still be considered as a two component fluid rather than a photon gas.

The equation of state directly follows from the self-consistency of Eqs.(3) and (4). Taking the scalar product of Eq.(4) with $`𝐮_\alpha `$ and comparing it with Eq.(3), we can derive $`d_\alpha lnP_\alpha /dt=z_\alpha d_\alpha G_\alpha /dt`$. After straightforward manipulation, the equation can be easily integrated to yield

$$P_\alpha =C\frac{K_2(z_\alpha )}{z_\alpha }\mathrm{exp}(z_\alpha G_\alpha )$$
(5)
where the arbitrary constant $`C`$ must be defined by the initial state. Using $`P_\alpha =n_\alpha T_\alpha /\gamma _\alpha `$ Eq.(5) reduces to the adiabatic equation of ”state”:

$$\frac{n_\alpha }{\gamma _\alpha }\frac{z_\alpha }{K_2(z_\alpha )}\mathrm{exp}(G_\alpha z_\alpha )=const$$
(6)
For nonrelativistic temperatures, Eq.(6) yields the usual result for a monoatomic ideal gas ($`n_\alpha ^r/T_\alpha ^{3/2}=const`$, where $`n_\alpha ^r=n_\alpha /\gamma _\alpha `$ is the density in the rest frame of fluid element) and for ultrarelativistic temperatures one obtains the equation of state for the photon gas ($`n_\alpha ^r/T_\alpha ^3=const`$). We would like to emphasize that the $`\gamma _\alpha `$ factor appearing in Eq.(6) is related to the coherent or directed motion of fluid elements whose origin may lie either in an initial macroscopic flow or in the motion imparted by intense EM radiation.

To complete the picture we must couple the plasma equations of motion with Maxwell equations:

$$c\times 𝐁=\frac{𝐄}{t}+4\pi 𝐉$$
(7)
$$c\times 𝐄=\frac{𝐁}{t}$$
(8)
$$𝐄=4\pi \rho $$
(9)
$$𝐁=0$$
(10)
where

$$\rho =\underset{\alpha }{}e_\alpha n_\alpha ,𝐉=\underset{\alpha }{}e_\alpha n_\alpha 𝐮_\alpha $$
(11)
are, respectively, the charge and current densities. The system of Eqs.(3)-(11) along with the continuity equation ( for each species)

$$\frac{n_\alpha }{t}+(n_\alpha 𝐮_\alpha )=0,$$
(12)
represents a closed set of equations which describe propagation of EM radiation in relativistic multicomponent plasmas.

The above system can be manipulated further to reveal interesting structural properties. To begin with, Eqs.(5)-(6) can be cast in the form

$$\frac{1}{n_\alpha }P_\alpha =\frac{m_{0\alpha }c^2}{\gamma _\alpha }G_\alpha $$
(13)
which, when substituted into Eq.(4), converts it to

$$\frac{}{t}(G_\alpha 𝐩_\alpha )+m_{0\alpha }c^2(G_\alpha \gamma _\alpha )=e_\alpha 𝐄+[𝐮_\alpha \times 𝛀_\alpha ]$$
(14)
where

$$𝛀_\alpha =\frac{e_\alpha }{c}𝐁+\times (G_\alpha 𝐩_\alpha )$$
(15)
is the so called generalized vorticity . Taking the curl of Eq.(14), we find that the evolution equation for $`𝛀_\alpha `$

$$\frac{𝛀_\alpha }{t}=\times [𝐮_\alpha \times 𝛀_\alpha ]$$
(16)
is of the standard vortex dynamics form. Although the system of Eqs.(14)-(16) can be traced to early publications (see for instance Ref.), their consequences are yet to be fully worked out. An immediate consequence, for instance, is the appropriate equivalent of Kelvin’s theorem- the flux of generalized vorticity $`\mathrm{\Omega }_\alpha `$ is frozen-in through a comoving area.

The system yields the following set of relativistic Beltrami- Bernoulli equations for equilibrium :

$$𝛀_\alpha =a_\alpha 𝐮_\alpha $$
(17)
$$G_\alpha \gamma _\alpha +\frac{e_\alpha \varphi }{m_{0\alpha }c^2}=const$$
(18)
where $`a_\alpha `$ are constants and $`\varphi `$ is the scalar potential. The relevance of these equilibria for astrophysics is the subject matter of a forthcoming paper.

For the current effort, we apply Eqs.(14)-(16) for wave processes in an unmagnetized plasma. From Eq.(16) it follows that if the generalized vorticity is initially zero ($`𝛀_\alpha =0`$) everywhere in space, it remains zero for all subsequent times. We assume that before the EM radiation is ”switched on”, the generalized vorticity of the system is zero. In this case the equation of motion may be written as:

$$\frac{}{t}𝚷_\alpha +m_{0\alpha }c^2\mathrm{\Gamma }_\alpha =e_\alpha 𝐄$$
(19)
where the temperature dependent momentum $`𝚷_\alpha `$ and $`\mathrm{\Gamma }_\alpha `$ are defined by:

$$𝚷_\alpha =G_\alpha 𝐩_\alpha $$
(20)
$$\mathrm{\Gamma }_\alpha =G_\alpha \gamma _\alpha =\sqrt{G_\alpha ^2+(𝚷_\alpha /m_{0\alpha }c)^2}$$
(21)
The condition of vanishing generalized vorticity connects $`𝚷_\alpha `$ with the magnetic field:

$$𝐁=\frac{c}{e_\alpha }\times 𝚷_\alpha $$
(22)
It is remarkable that in Eq.(19) the magnetic part of the Lorentz force is formally absent; this fact greatly simplifies analytical manipulations. It is equally remarkable that our equations which describe the dynamics of a hot relativistic plasma are structurally similar to equations used in the theoretical treatment of different aspects of ultrastrong laser interaction with a cold plasma . This similarity becomes even more evident when we study the interaction of short EM pulse with relativistic electron-ion plasmas. If the pulse is assumed to be shorter than the characteristic time for ion response (i.e. inverse of ion Langmuir frequency), the ion motion may be ignored, and the electric field may be found from the electron part of Eq.(19),

$$e𝐄=\frac{𝚷_𝐞}{t}m_{0e}c^2\mathrm{\Gamma }_e$$
(23)
Substituting this expression into Poisson’s Eq.(3) (which now reads as $`𝐄=4\pi e(n_{0i}n_e)`$, where $`n_{0i}`$ is the equilibrium ion density) we find the electron density

$$\frac{n_e}{n_{0i}}=1+\frac{1}{m_{0e}\omega _e^2}\frac{}{t}𝚷_𝐞+\frac{c^2}{\omega _e^2}\mathrm{\Delta }\mathrm{\Gamma }_e$$
(24)
where $`\omega _e=(4\pi e^2n_{0i}/m_{0e})^{1/2}`$ is the plasma frequency. Using Eqs.(22)-(24), Eq.(7) reduces to :

$$c^2\times \times 𝚷_𝐞+\frac{^2𝚷_𝐞}{t^2}+m_{0e}c^2\frac{\mathrm{\Gamma }_e}{t}+$$

$$+\omega _e^2\frac{𝚷_𝐞}{\mathrm{\Gamma }_e}\left[1+\frac{1}{m_{0e}\omega _e^2}\frac{}{t}𝚷_𝐞+\frac{c^2}{\omega _e^2}\mathrm{\Delta }\mathrm{\Gamma }_e\right]=0$$
(25)
which, along with the equation of state ($`z_e=m_{0e}c^2/T_e`$)

$$\frac{n_\alpha G_e}{\mathrm{\Gamma }_e}\frac{z_e}{K_2(z_e)}\mathrm{exp}(G_ez_e)=const,$$
(26)
constitutes the simplified system to which the entire set of Maxwell and relativistic hydrodynamic equations for an the electron -ion plasma have been reduced. An equation similar to Eq.(25) was derived in the cold plasma limit in Ref.. However, there are important differences : a) Due to the temperature dependence of $`G_e`$ in Eq.(21), the factor $`\mathrm{\Gamma }_e`$ and the momentum $`\mathrm{\Pi }_e`$ are no more related by simple relations as they do for a ”cold” case, and b) to incorporate the temperature variation in the system we must add the equation of state (26).

Though Eqs.(25)-(26) form a more complicated system than what we have for the cold plasma, we believe that many results obtained in the cold plasma limit can find appropriate analogies in the hot relativistic- temperature case. Detailed studies of the nonlinear dynamics of the electron-ion plasma is beyond the intended scope of the current paper and some new results will be presented separately elsewhere. In the next part of the current paper, we concentrate on the nonlinear dynamics of EM beams in e-p dominated plasmas.

## III The electron-positron dominated plasma

In this section we apply our general formulation to the problem of self-trapping of EM beams in pure e-p plasmas with relativistic temperatures. For notational convenience, we replace the subscripts (e) and (p) by superscripts $`()`$, and $`(+)`$. We assume that the equilibrium state of the plasma is characterized by an overall charge neutrality $`n_{\mathrm{}}^{}=n_{\mathrm{}}^+=n_{\mathrm{}}`$, where $`n_{\mathrm{}}^{}`$ and $`n_{\mathrm{}}^+`$ are the unperturbed number densities of the electrons and positrons in the far region of the EM beam localization. In most mechanisms for creating e-p plasmas, the pairs appear simultaneously and due to the symmetry of the problem it is natural to assume that $`T_{\mathrm{}}^{}=T_{\mathrm{}}^+=T_{\mathrm{}}`$, where $`T_{\mathrm{}}^{}`$ and $`T_{\mathrm{}}^+`$ are the respective equilibrium temperatures .

We shall assume that for the radiation field of interest, the plasma is underdense and transparent, i.e., $`ϵ=\omega _e/\omega 1`$, where $`\omega `$ is the mean frequency of EM radiation and $`\omega _e=(4\pi e^2n_{\mathrm{}}/m_{0e})^{1/2}`$ is the plasma frequency. Since both species are mobile, the e-p dynamics can not be reduced to just one vector equation similar to Eq.(25). We will display the entire set in terms of potentials ( the Coulomb gauge $`𝐀=0`$ will be used),

$$𝐄=\frac{1}{c}\frac{𝐀}{t}\varphi ,𝐁=\times 𝐀,$$
(27)
and the dimensionless quantities $`\stackrel{~}{t}=\omega t`$, $`\stackrel{~}{𝐫}=(\omega /c)𝐫`$, $`\stackrel{~}{T}^\pm =T^\pm /m_{0e}c^2`$, $`\stackrel{~}{𝐀}=e𝐀/(m_{0e}c^2)`$, $`\stackrel{~}{\varphi }=e\varphi /m_{0e}c^2`$, $`\stackrel{~}{𝚷}^\pm =𝚷^\pm /(m_{0e}c)`$, and $`\stackrel{~}{n}^\pm =n^\pm /n_{\mathrm{}}`$. Suppressing the label- ”tilde”, we may arrive at the dimensionless equations,

$$\frac{𝚷^\pm }{t}+\mathrm{\Gamma }^\pm =\frac{𝐀}{t}\varphi $$
(28)
$$\frac{^2𝐀}{t^2}\mathrm{\Delta }𝐀+\frac{}{t}\varphi ϵ^2(𝐉^+𝐉^{})=0$$
(29)
$$\mathrm{\Delta }\varphi =ϵ^2(n^{}n^+)$$
(30)
$$𝐀=0$$
(31)
$$\frac{n^\pm }{t}+𝐉^\pm =0$$
(32)
with $`𝐉^\pm =n^\pm 𝚷^\pm /\mathrm{\Gamma }^\pm `$ and $`\mathrm{\Gamma }^\pm =\sqrt{(G^\pm )^2+(𝚷^\pm )^2}`$. The species equation of state is:

$$\frac{n^\pm }{\mathrm{\Gamma }^\pm f(T^\pm )}=\frac{1}{\mathrm{\Gamma }_{\mathrm{}}^\pm f(T_{\mathrm{}})}$$
(33)
where

$$f(T^\pm )=\frac{K_2(1/T^\pm )T^\pm }{G^\pm }\mathrm{exp}[G^\pm /T^\pm ]$$
(34)
Of various techniques that could be invoked to investigate Eqs.(28)-(34) to study the self-trapping of high-frequency EM radiation propagating along the $`z`$\- axis, we choose the method presented in the excellent paper by Sun et al. . The method is based on the multiple scale expansion of the equations in the small parameter $`ϵ`$. Assuming that all variations are slow compared to the variation in $`\xi =zat`$, we expand all quantities $`Q=(𝐀,\varphi ,𝚷^\pm ,n^\pm ,\mathrm{})`$ as

$$Q=Q_0(\xi ,x_1,y_1,z_2)+ϵQ_1(\xi ,x_1,y_1,z_2)$$
(35)
where $`(x_1,y_1,z_2)=(ϵx,ϵy,ϵ^2z)`$ denote the directions of slow change, and $`a_1=(a^21)/ϵ^21`$. We further assume that the high-frequency EM field is circularly polarized,

$$𝐀_\mathrm{𝟎}=\frac{1}{2}(\widehat{𝐱}+i\widehat{𝐲})A\mathrm{exp}(i\xi /a)+c.c.$$
(36)
Here $`A`$ is the slowly varying envelope of the EM beam, $`\widehat{𝐱}`$ and $`\widehat{𝐲}`$ denote unit vectors, and $`c.c.`$ is the complex conjugate. We now give a short summary of the steps in the standard multiple- scale methodology (Ref.). To the lowest order in $`ϵ`$, we obtain the following. The transverse (to the direction of EM wave propagation $`z`$) component of Eq.(28) reduces to

$$𝚷_\mathrm{𝟎}^\pm =𝐀_\mathrm{𝟎}$$
(37)
and for the longitudinal components we get:

$$a\frac{\mathrm{\Pi }_{0z}^\pm }{\xi }+\frac{\mathrm{\Gamma }_0^\pm }{\xi }=(a)\frac{A_{0z}}{\xi }\frac{\varphi _0}{\xi }.$$
(38)
Equations (29)-(31) yield $`_\xi _{}\varphi _0=_\xi ^2\varphi _0=_\xi A_{z0}=0`$, where $`_{}`$ is the perpendicular Laplacian in $`(x_1,y_1)`$. These relations imply that $`\varphi `$ and $`A_{0z}`$ do not depend on the fast variable $`\xi `$. For the self-trapping problem, we can assume that $`A_{0z}=\mathrm{\Pi }_{0z}^\pm =0`$ . From Eq.(38), and from the lowest order continuity Eq. (32), we obtain : $`_\xi \mathrm{\Gamma }_0^\pm =_\xi n_0^\pm =0`$, i.e., $`\mathrm{\Gamma }_0^\pm `$ and $`n_0^\pm `$ also do not depend on the fast variable $`\xi `$.

To the next order (in $`ϵ`$), the transverse component of Eq.(28) reads :

$$a\frac{𝚷_{\mathrm{𝟏}}^{}{}_{}{}^{\pm }}{\xi }+_{}\mathrm{\Gamma }_0^\pm =(a)\frac{𝐀_\mathrm{𝟏}}{\xi }_{}\varphi _0$$
(39)
Averaging over the fast variable $`\xi `$ we obtain $`_{}\mathrm{\Gamma }_0^\pm =_{}\varphi _0`$ yielding the trivial solution $`\varphi _0=0`$ and $`\mathrm{\Gamma }_0^\pm =\mathrm{\Gamma }_0=const`$. Note that from Eqs.(30) and (33), we can deduce that $`n_0^+=n_0^{}=n_0`$ and $`T_0^+=T_0^{}=T_0`$.

Thus, as one would expect, the low frequency motion of the e-p plasma is driven by the ponderomotive pressure ($`𝚷_{\mathrm{𝟎}}^{}{}_{}{}^{2}`$) of the high-frequency EM field and this force, being same for the electrons and positrons, does not cause charge separation. It is also evident that due to the symmetry between the electron and positron fluids, their temperatures, being initially equal, will remain equal during the evolution of the system. The relation between the EM field and the temperature can be found using the equation $`\mathrm{\Gamma }_0=const`$ obtained above. Using Eqs.(36)-(37) and by choosing the $`const`$ by requiring that at infinity $`A0`$ and $`T_0T_{\mathrm{}}`$, we derive

$$G^2(T_0)=G^2(T_{\mathrm{}})|A|^2$$
(40)
It follows from Eq.(40) that the present hydrodynamical model, which describes the nonlinear waves in e-p plasma, is valid for $`|A|^2/G_{\mathrm{}}^21`$. In our opinion the origin of this restriction lies in the inadequacy of the basic model, and is not solely due to a failure of the perturbation technique used above. When this condition is violated, the EM waves are overturned and they will cause multistream motion of the plasma (i.e. wave breaking takes place). In such a situation, one must resort to kinetic description for studying the nonlinear wave motion. Notice however that the function $`G(T_{\mathrm{}})1`$ if $`T_{\mathrm{}}0`$ but rapidly increases with increase of $`T_{\mathrm{}}`$ thus providing room for $`|A|_{\mathrm{max}}`$ to reach from weak to relativistic values.

We are now ready to deal with the equation for the slowly varying envelope $`A`$ of the EM beam. To the lowest order in $`ϵ`$, one finds from Eq.(29)

$$a_1\frac{^2𝐀_\mathrm{𝟎}}{\xi ^2}_{}^2𝐀_\mathrm{𝟎}2\frac{^2𝐀_\mathrm{𝟎}}{\xi z_2}+2\frac{n_0(T_0)}{G_{\mathrm{}}}𝐀_\mathrm{𝟎}=0$$
(41)
For deriving this equation, we used the relation $`\mathrm{\Gamma }_0=\sqrt{G^2(T_0)+|A|^2}=G_{\mathrm{}}`$, and

$$n_0(T_0)=\frac{f(T_0)}{f(T_{\mathrm{}})}$$
(42)
which follows from Eq.(33). Substituting Eq.(36) into Eq.(41) we find:

$$2i\frac{A}{z}+_{}^2A+\frac{2}{G_{\mathrm{}}}[1n(T_0)]A=0$$
(43)
where subscripts for variables ($`z_2,x_1,y_1`$) are dropped for simplicity. We also assumed without loss of generality that $`(a^21)/ϵ^2a^2=2/G_{\mathrm{}}`$, which in dimensional units coincides with the linear dispersion relation of the EM wave in an e-p plasma, namely: $`\omega ^2=2\omega _e^2/G_{\mathrm{}}+k^2c^2`$ provided that $`a=\omega /kc`$.

Thus, the dynamics of EM beams in hot relativistic e-p plasma has become accessible within the context of a generalized nonlinear Schrödinger equation (NSE) (43).

## IV The self-trapped beams in e-p plasma

In this section we seek the localized 2D soliton solutions of Eq.(43), and analyze the stability of such solutions. Making the self-evident re-normalization of variables $`zzG_{\mathrm{}}`$ , $`r_{}r_{}\sqrt{G_{\mathrm{}}/2}`$, Eq.(43) can be written as:

$$i\frac{A}{z}+_{}^2A+\mathrm{\Psi }A=0$$
(44)
where $`\mathrm{\Psi }=1n_0(T_0)`$ represents the generalized nonlinearity. The companion equation (40) can be viewed as a transcendental algebraic relation between $`T_0`$ and $`|A|^2`$. Thus we conclude that $`\mathrm{\Psi }`$ is a function of $`|A|^2`$ ($`\mathrm{\Psi }=\mathrm{\Psi }(|A|^2)`$. We note that Eq.(44) can be written in the Hamiltonian form $`iA_z=\delta H/\delta A^{}`$, where $`H=𝑑𝐫_{}[|_{}A|^2F(|A|^2)]`$ and $`F(t)=_0^t\mathrm{\Psi }(t^{^{}})𝑑t^{^{}}`$. The Hamiltonian structure implies that Eq.(44) conserves the Hamiltonian $`H`$ in addition to the power (”photon number”) $`N=𝑑𝐫_{}|A|^2`$.

Unfortunately, it is not possible, in general, to derive an explicit analytic relation $`\mathrm{\Psi }=\mathrm{\Psi }(|A|^2)`$ for arbitrary value of $`T_{\mathrm{}}`$. Some qualitative deductions readily follow. Equation (40) shows that the presence of EM radiation reduces the temperature $`T_0`$. Since $`df(T_0)/dT_0>0`$, from Eq.(42) we conclude that the plasma density is also reduced in the region of the EM field localization which is in accordance with adiabatic motion of the plasma. For higher strength of the EM field, a complete expulsion of plasma i.e. plasma cavitation can take place ($`n_00`$); this has been predicted in Ref.. Thus the nonlinearity function $`\mathrm{\Psi }`$ shows a saturating character with the increase of EM field strength (note that present model is valid provided $`|A|^2/(G_{\mathrm{}}^21)1`$). To illustrate, we exhibit in Fig.1 a plot of $`\mathrm{\Psi }`$ versus $`|A|^2`$ for $`T_{\mathrm{}}=0.1`$. One can see that the nonlinearity function indeed saturates at high intensity.

For small temperatures, we can even obtain an analytic expression for the function $`\mathrm{\Psi }`$. Remembering $`T_0T_{\mathrm{}}`$, assuming $`T_{\mathrm{}}1`$, and by using Eq.(42) along with the asymptotic expansions $`G_0(1+5T_0/2)`$ and $`f(T_0^3)`$, we derive for the nonlinearity function:

$$\mathrm{\Psi }=1\left(1\frac{|A|^2}{5T_{\mathrm{}}}\right)^{3/2}$$
(45)
Equations (44)-(45) admit a stationary, nondiffracting axially symmetric solution of the form $`A/\sqrt{5T_{\mathrm{}}}=U(r)\mathrm{exp}(i\lambda z)`$ where $`r=(x^2+y^2)^{1/2}`$ and $`\lambda `$ is the nonlinear wave-vector shift. The radially dependent envelope $`U(r)`$ obeys an ordinary nonlinear differential equation:

$$\frac{d^2U}{dr^2}+\frac{1}{r}\frac{dU}{dr}\lambda U+\mathrm{\Psi }(U^2)U=0$$
(46)
where $`\mathrm{\Psi }=1(1U^2)^{3/2}`$. This equation corresponds to a boundary value problem with the boundary conditions: $`U`$ has its maximum $`U_m`$ at $`r=0`$, and $`U0`$ as $`r\mathrm{}`$. We remind the reader that it has been shown in a seminal paper of Vakhitov and Kolokolov  that such solutions exist for arbitrary saturating nonlinearity functions $`\mathrm{\Psi }`$, provided that the eigenvalue $`\lambda `$ satisfies $`0<\lambda <\mathrm{\Psi }_m`$, where $`\mathrm{\Psi }_m`$ is a maximal value of the nonlinearity function. Equation (46) admits an infinity of discrete bound states characterized by $`j=0,1,2\mathrm{}`$ zeros at finite $`r`$. We consider only the lowest-order nodeless solution of Eq.(46), i. e. ”ground state” that is positive and monotonically decreasing with increasing $`r`$. In the asymptotic region the solution must decay as $`U_r\mathrm{}\mathrm{exp}(\sqrt{\lambda }r)/\sqrt{\lambda r}`$. Our nonlinearity function $`\mathrm{\Psi }`$ has a maximum $`\mathrm{\Psi }_m=1`$ found at $`U_m(=1)`$, i. e. at the maximally allowed strength of the field. As a consequence the upper bound of the propagation constant $`\lambda _c`$ must satisfy $`\lambda _c<\mathrm{\Psi }_m`$. Numerical simulations show that the amplitude of the ground state solution $`U_m=U(r=0,\lambda )`$ is a growing function of $`\lambda `$ (see Fig.2) and it acquires its maximum value ($`=1`$) at $`\lambda =\lambda _c0.29`$.

The solution represents a trapped, localized EM solitary beam. The beam becomes wider at low amplitudes.

The stability of the solitonic solutions can be investigated using the criterion of Vakhitov and Kolokolov - the soliton is stable against small, arbitrary perturbations if $`dN/d\lambda >0`$, where $`N`$ is the power of the trapped mode:

$$N(\lambda )=2\pi _0^{\mathrm{}}𝑑rrU^2(r,\lambda )$$
(47)
In Fig.3 we plot the numerically obtained solutions of $`N`$ for various $`\lambda `$. Since the curve has positive slope everywhere, the corresponding ground state solution is stable for $`0<\lambda <\lambda _c`$. Notice that the power of the solitary beam always exceeds a certain critical value $`N>N_c7.8`$ . We also know that $`N`$ must be bounded from above ($`NN_m10.5`$).

For arbitrary temperatures, explicit form of $`\mathrm{\Psi }=\mathrm{\Psi }(|A|^2)`$ can not be found. However, due to its saturating character, solutions with properties similar to the small temperature case (which can be explicitly solved) could be expected. Using relations (34), (40) and (42), we numerically find a stationary solution of Eq.(44) for arbitrary $`T_{\mathrm{}}`$. For convenience we use following representation of vector potential $`A/A_c=U(r)\mathrm{exp}(i\lambda z)`$, where $`A_c=(G_{\mathrm{}}^21)^{1/2}`$. Though the maximum value of $`U`$ is still restricted by the condition $`0<U_m1`$, the amplitude of vector potential $`A_m`$ can reach a considerable value. For ultrarelativistic temperatures, $`T_{\mathrm{}}>>1`$, we have $`A_c=\sqrt{15T_{\mathrm{}}}>>1`$ and since $`0<A_mA_c`$ the soliton solution with ultrarelativistic strength of EM field is possible. Here we present results of simulations for $`T_{\mathrm{}}=1`$ (i.e. $`T_{\mathrm{}}0.5MeV`$ ).

In Fig.4 we plot the amplitude of the ground state solution $`U_m`$ versus the propagation constant $`\lambda `$. The solution exists provided $`0<\lambda <\lambda _c(0.22)`$. The profiles of the field $`U(r)`$ the plasma density $`n_0(r)`$ and the temperature $`T_0(r)`$ are exhibited in Fig.5 for $`\lambda =0.1`$.

One can see that in the region of field localization, the plasma temperature and density is reduced. Similar plots could be obtained for all allowed values of $`\lambda `$. When $`\lambda \lambda _c`$, plasma cavitation takes place, i.e. at $`r=0`$ the plasma density and temperature tends down to zero. Appearance of zero temperature is not surprising since the corresponding region is the ”plasma vacuum”; all particles are gone away.

The dependence of $`N`$ on $`\lambda `$ is presented in Fig.6. One can see that the curve $`N=N(\lambda )`$ has a positive slope and according to Vakhitov and Kolokolov criterion, the corresponding solitary solutions are stable against small perturbations.

The detailed dynamics of arbitrary field distribution must be studied by direct simulations of Eq.(44). We can learn much , however, from the recent extensive elucidations of the dynamical properties of the solutions of NSE with saturating nonlinearity. It seems that the general features of evolution are not sensitive to the details of the saturating nonlinearity (see for instance  and references therein). For all such systems the beam will monotonically diffract if the beam power is below a critical value ($`N<N_c`$), and it will be trapped if $`N>N_c`$ and the Hamiltonian $`H<0`$ is negative. In the latter case, the beam parameters will oscillate near the equilibrium, ground state values. This oscillations monotonically decrease with increase $`z`$ due to the appearance of the radiation spectrum. For larger $`z`$, the oscillations are damped out, and the formation of the soliton in its ground state takes place. If the initial profile of the beam is close to the equilibrium one, the beam quickly reaches the ground-state equilibrium, and propagates for a long distance without distortion of its shape. The initial beam, even when its parameters (i.e. amplitude, effective width and phase) are quite far from equilibrium, will either focus or defocus to the ground state, exhibiting damped oscillations around it. As a consequence the beam reaches an equilibrium with its final power slightly smaller than the initial. Such an evolutionary scenario may not hold for very intense beams with $`N>>N_c`$; the beam may then break up into filaments due to a modulation instability. However, each filament, will tend to evolve towards its own equilibrium state corresponding to the power it carries. Thus, the ground-state equilibrium seems to be an attractor.

Our own studies indicate that Eq.(44), with the nonlinearity particular to the problem at hand, reproduces the general expected behavior described above. However we find that the soliton formation requires the initial beam power to be in the range $`N_c<N<N_m`$. For $`N>N_m`$, the multistream motion of the plasma prevents the system from settling in a steady state.

## V conclusions

We have investigated the nonlinear propagation of strong EM radiation in a relativistic, unmagnetized two-fluid plasma. The treatment is fully relativistic- in the coherent or directed motion as well as in the random or thermal motion of the plasma particles. The assumption that prior to the switching of the field-plasma interaction, the generalized vorticity is zero, greatly simplifies the system of relativistic fluid equations. In particular, in the electron-ion dominated plasma, under well defined conditions the system of Maxwell -fluid equation ( Eqs.(25)-(26)) turns out to be structurally similar to the one obtained for a cold plasma. Consequently we would expect that results already established in cold plasma limit can find appropriate analogy in the hot plasma case.

We presented a somewhat detailed study of EM beam propagation in transparent e-p plasmas. Applying a reductive perturbation technique, the system of relativistic Maxwell- fluid equations is reduced to a 2D nonlinear Schrödinger equation with a saturating nonlinearity. We found that if the strength of the EM field amplitude is below the wave breaking limit, the beam can enter the self-trapped regime resulting in the formation of stable, self-guided $`2D`$ solitonic structures. The beam-trapping owes its origin to the thermal pressure (which opposes the ponderomotive pressure) - Naturally such trapping regimes are not accessible in the relativistic but cold plasma limit. In the region of beam trapping, the plasma density as well as its temperature is reduced and under certain conditions these parameters can be reduced considerably (i.e. plasma cavitation takes place).

The fact that relativistically hot e-p plasmas are capable of sustaining high amplitude localized structures of high amplitude electromagnetic fields should be a result of considerable importance to an understanding of the complex radiative properties of different astrophysical objects where such plasmas are known to exist.

The work of Z Y is partially supported by Toray Science Foundation, the work of V I B was partially supported by the INTAS Georgian CALL-97, and S M M’s work was supported by the U.S. Department of Energy Contract No. DE-FG03-96ER-54346.


format: text

text:
# Non-Gaussian Correlations in the McLerran-Venugopalan Model

## I Introduction

With the first collisions at the Brookhaven Relativistic Heavy Ion Collider, the experimental frontier in both energy and density has advanced. One of the goals of this machine is to search for the so-called quark-gluon plasma (QGP), a deconfined state of QCD where quarks and gluons form a sea of thermalized particles occupying a volume of space significantly larger than that of a typical hadron. Theoretically, whether and when this occurs depends on the initial conditions of the system. That is, what is the distribution of the quarks and gluons immediately before the nuclei collide?

One particularly fruitful approach to the determination of parton distribution functions in nuclei has been the McLerran-Venugopalan (MV) model. The key observation underpinning the MV model is that for large enough nuclei and for small enough values of the longitudinal momentum fraction $`x__F`$, a new hard scale $`\kappa ^2`$, corresponding to the large color charge per unit transverse area, enters into the problem. Consequently, when $`\kappa ^2\mathrm{\Lambda }_{\mathrm{QCD}}^2`$, it is argued that a classical treatment ought to provide a reasonable approximation to the gluon distribution. Three interrelated factors contribute to making the small-$`x__F`$ region conducive to a classical description. First, a large number of color charges contribute to the source, implying that the vector potential takes on large values within the nucleus. Second, the large vector potential corresponds to the existence of a large number of gluons, inviting us to apply the Weiszäcker-Williams technique. And, third, the new large scale appears in the running strong coupling constant, $`\alpha _s(\kappa ^2)`$, which is small for sufficiently large $`\kappa ^2`$. In the classical treatment advocated in Refs. , quantum mechanical expectation values are replaced by averages over a suitably chosen ensemble of color sources whose statistical weight function, $`W[\rho ]`$, was argued to take on a Gaussian form.

Quantum corrections to the MV model were first computed in Ref. : they were found to be large due to contributions of order $`\alpha _s\mathrm{ln}(1/x__F)`$. The subsequent re-examination of the foundations of the MV model performed in Ref.  led to the formulation of a set of renormalization group equations (RGE) describing the evolution of the weight function $`W[\rho ]`$ as one changes the separation scale between hard and soft partons. Early on in the development of the RGE it was conjectured  that the MV model could be viewed as an effective theory at small $`x__F`$, with the effects of the quantum corrections absorbed into a renormalization of the weight function $`W[\rho ]`$. This conjecture has recently been proven by the calculation in Ref. . While explicit solutions to the evolution equations for $`W[\rho ]`$ have not yet been obtained, it is known that a Gaussian distribution is not a solution to these equations.

While these developments were taking place, we set out to address a different issue, namely the poor infrared behavior of the correlation functions in the MV model. The two-point vector potential obtained in Ref.  grows like $`(x^2)^{x^2}`$ at large distances ($`x\mathrm{\Lambda }_{\mathrm{QCD}}^1`$), signaling the onset of non-perturbative effects associated with confinement. In Ref.  we observed that since individual nucleons do not exhibit a net non-zero color charge, there should not be any long-range ($`\mathrm{\Lambda }_{\mathrm{QCD}}^1`$) correlations between quarks. This requirement of color neutrality was cast into the form of a mathematical constraint on the two-point charge density correlation function. We found that the infrared divergence appearing in the MV model is completely absent when the color neutrality condition is enforced.

With the infrared divergences under control, it became clear that, at least in the classical treatment of Ref. , the $`x__F`$ dependence of the gluon distribution remained trivial, with $`x__FdN/dx__F`$ independent of $`x__F`$. There are two fairly obvious possible sources of non-trivial $`x__F`$ dependence, depending on the value of $`x__F`$ under consideration. First, we can imagine taking $`x__F`$ to be somewhat larger than the values specified in the MV model. In this situation, the longitudinal resolution of the gluons becomes good enough to start to probe the longitudinal structure of the Lorentz-contracted nucleus which they see. This avenue of investigation was pursued in Ref. , where we developed a fully three-dimensional treatment of the classical gluon field of a large nucleus. Our view of the nucleus (radius $`R`$) as containing color-neutral nucleons (each of radius $`a`$) played a vital role in the calculation. Taking advantage of the smallness of the ratio $`a/R`$ for a large nucleus, and assuming that the nucleus possesses a spherically symmetric distribution of color charge in its rest frame, we found that the additional $`x__F`$ dependence manifests itself by the appearance of the combination $`𝒒^2+(x__Fm)^2`$ in some (but not all) of the functions parametrizing the final result (see Eqs. (5.18)–(5.20) of Ref. ). Here $`𝒒`$ is the transverse momentum of the gluon and $`m`$ is the nucleon mass. Over most of the region where our extended treatment is valid, $`(x__Fm)^2𝒒^2`$, implying that these corrections are small.

The second possible source of non-trivial $`x__F`$ dependence, namely quantum corrections proportional to $`\alpha _s\mathrm{ln}(1/x__F)`$, was proposed  even before the first calculation of these corrections in Ref.  verified that they do indeed contain the necessary logarithms. The idea is that a resummation of such powers to all orders would convert the $`1/x__F`$ appearing in the leading-order expression for $`dN/dx__F`$ into $`1/x__F^{1+C\alpha _s}`$ for some $`C`$. One of the goals of the studies conducted in Refs.  may thus be phrased as the determination of $`C`$. A priori it would appear that in order to determine the gluon distribution in the interesting nonlinear regime, we would require a complete solution to the evolution equations for $`W[\rho ]`$ developed in Refs. . And, furthermore, since a Gaussian is not a solution to these equations, it would seem that the general form of the gluon distribution obtained in Refs.  could be fundamentally altered once the non-Gaussian correlations are taken into account. Such an outcome would be at odds with related calculations based on different techniques  which already possess general agreement with the MV model. The purpose of this paper is to investigate this issue. Once again it will prove valuable to incorporate the physics associated with confinement (i.e. the color neutrality of the nucleons) into the discussion. We will demonstrate that to leading order in the strong coupling and $`a/R`$ it is sufficient to consider only the (renormalized) Gaussian contributions to the gluon number density. All contributions from the non-Gaussian portion of $`W[\rho ]`$ are suppressed by additional factors of $`\alpha _s`$ and/or $`a/R`$. Thus, the MV model remains in general agreement with the various other approaches to small-$`x__F`$ physics referred to above, even if we employ a non-Gaussian weight function.

In light of the arguments in favor of a Gaussian form for $`W[\rho ]`$ put forth in Ref. , the reader may be tempted to think that our conclusion is a trivial consequence of the central limit theorem. However, as we will demonstrate, the central limit theorem does not apply to this case. Instead, the additional factors which suppress the non-Gaussian terms relative to their Gaussian counterparts come about because of the fact that the individual nucleons are color neutral objects which do not have significant color-charge correlations at distances much greater than $`a`$.

The remainder of this paper is organized into two main sections. In the first, we begin with a review of the central limit theorem to establish the notation used in the subsequent discussion. We will pay particular attention to the conditions under which the central limit theorem is valid, and then show that for a large nucleus these conditions are violated. The essential observation is that any non-trivial longitudinal correlations, no matter how short range, spoil the Gaussian form of the weight function in the very large nucleus ($`A^{1/3}\mathrm{}`$) limit. In the second half of this paper we will return to the three-dimensional framework for the MV model developed in Ref.  and consider the additional terms which would be generated in the presence of non-Gaussian contributions to $`W[\rho ]`$. Assuming only that the nucleus consists of color-neutral nucleons of radius $`a`$ which possess no nontrivial correlations at separations much larger than $`a`$, we will demonstrate that no contributions beyond those already considered in Refs.  are generated to leading order in $`\alpha _s`$ and $`a/R`$. Thus, for a large enough nucleus, it is sufficient to know only the two-point function associated with $`W[\rho ]`$, even when the quantum corrections are taken into account. We close with a few words about the implications of our result.

## II Gaussian and Non-Gaussian Distributions

A Gaussian (or Normal) distribution in a random variable $`𝒕`$ usually comes about because of the central limit theorem. In this section we review the conditions under which this theorem is valid, and argue that all of these conditions are not, in general, satisfied by the color-charge distribution of a large nucleus, even when we omit radiative corrections. Hence, the statistical weight function $`W[\rho ]`$ appearing in the MV model is generically not expected to be a Gaussian, even classically.

### A Generating Functions

Let $`𝝉`$ be a $`D`$-dimensional random variable, with normalized distribution $`P(𝝉)`$. For what follows it is more convenient to specify it by its moments:

$$\tau _i\tau _j\mathrm{}\tau _kd^D\tau \tau _i\tau _j\mathrm{}\tau _kP(𝝉).$$
(1)
These two descriptions are equivalent because the Fourier transform of $`P(𝝉)`$, namely

$$\stackrel{~}{P}(𝝀)=d^D\tau e^{i𝝀𝝉}P(𝝉),$$
(2)
has a Taylor series expansion in terms of the moments:

$$\stackrel{~}{P}(𝝀)=1+\underset{n=1}{\overset{\mathrm{}}{}}\frac{(i)^n}{n!}(𝝉𝝀)^n.$$
(3)
The moments defined by Eq. (1) are reducible in the following sense. Suppose that the distribution $`P(𝝉)`$ factorizes as

$$P(𝝉)=\underset{i=1}{\overset{D}{}}p(\tau _i),$$
(4)
where $`p(\tau _i)`$ is some single variable distribution, i.e. suppose that the components of the vector $`𝝉`$ are completely unrelated to each other. In this case, the moment $`\tau _i\tau _j\mathrm{}\tau _k`$ factorizes into $`\tau _i\tau _j\mathrm{}\tau _k`$. More generally, however, the components of $`𝝉`$ will not be completely independent of each other. Then, the factorization just described will fail: for example, $`\tau _i\tau _j\tau _i\tau _j`$. The difference between these two quantities provides a measure of the amount of interdependence between the components of $`𝝉`$. This information may be conveniently organized in terms of the (irreducible) cluster moments $`\tau _i\tau _j\mathrm{}\tau _k`$, which are recursively defined through the relations

$`\tau _i`$ $`=`$ $`\tau _i,`$ (6)
$`\tau _i\tau _j`$ $`=`$ $`\tau _i\tau _j+\tau _i\tau _j,`$ (8)
$`\tau _i\tau _j\tau _k`$ $`=`$ $`\tau _i\tau _j\tau _k+\tau _i\tau _j\tau _k+\tau _j\tau _i\tau _k+\tau _k\tau _i\tau _j+\tau _i\tau _j\tau _k,`$ (10)
$`\tau _i\tau _j\tau _k\tau _l`$ $`=`$ $`\tau _i\tau _j\tau _k\tau _l`$ (12)
$`+\tau _i\tau _j\tau _k\tau _l+\tau _j\tau _i\tau _k\tau _l+\tau _k\tau _i\tau _j\tau _l+\tau _l\tau _i\tau _j\tau _k`$ (14)
$`+\tau _i\tau _j\tau _k\tau _l+\tau _i\tau _k\tau _j\tau _l+\tau _i\tau _l\tau _j\tau _k`$
$`+\tau _i\tau _j\tau _k\tau _l+\tau _i\tau _k\tau _j\tau _l+\tau _i\tau _l\tau _j\tau _k`$ (16)
$`+\tau _j\tau _k\tau _i\tau _l+\tau _j\tau _l\tau _i\tau _k+\tau _k\tau _l\tau _i\tau _j`$
$`+\tau _i\tau _j\tau _k\tau _l.`$ (17)
The generalization to higher orders is obvious. It turns out that the cluster moments are generated by $`\mathrm{ln}\stackrel{~}{P}(𝝀)`$, in the same way that the ordinary moments are generated by $`\stackrel{~}{P}(𝝀)`$:

$$\mathrm{ln}\stackrel{~}{P}(𝝀)=\underset{n=1}{\overset{\mathrm{}}{}}\frac{(i)^n}{n!}(𝝉𝝀)^n.$$
(18)
A simple example is provided by the following Gaussian distribution:

$$P(𝝉)=\frac{\mathrm{exp}(𝝉^2)}{\pi ^{D/2}}.$$
(19)
All of the odd moments (reducible and irreducible) of Eq. (19) vanish. Of the irreducible moments, only the two-point function $`\tau _i\tau _j`$ is nonvanishing. All other cluster moments are zero, and the even reducible moments are expressible entirely in terms $`\tau _i\tau _j`$.

### B Central Limit Theorem

Suppose the random variable $`𝝉`$ is independently sampled $`N`$ times, yielding the values $`𝝉^{(p)}(p=1,2,\mathrm{},N)`$. The central limit theorem asserts that when $`N1`$, the mean value of these measurements,

$$𝒕\frac{1}{N}\underset{p=1}{\overset{N}{}}𝝉^{(p)},$$
(20)
obeys a Gaussian distribution with

$$t_i=\tau _i=\tau _i$$
(22)
and

$$t_it_j=\frac{\tau _i\tau _j}{N}.$$
(23)
In other words, the distribution of $`𝒕`$ is given by

$$W(𝒕)=\sqrt{\frac{det\mathrm{\Omega }}{(2\pi )^D}}\mathrm{exp}\left[\frac{1}{2}(t_i\tau _i)\mathrm{\Omega }_{ij}(t_j\tau _j)\right],$$
(24)
where

$$\left(\mathrm{\Omega }^1\right)_{ij}\frac{\tau _i\tau _j}{N}.$$
(25)
The proof of the central limit theorem relies crucially on the independence of samplings: the joint distribution $`Q`$ of the sampled variables $`𝝉^{(p)}`$ must be factorizable, i.e.

$$Q(𝝉^{(1)},𝝉^{(2)},\mathrm{},𝝉^{(N)})=\underset{p=1}{\overset{N}{}}P(𝝉^{(p)}).$$
(26)
In terms of $`Q`$, the probability $`W(𝒕)`$ for finding the mean value $`𝒕`$ is given by

$$W(𝒕)=d^D𝝉^{(1)}d^D𝝉^{(2)}\mathrm{}d^D𝝉^{(N)}Q(𝝉^{(1)},𝝉^{(2)},\mathrm{},𝝉^{(N)})\delta ^D\left(𝒕\underset{p=1}{\overset{N}{}}𝝉^{(p)}/N\right).$$
(27)
Using the distribution given in (26), the moment generating function is simply

$`\stackrel{~}{W}(𝝀)`$ $`=`$ $`{\displaystyle d^D𝒕e^{i𝝀𝒕}W(𝒕)}`$ (28)
$`=`$ $`\left[\stackrel{~}{P}\left(𝝀/N\right)\right]^N.`$ (29)
This leads to the cluster moment generating function

$`\mathrm{ln}\stackrel{~}{W}(𝝀)`$ $`=`$ $`N\mathrm{ln}\stackrel{~}{P}\left(𝝀/N\right)`$ (30)
$`=`$ $`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(i)^n}{n!}}{\displaystyle \frac{(𝝉𝝀)^n}{N^{n1}}},`$ (31)
where we have applied Eq. (18) to obtain the second line. Now consider the $`N1`$ limit. If we approximate (31) by just the first term, we obtain $`W(𝝀)=\mathrm{exp}(i𝝀𝝉)`$, corresponding to $`W(𝒕)=\delta ^D(𝒕𝝉)`$. If instead, we include the first two terms in the infinite series, we get

$$\stackrel{~}{W}(𝝀)=\mathrm{exp}\left(i𝝀𝝉\lambda _i\lambda _j\tau _i\tau _j/2N\right),$$
(32)
corresponding to the Gaussian distribution (24). Deviations from this Gaussian are suppressed by the additional powers of $`1/N`$ present in the higher order terms.

### C Non-Gaussian Distribution

The proof of the central limit theorem depends critically on the assumption of independent sampling. If (26) is violated, then the proof fails. Even a small correlation between successive samplings is sufficient to destroy the conclusion, as the following example illustrates.

Suppose we introduce a “nearest-neighbor” correlation by replacing (26) with

$$Q(𝝉^{(1)},𝝉^{(2)},\mathrm{},𝝉^{(N)})=\left[1+\underset{m=1}{\overset{N1}{}}\frac{f(𝝉^{(m)},𝝉^{(m+1)})}{P(𝝉^{(m)})P(𝝉^{(m+1)})}\right]\underset{p=1}{\overset{N}{}}P(𝝉^{(p)}).$$
(33)
To ensure that (33) is properly normalized, we further assume $`f(𝝉^{(m)},𝝉^{(m+1)})`$ to be symmetric in its arguments, and that

$$d^D𝝉^{(m)}f(𝝉^{(m)},𝝉^{(m+1)})=0.$$
(34)
The moment generating function for the distribution (27) is now

$$\stackrel{~}{W}(𝝀)=\left[\stackrel{~}{P}(𝝀/N)\right]^N\left\{1+\frac{(N1)\stackrel{~}{f}(𝝀)}{\left[\stackrel{~}{P}(𝝀/N)\right]^2}\right\},$$
(35)
where

$$\stackrel{~}{f}(𝝀)d^D𝝉^{(1)}d^D𝝉^{(2)}e^{i𝝀[𝝉^{(1)}+𝝉^{(2)}]}f(𝝉^{(1)},𝝉^{(2)}).$$
(36)
Note that $`\stackrel{~}{f}(\mathrm{𝟎})=0`$. Unlike Eq. (31), the irreducible moments generated by

$$\mathrm{ln}\stackrel{~}{W}(𝝀)=N\mathrm{ln}\stackrel{~}{P}(𝝀/N)+\mathrm{ln}\left\{1+\frac{(N1)\stackrel{~}{f}(𝝀/N)}{\left[\stackrel{~}{P}(𝝀/N)\right]^2}\right\}$$
(37)
do not approach a linear function of $`𝝀`$ as $`N\mathrm{}`$. In fact, we have

$$\mathrm{ln}\stackrel{~}{W}(𝝀)=i𝝀𝝉+\mathrm{ln}\left[1i𝝀\mathit{\varphi }\right]+𝒪(1/N),$$
(38)
where

$$\varphi _ji\frac{}{\lambda _j}\stackrel{~}{f}(𝝀)|_{\lambda =0}.$$
(39)
For $`\mathit{\varphi }0`$, Eq. (38) is nonlinear in $`𝝀`$. This distribution possesses non-zero values for all of its irreducible moments, independent of $`N`$. Including the $`𝒪(1/N)`$ contributions does not change this conclusion. The corresponding distribution $`W(𝒕)`$ is no longer a pure Gaussian: the central limit theorem does not hold in this case. While it is straightforward to determine the (new) form of the distribution function for this particular example, in general this task will not be so easy. Fortunately, all we need to know for what follows is the fact that the resulting distribution is fundamentally non-Gaussian when the successive measurements are (even mildly) correlated.

It is useful to see a bit more about how the correlators change when we relax the assumption of independent measurements. With independent sampling, as in Eq. (26), correlators such as $`\tau _i^{(p)}\tau _j^{(q)}\mathrm{}\tau _k^{(s)}`$ factorize into groups with the same superscripts. For example,

$`\tau _i^{(p)}\tau _j^{(q)}\tau _k^{(s)}\tau _l^{(s)}`$ $`=`$ $`{\displaystyle d^D𝝉^{(1)}\mathrm{}d^D𝝉^{(N)}\tau _i^{(p)}\tau _j^{(q)}\tau _k^{(s)}\tau _l^{(s)}P(𝝉^{(1)})\mathrm{}P(𝝉^{(N)})}`$ (40)
$`=`$ $`\tau _i\tau _j\tau _k\tau _l,`$ (41)
assuming that $`p`$, $`q`$, and $`s`$ all take on distinct values. Of particular interest in connection with the MV model will be the special case

$$\tau _i^{(p)}\tau _j^{(q)}=\delta ^{pq}\tau _i\tau _j+(1\delta ^{pq})\tau _i\tau _j.$$
(42)
In contrast, if the measurements are not independent, such as in (33), then the factorizations obtained in Eqs. (41) and (42) are no longer valid. Additional non-factorizable correlators involving different superscripts will, in general, be present, signaling the breakdown of the central limit theorem.

### D Color Charge Distribution in a Large Nucleus

Let $`\rho ^a(x_{},𝒙)`$ be the color charge density of a large nucleus with color $`a`$ and atomic weight $`A1`$, at the transverse position $`𝒙`$ and the longitudinal position $`x_{}`$ (a precise definition of $`x_{}`$ is given in Eq. (49) below). Let

$$\rho _2^a(𝒙)_{\mathrm{}}^{\mathrm{}}𝑑x_{}\rho ^a(x_{},𝒙)$$
(43)
be the accumulated two-dimensional charge density. On account of Lorentz contraction of the longitudinal size of the relativistic nucleus, this is expected to be the relevant quantity governing the gluon distribution, provided that only valence quarks are included in $`\rho `$, and provided that the longitudinal wavelength of the gluon is much bigger than the longitudinal size of the relativistic nucleus. These are two of the basic assumptions underlying the MV model.

The typical size of $`\rho _2^a(𝒙)`$ for a nucleus is of order $`A^{1/3}`$ times the typical size of $`\rho _2^a(𝒙)`$ for a single nucleon. If we identify $`A^{1/3}`$ with the number of measurements $`N`$ of the previous subsections, and the accumulated charge density $`\rho _2^a(𝒙)/A^{1/3}`$ with the mean value $`t_i`$ (with $`i\{a,𝒙\}`$), then the central limit theorem would force $`W(𝒕)`$ (that is, $`W[\rho _2(𝒙)]`$) to be Gaussian in the limit of large $`A^{1/3}`$, provided that the conditions needed to prove the theorem are obeyed. With this mapping, the superscript $`p`$ in the random variable $`\tau _i^{(p)}`$ corresponds to $`x_{}`$, and $`\tau _i^{(p)}`$ corresponds to $`\rho ^a(x_{},𝒙)`$. A summary of these connections appears in Table I.

As discussed in Sec. II C, a necessary condition for the central limit theorem to hold is that the correlators $`\tau _i^{(p)}\tau _j^{(q)}\mathrm{}\tau _k^{(s)}`$ must factorize into groups with identical superscripts. In particular, consider the translation of Eq. (42) into the MV model. The color neutrality condition tells us that $`\rho ^a(x_{},𝒙)=0`$. Therefore, we are left with

$`\rho ^a(x_{},𝒙)\rho ^b(x_{}^{},𝒙^{})`$ $`=`$ $`\rho ^a(x_{},𝒙)\rho ^b(x_{}^{},𝒙^{})`$ (44)
$`=`$ $`\delta (x_{}x_{}^{})\rho _2^a(𝒙)\rho _2^b(𝒙^{}).`$ (45)
This is effectively the form of correlator argued for in Ref. , and employed extensively in the literature .

However, there are two reasons why it is desirable to go beyond the Gaussian approximation represented by Eq. (45). First, at very small $`x__F`$ the quantum corrections in the MV model become large. If we wish to incorporate these corrections via the RGE analysis of Refs. , we must go beyond the Gaussian approximation, since a Gaussian is not a solution to these equations. Second, if we go to somewhat larger values of $`x__F`$, where the quantum corrections may not be too large, the presence of longitudinal correlations between quark charges within a nucleon and the fact that the longitudinal size of a nucleon is not exactly zero begin to have an effect, even classically. In Ref. , we employ the form

$$\rho ^a(x_{},𝒙)\rho ^b(x_{}^{},𝒙^{})\delta ^{ab}\left[\delta (x_{}x_{}^{})\delta ^2(𝒙𝒙^{})C(x_{}x_{}^{};𝒙𝒙^{})\right].$$
(46)
where $`C(x_{}x_{}^{};𝒙𝒙^{})`$ is a reasonably smooth function parameterizing the mutual correlations between pairs of quarks. This function enforces the color neutrality condition. The nonfactorizability of Eq. (46) tells us that for a nucleus with a non-zero longitudinal thickness we should not expect any single layer to be color neutral by itself. Equivalently, we may view the color neutrality condition as forcing us to have non-trivial longitudinal correlations. In Ref. , all of the higher-order cluster moments of the distribution corresponding to Eq. (46) are taken to vanish, implying that $`W[\rho ^a(x_{};𝒙)]`$ is Gaussian. However, Eq. (45) is violated, and the corresponding 2-dimensional distribution $`W[\rho _2^a(𝒙)]`$ is not a Gaussian, because the central limit theorem does not apply to this computation. Nevertheless, the result obtained in Ref.  remarkably turns out to have exactly the same form as if $`W[\rho _2^a(𝒙)]`$ had been Gaussian: for sufficiently small $`x__F`$, it matches the results of Refs. . In the next section we will see that this outcome is quite general, by considering the consequences of allowing for a non-Gaussian weight function.

## III Effect of Non-Gaussian Contributions

We have just demonstrated that we do not expect the statistical weight function $`W[\rho _2(𝒙)]`$ to have a Gaussian form, even in the purely classical case. On the other hand, a Gaussian has been frequently employed in the literature dealing with the MV model . In this section we will argue that for a large enough nucleus, the results which would be obtained with a non-Gaussian weight function are identical to those which follow from a Gaussian weight function to leading order in $`\alpha _s`$ and $`A^{1/3}`$. That is, the new contributions introduced by the non-Gaussian weight are suppressed by additional factors of $`\alpha _s`$ and/or do not contain as many factors of $`A^{1/3}`$. Our conclusion is a consequence of the color neutrality of the nucleons: all non-trivial correlations are limited to length scales of order $`a\mathrm{\Lambda }_{\mathrm{QCD}}^1`$. Beyond the color neutrality of individual nucleons, we will assume nothing about the form of $`W[\rho (\stackrel{}{x})]`$. Our discussion is framed in terms of the three-dimensional extension of the MV model introduced in Ref. , with the addition of non-zero higher-point cluster moments.

### A Diagrammatic Representation of the Gluon Number Density

Imagine a very large nucleus as viewed in its rest frame. The current corresponding to this situation may be written in the simple form

$$J_r^0=\rho (z_r,𝒙_r);J_r^1=J_r^2=J_r^3=0,$$
(47)
where we employ the subscript “$`r`$” to denote rest frame quantities. In this frame the Yang-Mills equations for the vector potential possess the “obvious” time-independent Coulomb solution. Furthermore, since only $`A_r^00`$, we have $`_0A_r^0=A_r=0`$: the Coulomb solution is the same as the covariant gauge solution in this frame. Thus, we conclude that even when we boost to the lab frame, it is natural to begin with the covariant gauge solution for the vector potential. Indeed, it has been observed that although the expression for the gluon number density is most easily written in terms of the vector potential in the light-cone gauge (see Eq. (57) below), it is nevertheless easiest to work in terms of the covariant gauge expression for the current.

In the lab frame, where the nucleus is moving along the $`+z`$ axis with speed $`\beta `$, we take the source to be of the form

$$J^+=\frac{1}{\epsilon }\rho (x_{},𝒙);J^{}=\frac{\epsilon }{2}J^+;𝑱=\mathrm{𝟎}.$$
(48)
The longitudinal coordinate $`x_{}`$ is defined by

$$x_{}\frac{1}{\epsilon }x^{}\frac{\epsilon }{2}x^+.$$
(49)
The parameter

$$\epsilon \sqrt{\frac{2(1\beta )}{1+\beta }}$$
(50)
measures how close the source is to being exactly on the light cone.<sup>*</sup><sup>*</sup>*We define the light-cone coordinates to be $`x^\pm =x_{}=(x^0\pm x^3)/\sqrt{2}.`$ The transverse coordinates $`x^1`$ and $`x^2`$ form a two-vector which we write in bold-face: $`𝒙`$. Our metric has the signature $`(,+,+,+)`$. Thus, the scalar product in light-cone coordinates reads $`q_\mu x^\mu =q^+x^{}q^{}x^++𝒒𝒙.`$ As explained in Ref. , we prefer to work with natural (order unity) quantities, and keep track of all small and large parameters explicitly through $`\epsilon `$. At the end of the calculation we let $`\epsilon 0`$. Interestingly, in terms of $`x_{}`$, the function describing the source is still spherical: the Lorentz contraction that shrinks $`x^{}`$ towards zero in the lab frame is exactly compensated by the factor $`1/\epsilon `$. Reflecting this fact, we will frequently use the notation $`\stackrel{}{x}(x_{};𝒙)`$.

As explained above, it is natural to first solve the Yang-Mills equations for the vector potential in the covariant gauge, using the covariant gauge source (48), and to subsequently transform the result into the light-cone gauge $`A^+=0`$. The result of this procedure may be written in the form

$`A^j(q_{};𝒙)={\displaystyle \frac{g}{iq_{}}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}(ig^2)^{m1}`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}d^mx_{}^{}{}_{}{}^{}\mathrm{exp}(iq_{}x_{}^{}{}_{1}{}^{}){\displaystyle d^3\stackrel{}{\xi }_1^jG(x_{}^{}{}_{1}{}^{}\xi _{}^{}{}_{1}{}^{};𝒙𝝃_1)}`$ (51)
$`\times `$ $`\left({\displaystyle \underset{k=2}{\overset{m}{}}}{\displaystyle d^3\stackrel{}{\xi }_kG(x_{}^{}{}_{k}{}^{}\xi _{}^{}{}_{k}{}^{};𝒙𝝃_k)}\right)\mathbf{[}\mathbf{[}\rho (\stackrel{}{\xi }_1)\rho (\stackrel{}{\xi }_2)\mathrm{}\rho (\stackrel{}{\xi }_m)\mathbf{]}\mathbf{]},`$ (52)
where $`q_{}`$ is the momentum conjugate to $`x_{}`$, i.e.

$$A^j(q_{};𝒙)_{\mathrm{}}^{\mathrm{}}𝑑x_{}\mathrm{exp}(iq_{}x_{})A^j(x_{},𝒙).$$
(53)
The power counting rules we will present in Sec. III B are simplified by using $`A(q_{};𝒙)`$ instead of $`A(x_{},𝒙)`$.

Although Eq. (52) appears complicated at first glance, it is easily understood in terms of the diagrammatic representation introduced in Ref. . The light-cone gauge vector potential is a non-linear function of the source, containing all possible “powers” of $`\rho `$. Note that the series begins at order $`g`$, and that each additional occurrence of the source adds a factor of $`g^2`$. Fig. 1 illustrates the first few diagrams corresponding to the series. The $`m`$th order diagram consists of $`m`$ copies of the source, represented by the circled crosses. In Eq. (52) they appear as the nested multiple commutator

$$\mathbf{[}\mathbf{[}\rho (\stackrel{}{\xi }_1)\rho (\stackrel{}{\xi }_2)\mathrm{}\rho (\stackrel{}{\xi }_m)\mathbf{]}\mathbf{]}\mathbf{\left[}\mathbf{\right[}\mathbf{[}\mathrm{}\mathbf{[}\rho (\stackrel{}{\xi }_1),\rho (\stackrel{}{\xi }_2)\mathbf{]},\rho (\stackrel{}{\xi }_3)\mathbf{]},\mathrm{}\mathbf{]},\rho (\stackrel{}{\xi }_m)\mathbf{]}.$$
(54)
To each source we attach a propagator (Green’s function) connecting the source point $`(\xi _{}^{}{}_{k}{}^{},𝝃_k)`$ to the point $`(x_{}^{}{}_{k}{}^{},𝒙)`$:

$$G(x_{}^{}{}_{k}{}^{}\xi _{}^{}{}_{k}{}^{};𝒙𝝃_k)=\frac{1}{4\pi }\frac{1}{\sqrt{(x_{}^{}{}_{k}{}^{}\xi _{}^{}{}_{k}{}^{})^2+(𝒙𝝃_k)^2}}.$$
(55)
Here $`𝒙`$ is the transverse position at which we wish to know the vector potential. The index labelling the uppermost propagator on the diagram represents the derivative indicated in Eq. (52). Finally, we have an ordered integration over the longitudinal variables $`x_{}^{}{}_{k}{}^{}`$:

$$_{\mathrm{}}^{\mathrm{}}d^mx_{}^{}{}_{}{}^{}_{\mathrm{}}^{\mathrm{}}𝑑x_{}^{}{}_{1}{}^{}_{\mathrm{}}^{x_{}^{}{}_{1}{}^{}}𝑑x_{}^{}{}_{2}{}^{}\mathrm{}_{\mathrm{}}^{x_{}^{}{}_{m1}{}^{}}𝑑x_{}^{}{}_{m}{}^{}.$$
(56)
Mathematically, this integration results from the gauge transformation to the light-cone gauge from the covariant gauge. Physically this integration corresponds to the final state rescatterings which would be present in a computation of the gluon number density based entirely on the covariant gauge. The ordered integration is represented by the dots on the vertical line. The dots are to slide up and down the entire length of the line without passing one another. Although we have written all of the integrations over an infinite range, in practice the source provides non-zero contributions only over a region of size $`R`$, the nuclear radius.

The next step is to connect the gluon number density to the two-point correlation function for the vector potential:

$`{\displaystyle \frac{dN}{dq_{}d^2𝒒}}`$ $``$ $`{\displaystyle \frac{q_{}}{4\pi ^3}}{\displaystyle d^2𝒙d^2𝒙^{}e^{i𝒒(𝒙𝒙^{})}A_i^a(q_{};𝒙)A_i^a(q_{};𝒙^{})}.`$ (57)
Intuitively, the result in Eq. (57) may be understood by envisioning the expansion of the vector potential in terms of creation and annihilation operators and recognizing that $`AA`$ contains the number operator. The vector potential appearing in Eq. (57) is in the light-cone gauge. This choice reflects the fact that the intuitive picture of the parton model is most transparently realized in the light-cone gauge. Strictly speaking, the angled brackets in Eq. (57) represent a quantum-mechanical expectation value. In the MV model, we make a classical approximation to this quantity by performing an ensemble average with an appropriate weighting function $`W[\rho (\stackrel{}{x})]`$.

Given a specific form for $`W[\rho (\stackrel{}{x})]`$, we may (in principle) evaluate Eq. (57) by inserting two copies of (52), performing the required average term-by-term, and summing the resulting series. Diagrammatically, the quantity $`AA`$ may be represented by drawing all possible pairs of diagrams for a single $`A`$, with the understanding that all possible contractions should be performed (see Figs. 24 for the first three orders in this expansion). For a Gaussian weight function, only pairwise contractions appear. On the other hand, for a non-Gaussian distribution, contractions connecting three or more sources must also be considered. Because the nucleons are neutral, $`\rho (\stackrel{}{x})=0`$. Hence there must be no uncontracted sources left over in any diagram.

### B Power Counting Rules

To proceed further, we must make some reasonable assumptions about the correlation functions corresponding to the moments of $`W[\rho (\stackrel{}{x})]`$. Let us parameterize the $`m`$-point correlation function ($`m>1`$) by

$$\rho (\stackrel{}{x}_1)\rho (\stackrel{}{x}_2)\mathrm{}\rho (\stackrel{}{x}_m)𝒮_m\left(\frac{1}{m}\underset{j=1}{\overset{m}{}}\stackrel{}{x}_j\right)𝒟_m(\stackrel{}{x}_1\stackrel{}{x}_2;\stackrel{}{x}_2\stackrel{}{x}_3;\mathrm{};\stackrel{}{x}_{m1}\stackrel{}{x}_m).$$
(58)
The exact choice made for the $`m1`$ difference coordinates is arbitrary and unimportant to our argument. What does matter is the fact that the source corresponds to a large nucleus of radius $`R`$ which is in turn composed of $`A`$ nucleons, each of radius $`a`$. The nucleons themselves are color-neutral. Because of confinement, we do not expect the field to be correlated at distances much greater than $`a`$: what is happening inside one nucleon is largely independent of what is happening inside of the others. Thus, the function $`𝒟_m`$ ought to be small unless $`|\stackrel{}{x}_i\stackrel{}{x}_j|a`$ for all pairs of points. Furthermore, the center-of-mass coordinate ought to point at a position somewhere inside the nucleus (i.e. it should have a magnitude $`R`$) in order for $`𝒮_m`$ to take on a non-negligible value.

The above physical considerations are sufficient to allow us to determine the order of magnitude of an arbitrary diagram in terms of powers of $`\alpha _s`$ and $`R/a`$ ($`A^{1/3}`$). The powers of the coupling are simple. Recall that each of the diagrams for the gluon number density (57) are formed by gluing together two copies of the expansion for the vector potential (52). Since (52) contains only odd powers of $`g`$, the diagram representing the contribution to (57) with a total of $`j`$ sources contains the factor $`g^{2j2}`$ (or, equivalently $`\alpha _s^{j1}`$).

Now let us consider the integrations which go into the computation of the contributions to $`A_i^a(q_{};𝒙)A_i^a(q_{};𝒙^{})`$. The integrand for a given diagram will contain several propagators plus factors of $`𝒮_m`$ and $`𝒟_m`$, depending upon how the sources are contracted. As described above, $`𝒟_m`$ contains the length scale $`a`$ whereas $`𝒮_m`$ contains the length scale $`R`$. The Green’s functions contain no intrinsic scale of their own. What we would like to know is how many powers of $`R/aA^{1/3}`$ are generated when we perform the required integrations, which range over all possible locations of the sources as well as the entire length of the two vertical lines.

To first approximation, the function $`𝒮_m`$ in Eq. (58) is essentially a constant when we stay well inside the large nucleus since the individual nucleons are identical insofar as the strong interaction is concerned. Thus, the quantity $`A_i^a(q_{};𝒙)A_i^a(q_{};𝒙^{})`$ computed from the expansion given in Eq. (52) is of the form

$$(R/a)^cF(q_{};𝒙𝒙^{})(A^{1/3})^cF(q_{};𝒙𝒙^{})$$
(59)
where $`F`$ is a dimensionless function. For a given diagram, the power $`c`$ is equal to the number of independent clusters formed when the vertical lines are removed. Two clusters are independent if they can be represented by a planar diagram (no crossings). For example, $`c=1`$ for the diagrams in Figs. 2, 3, and 4(f)–4(k), whereas $`c=2`$ for the diagrams in Figs. 4(a)–4(e).

The origin of the factor of $`A^{1/3}`$ for each independent cluster is the divergence which appears in some of the integrations when we take $`𝒮_m=const.`$: these are the integrations which get cut off at the scale $`R`$ instead of the scale $`a`$. To see this, imagine that we have changed the integration variables to a set of sum and difference coordinates. The sum (center-of-mass) coordinates describe the positions of the clusters (whether or not these clusters are independent).Note that transverse position $`(𝒙+𝒙^{})/2`$ of the clusters is not integrated over when computing $`A_i^a(q_{};𝒙)A_i^a(q_{};𝒙^{})`$. Hence, the only center-of-mass integrations present are longitudinal in nature. The difference coordinates describe cluster-cluster separations as well as the internal separations of the components of each cluster. All of the difference integrals are finite and so inherit the scale $`a`$ associated with $`𝒟_m`$. On the other hand, the center-of-mass integrations are insensitive to $`a`$: so long as the cluster is located somewhere inside the nucleus, the integrand can be significant. Hence, these integrations produce factors of $`R`$ instead of $`a`$. In the absence of the ordering of the locations of the points on the vertical lines, then, we would obtain a factor $`(R/a)^k`$ for $`k`$ clusters. What the ordering does is to force clusters which cross each other in nonplanar diagrams to maintain a relative separation which does not greatly exceed $`a`$. Thus, if only $`c`$ clusters are free to move independently of each other, the factor we obtain will be only $`(R/a)^c(A^{1/3})^c`$.

Finally, we note that every diagram will pick up an additional factor $`(R/a)^2A^{2/3}`$ when the Fourier transform indicated in Eq. (57) is performed. The origin is the same as above: the integration over $`𝒙𝒙^{}`$ gets cut off at the scale $`a`$ whereas the integration over $`(𝒙+𝒙^{})/2`$ is cut off at the scale $`R`$. Assembling all of the pieces, we have the result that the diagram containing $`j`$ sources and $`c`$ independent clusters contains the factor $`\alpha _s^{j1}(A^{1/3})^{c+2}`$.

With our power-counting rules in hand, we may now classify the relative contributions of the diagrams in Figs. 24. The lowest order diagram is given in Fig. 2, and contains two sources and a single cluster. Hence, its contribution is of order $`\alpha _sA`$. The third-order diagrams of Fig. 3 have a total of three sources tied together as a single cluster. Their contribution is therefore order $`\alpha _s^2A`$ and is suppressed by a power of $`\alpha _s`$ relative to the leading order diagram. Finally, let us look at the fourth-order diagrams of Fig. 4. According to our power-counting rules, diagrams (a)–(e) all contribute at the $`\alpha _s^3A^{4/3}`$ level (four sources and two independent clusters), whereas diagrams (f)–(k) produce only $`\alpha _s^3A`$. We imagine that we have a very large nucleus such that $`\alpha _s^2A^{1/3}`$ is of order 1 or greater. Thus, we consider diagrams (a)–(e) to be as important as the leading order, whereas diagrams (f)–(k) are subleading, being suppressed by one factor of $`a/RA^{1/3}`$.

The generalization to higher orders is obvious. Clearly, the leading diagrams are all ladders: they contain only two-point contractions which do not cross, thus producing the maximum number of independent clusters for a given number of sources. All other diagrams are suppressed: either they contain extra factors of $`\alpha _s`$ or they have fewer factors of $`A^{1/3}`$. Therefore, the $`A`$ and $`\alpha _s`$ dependence of the leading contributions to the final result must look like

$$x__F\frac{dN}{dx__Fd^2𝒒}=c_0\alpha _sA\left(1+c_1\alpha _s^2A^{1/3}+c_2\alpha _s^4A^{2/3}+\mathrm{}\right).$$
(60)
But this is exactly the form obtained by expanding the expressions contained in Eqs. (5.18)–(5.20) of Ref. . In fact, the set of diagrams described above is precisely the set of diagrams considered in Ref. . Moreover, it is also the set of diagrams whose two-dimensional reduction corresponds to those terms which were retained in Refs. . These contributions are easily resummed into a relatively simple “exponential”. Thus, we see that even if the weight function $`W[\rho (\stackrel{}{x})]`$ is non-Gaussian, we obtain the same result to leading order as if we had chosen to use a Gaussian weight function instead.

## IV Implications

We have just demonstrated that the non-Gaussian portions of the statistical weight function $`W[\rho (\stackrel{}{x})]`$ in the MV model do not contribute to the gluon number density at leading order. On the other hand, the results of Ref.  suggest that the quantum corrections to the MV model may be incorporated by solving the RGE for the weight function and using this (non-Gaussian) result to perform the ensemble averages required to compute the gluon number density. A synthesis of these two conclusions has the following implication: it is sufficient to solve the RGE for the new value of the two point charge-density correlation function, and to use this renormalized correlator as defining an effective Gaussian weight function to be input to the MV model. The separation scale between hard and soft partons appearing in Ref.  may be recast as a dependence of the renormalized correlator on $`x__F`$. In terms of Eq. (46), this dependence can potentially show up in the detailed shape of the smooth function $`C(\stackrel{}{x}\stackrel{}{x}^{})`$, as well as in the prefactors which were suppressed in writing down Eq. (46). In terms of the expansion (60), this means that the $`c_i`$’s would be $`x__F`$-dependent. Since the quantum corrections presumably incorporate the physics of the Yukawa cloud believed to surround the nucleons, it is possible that the detailed relation between $`a`$ and $`R`$ could change somewhat, again in an $`x__F`$-dependent way.

In Ref.  it was argued that not only should the RGE-improved MV model apply to nuclei at sufficiently small values of $`x__F`$, but that it should also apply to hadrons as well, at even smaller values of $`x__F`$. Our argument does not apply to this situation. For a hadron we would effectively have $`a/R=1`$, meaning that all of the suppressed contributions we have dropped are no longer unimportant. In this situation the full non-Gaussian solution to the RGE has relevance and must be taken into account.

###### Acknowledgements.

High energy physics research at McGill University is supported in part by the Natural Sciences and Engineering Research Council of Canada and the Fonds pour la Formation de Chercheurs et l’Aide à la Recherche of Québec. High energy physics research at the University of Arizona is supported by the U.S. Department of Energy under contract number DE-FG02-95ER40906.