Quantum chromodynamics at high energyand statistical physics

# Quantum chromodynamics at high energy and statistical physics

S. Munier Centre de physique théorique, École Polytechnique, CNRS, Palaiseau, France
###### Abstract

When hadrons scatter at high energies, strong color fields, whose dynamics is described by quantum chromodynamics (QCD), are generated at the interaction point. If one represents these fields in terms of partons (quarks and gluons), the average number densities of the latter saturate at ultrahigh energies. At that point, nonlinear effects become predominant in the dynamical equations. The hadronic states that one gets in this regime of QCD are generically called “color glass condensates”.

Our understanding of scattering in QCD has benefited from recent progress in statistical and mathematical physics. The evolution of hadronic scattering amplitudes at fixed impact parameter in the regime where nonlinear parton saturation effects become sizable was shown to be similar to the time evolution of a system of classical particles undergoing reaction-diffusion processes. The dynamics of such a system is essentially governed by equations in the universality class of the stochastic Fisher-Kolmogorov-Petrovsky-Piscounov equation, which is a stochastic nonlinear partial differential equation. Realizations of that kind of equations (that is, “events” in a particle physics language) have the form of noisy traveling waves. Universal properties of the latter can be taken over to scattering amplitudes in QCD.

This review provides an introduction to the basic methods of statistical physics useful in QCD, and summarizes the correspondence between these two fields and its theoretical and phenomenological implications.

###### keywords:
Quantum chromodynamics, color dipole model, color glass condensate, stochastic fronts, traveling waves, reaction-diffusion 13.60.Hb, 12.38.-t

## 1 Introduction to high energy scattering in QCD

The study of quantum chromodynamics in the high-energy regime has taken a new soar in the last 15 years with the wealth of experimental data that have been collected, first at the electron-proton collider DESY-HERA, and then at the heavy-ion collider RHIC. More energy in the collision enables the production of objects of higher mass in the final state, and thus the discovery of new particles. But higher energies make it also possible to observe more quantum fluctuations of the incoming objects, that is to say, to study more deeply the structure of the vacuum.

The well-established microscopic theory which describes the interactions of hadronic objects is quantum chromodynamics (QCD). (For a comprehensive textbook, see Ref. [1]). There are not many known analytical approaches to QCD, except perturbative expansions of observables in powers of the strong coupling constant which, thanks to asymptotic freedom, is justified for carefully chosen observables in special kinematical regimes. But fixed-order calculations in QCD are known to usually have a very limited range of applicability. This is because in the evaluation of Feynman graphs, the coupling constant always comes with “infrared” and “collinear” logarithms that are related to the phase space that is available to the reaction, that is to say, to kinematics. Resumming part of these logarithms is mandatory. All of them is too difficult. The question is to carefully select the dominant ones, and this is not at all easy.

At the HERA collider, electrons or positrons scattered off protons at the center-of-mass energy , exchanging a photon of virtuality . Through the scattering, one could probe partonic fluctuations of the proton (made of quarks and gluons) of transverse momenta , and longitudinal momentum fractions .

For a long time, the dominant paradigm had been that the collinear logarithms , that become large when is large compared to the QCD confinement scale , were the most important ones. As a matter of fact, searches for new particles or for exotic physics require to scrutinize matter at very small distances, and hence very large have to be considered. Perturbative series of powers of have to be fully resummed. The equation that performs this resummation is the celebrated Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equation [2, 3, 4].

However, once HERA had revealed its ability to get extremely good statistics in a regime in which is moderate (from to ) and very small (down to ) it became clear that infrared logarithms () could show up and even dominate the measured observables. The resummation of the series of infrared logs is performed by the Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation [5, 6, 7]. The series (with appropriate coefficients) is the leading order (LO), while the series is the next-to-leading order (NLO), which has also been computed [8, 9]. The BFKL equation is a linear integro-differential equation.

At ultrahigh energy, the bare BFKL equation seems to violate the Froissart bound, that states that total hadronic cross sections cannot rise faster than . The latter is a consequence of the unitarity of the probability of scattering. The BFKL equation predicts a power rise with the energy of the form , where is positive and quite large ( to according to the effective value of that is chosen). The point at which the BFKL equation breaks down depends on the value of the typical transverse momentum which characterizes the observable (It is the photon virtuality in the case of deep-inelastic scattering). One may define the energy-dependent saturation scale in such a way that the BFKL equation holds for . For , the probability for scattering to take place is of order 1, and for , it would be larger than 1 if one trusted the BFKL equation. The saturation scale is a central observable, which we shall keep discussing in this review: It signs the point at which the linear (BFKL) formalism has to be corrected for nonlinear effects. The regime in which nonlinearities manifest themselves is a regime of strong color fields, sometimes called the color glass condensate (For the etymology of this term, see e.g. the lectures of Ref. [10]; for a review, see Ref. [11]).

The fact that unitarity is violated is not only due to the lack of a hadronic scale in the BFKL equation, which is a perturbative equation; Introducing confinement in the form of a cutoff would not help this particular problem. It simply means that still higher orders are needed. The NLO corrections to the BFKL kernel indeed correct this behavior in such a way that the description of the HERA data in the small- regime is possible by the BFKL equation. However, these corrections are not enough to tame the power-like growth of cross sections as predicted by the LO BFKL equation. It seems that a resummation of contributions of arbitrary order would be needed.

New equations were proposed well before the advent of colliders able to reach this regime. Gribov-Levin-Ryskin wrote down a model for the evolution of the hadronic scattering cross sections in the early 80’s [12, 13], and Mueller and Qiu derived a similar equation from QCD a bit later [14]. These equations are integral evolution equations with a nonlinear term, which basically takes into account parton saturation effects, that is to say, recombination or rescattering. The latter cannot be described in a linear framework such as the BFKL formalism. Subsequently, more involved QCD evolution equations were derived from different points of view. In the 90’s, McLerran and Venugopalan [15, 16, 17] proposed a first model, mainly designed to approach heavy-ion collisions. Subsequently, Balitsky [18], Jalilian-Marian, Iancu, McLerran, Weigert, Leonidov and Kovner (B-JIMWLK) [19, 20, 21, 22, 23] worked out QCD corrections to this model, and got equations that reduce to the BFKL equation in the appropriate limit. Technically, these equations actually have the form of an infinite hierarchy of coupled integro-differential equations (in Balitsky’s formulation [18]), of a functional renormalizaton group equation, or alternatively, of a Langevin equation (in Weigert’s formulation [23]). A much simpler equation was derived in 1996 by Balitsky [18] and rederived by Kovchegov in 1999 [24, 25] in a very elegant way within a different formalism. The obtained equation is called the Balitsky-Kovchegov equation (BK). The latter derivation was based on Mueller’s color dipole model [26], which proves particularly suited to represent QCD in the high energy limit.

The exciting feature of this kinematical regime of hadronic interaction from a theoretical point of view is that the color fields are strong, although, at sufficiently high energies, the QCD coupling is weak, authorizing a perturbative approach, and thus analytical calculations. In such strong field regime, nonlinear effects become crucial. But the conditions of applicability of the different equations that had been found had never been quite clear. Anyway, these equations are extremely difficult to solve, which had probably been the main obstacle to more rapid theoretical developments in the field until recently.

Furthermore, for a long time, the phenomenological need for such a sophisticated formalism was not obvious, since linear evolution equations such as the DGLAP equation were able to account for almost all data. But Golec-Biernat and Wüsthoff showed that unitarization effects may have already been seen at HERA [27, 28]. Their model predicted, in particular, that the virtual photon-proton cross section should only depend on one single variable , made of a combination of the transverse momentum scale (fixed by the virtuality of the photon ) and . This phenomenon was called “geometric scaling” [29]. It was found in the HERA data (see Fig. 1): This is maybe one of the most spectacular experimental result from HERA in the small- regime.

This observation has triggered many phenomenological and theoretical works. Soon after its discovery in the data, geometric scaling was shown to be a solution of the Balitsky-Kovchegov (BK) equation, essentially numerically, with some analytical arguments (see e.g. [31, 32, 33, 34]). The energy dependence of the saturation scale was eventually precisely computed by Mueller and Triantafyllopoulos [35]. Later, it was shown that the BK equation is actually in the universality class of the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation [36, 37], and geometric scaling was found to be implied by the fact that the latter equation admits traveling wave solutions [38].

A first step beyond the BK equation, in the direction of a full solution to high energy QCD, was taken by Mueller and Shoshi in 2004 [39]. Actually, they did not solve the B-JIMWLK equations, but instead, they solved the linear BFKL equation with two absorptive boundary conditions, which they argued to be appropriate to represent the expected nonlinearities. Geometric scaling violations were found from their calculation, which should show up at any energies.

Subsequently, it was shown that high-energy QCD at fixed coupling is actually in the universality class of reaction-diffusion processes, studied in statistical physics, whose dynamics may be encoded in equations similar to the stochastic FKPP equation [40]. The Mueller-Shoshi solution was shown to be consistent with solutions to the latter equation. So high-energy QCD seems to be in correspondence with disordered systems studied in statistical physics. This correspondence has provided a new understanding of QCD in the high-energy regime, and it has proven very useful to find more features of high-energy scattering. The obtained results go beyond a solution to the B-JIMWLK equation, which in fact, thanks to the new picture, is seen to be incomplete.

Scope

The goal of this review article is to summarize the main ideas behind this conjectured correspondence between scattering at high-energy in QCD and some processes studied in statistical physics, as well as to introduce the QCD reader to the useful technical tools borrowed from statistical physics. We also feel that there is a cultural gap to be filled between statistical physics and particle physics. Indeed, statistical physicists are used to build simple toy models which contain the interesting physics, and whose main properties are likely to be independent of the details of the model, i.e. universal. In QCD, since the theory is well-established, we are often reluctant to give up some of its features to work out exact results in a toy model. One of our aims is to convince the reader that such a way of thinking is efficient in the case of high-energy QCD, by exhibiting results for QCD scattering amplitudes, obtained by looking for the universality class of the considered process, and that are believed to be exact.

Over the last few years, several hundreds of papers have appeared related to this subject, mainly issued from a very active though restricted community. Obviously, we cannot give a complete account of this abundant literature. As a matter of fact, some important recent developments had to be left out, for which we shall only provide references for the interested reader who might want to deepen his study in these directions. Concerning the correspondence itself, we do not attempt to establish a definite stochastic nonlinear evolution equation for QCD amplitudes, for to our judgement, this research line is not mature enough yet: A better understanding of the very saturation mechanism at work in QCD is definitely needed before one may come to this issue. Furthermore, it is not clear to us that a stochastic formulation would be a technical progress, since there are not many known methods to handle complicated stochastic equations. We feel that the same is true for the search for effective actions that would include so-called Pomeron loops. We also do not address the developments based on the boost-invariance symmetry that scattering amplitudes should have: This would drive us too far off the main focus of this review. As for more phenomenological aspects, we only discuss the basic features of total cross sections without attempting to address other observables such as diffraction. We do also not address the issue of next-to-leading effects such as the running of the QCD coupling. This discussion, though crucial if one wants to make predictions for actual colliders, would probably only be technical in its nature: There is no conceptual difference between the fixed coupling and the running coupling case. Here, only basic phenomenological facts brought about by this new understanding of high-energy QCD are addressed, namely geometric scaling and diffusive scaling.

Outline

The outline goes as follows. The next section is devoted to describing scattering in QCD from a -channel point of view, relying essentially on the parton model or rather on an interpretation useful in the high-energy limit, the color dipole model. Once this picture is introduced, it is not difficult to understand the correspondence with reaction-diffusion processes occuring in one spatial dimension, whose dynamics is captured by equations in the university class of the Fisher-Kolmogorov-Petrovsky-Piscounov (FKPP) equation. We then explain how traveling waves appear in this context. In Sec. 3, we study in greater detail a toy model for which many technics (field theory, statistical methods) may be worked out completely. This model however ignores spatial dimensions, and thus, does not account for traveling waves. We summarize the state-of-the-art research on equations in the universality class of the FKPP equation in Sec. 4. Finally, we come back to QCD, discussing the relevance of one-dimensional-like models in the FKPP class, and showing how noisy traveling waves may show up in the actual data.

## 2 Hadronic interactions in a s-channel picture and analogy with reaction-diffusion processes

In this section, we shall introduce the physical picture of high-energy scattering in the parton model. In particular, the color dipole model [26] is described since it is particularly suited to address high-energy scattering, especially close to the regime in which nonlinear effects are expected to play a significant role. In a second part, we shall argue that high-energy scattering is a peculiar reaction-diffusion process.

### 2.1 Parton model and dipoles

#### 2.1.1 General picture

For definiteness, let us consider the scattering of a hadronic probe off a given target, in the restframe of the probe and at a fixed impact parameter, that is to say, at a fixed distance between the probe and the center of the target in the two-dimensional plane transverse to the collision axis. In the parton model, the target interacts through one of its quantum fluctuations, made of a high occupancy Fock state if the energy of the reaction is sufficiently high (see Fig. 2a). As will be understood below, the probe effectively “counts” the partons in the current Fock state of the target whose transverse momenta (or sizes ) are of the order of the one that characterizes the probe: Roughly speaking, the amplitude for the scattering off this particular partonic configuration is proportional to the number of such partons.

The observable that is maybe the most sensitive to quantum fluctuations of a hadron is the cross section for the interaction of a virtual photon with a hadronic target such as a proton or a nucleus. The virtual photon is emitted by an electron (or a positron). What is interesting with this process, called “deep-inelastic scattering”, is that the kinematics of the photon is fully controlled by the measurement of the scattered electron. The photon can be considered a hadronic object since it interacts through its fluctuations into a quark-antiquark state. The latter form color dipoles since although both the quark and the antiquark carry color charge, the overall object is color neutral due to the color neutrality of the photon. The probability distribution of these fluctuations may be computed in quantum electrodynamics (QED). Subsequently, the dipole interacts with the target by exchanging gluons. The dipole-target cross section factorizes at high energy. One typical event is depicted in Fig. 2a.

Dipole models [41, 42] have become more and more popular among phenomenologists since knowing the dipole cross section enables one to compute different kinds of observables. Like parton densities, the latter is a universal quantity, that may be extracted from one process and used to predict other observables. Different phenomenological models may be tried for the dipole cross section. QCD evolution equations may even be derived, as we shall discover below. An accurate recent study of the foundations of dipole models may be found in Ref. [43, 44].

In QCD, the state of a hadronic object, encoded in a set of wave functions, is built up from successive splittings of partons starting from the valence structure. This is visible in the example of Fig. 2a: The quark and the antiquark that build up, in this example, the target in its asymptotic state each emit a gluon, which themselves emit, later on in the evolution, other gluons. As one increases the rapidity by boosting the target, the opening of the phase space for parton splittings makes the probability for high occupation numbers larger. Indeed, the probability to find a gluon that carries a fraction (up to ) of the momentum of its parent parton (which may be a quark or a gluon) is of order for small . There is a logarithmic singularity in , meaning that emissions of very soft gluons (small ) are favored if they are allowed by the kinematics. The splitting probability is of order 1 when the total rapidity of the scattering is increased by roughly , where the convenient notation has been introduced. Only splittings of a quark or of a gluon into a gluon exhibit the singularity. Therefore, at large rapidities, gluons eventually dominate the partonic content of the hadrons.

The parton model in its basic form, where the fundamental objects of the theory (quarks and gluons) are directly considered, is not so easy to handle in the high-energy regime. One may considerably simplify the problem by going to the limit of large number of colors (), in which a gluon may be seen as a zero-size quark-antiquark pair. Then, color-neutral objects become collections of color dipoles, whose endpoints consist in “half gluons” (see Fig. 2b). There is only one type of objects in the theory, dipoles, which simplifies very much the picture. Furthermore, going to transverse coordinate space (instead of momentum space, usually used in the DGLAP formalism) by trading the transverse momenta of the gluons for the sizes of the dipoles (through an appropriate Fourier transform) brings another considerable simplification. Indeed, the splittings that contribute to amplitudes in the high-energy limit are the soft ones, for which the emitted gluons take only a small fraction of the momentum of their parent (the latter being very large). Therefore, the positions of the gluons (and thus of the edges of the dipoles) in the plane transverse to the collision axis are not modified by subsequent evolution once the gluons have been created. Thus, the evolution of each dipole proceeds through completely independent splittings to new dipoles.

We will now see how this picture translates into a QCD evolution equation for scattering amplitudes, first in the regime in which there are no nonlinear effects. In a second step, we will try and understand how to incorporate the latter.

#### 2.1.2 BFKL equation from the dipole model

The building up of the states of each hadron is specified by providing the splitting rate of a dipole whose endpoints have transverse coordinates into two dipoles and as the result of a gluon emission at position . It is computed in perturbative QCD and reads [26]

 dPd(¯αy)((x0,x1)→(x0,x2),(x2,x1))=|x0−x1|2|x0−x2|2|x1−x2|2d2x22π. (1)

Dipole splittings are independent. After some rapidity evolution starting from a primordial dipole, one gets a chain of dipoles such as the one depicted in Fig. 3.

The elementary scattering amplitude for one projectile dipole off a target dipole is independent of the rapidity and reads [26]

 Tel((x0,x1),(z0,z1))=π2α2s2log2|x0−z1|2|x1−z0|2|x0−z0|2|x1−z1|2. (2)

If the target is an evolved state at rapidity , then it consists instead in a distribution of dipoles. The (forward elastic) scattering amplitude is just given by the convolution of and , namely

 A(y,(x0,x1))=∫d2z02πd2z12πT% el((x0,x1),(z0,z1))n(y,(z0,z1)). (3)

Let us examine the properties of . To this aim, it is useful to decompose the coordinates of the dipoles in their size vector (resp. ) and impact parameter (resp. ). In the limit in which the relative impact parameters of the dipoles is very large compared to their sizes, we get the simplified expression

 Tel(ra,rb,b)∼|ra|,|rb|≪|b|α2sr2ar2bb4, (4)

and thus the scattering amplitude decays fast as a function of the relative impact parameter. If instead the relative impact parameter is small (of the order of the size of the smallest dipole), we get

 Tel(ra,rb,b)∼|ra|,|rb|∼|b|α2sr2, (5)

where , , and the integration over the angles has been performed.

Equation (4) means that the dipole interaction is local in impact parameter: It vanishes as soon as the relative distance of the dipoles is a few steps in units of their size. Eq. (5) shows that only dipoles whose sizes are of the same order of magnitude interact. These properties are natural in quantum mechanics. Thus the amplitude in Eq. (3) effectively “counts” the dipoles of size of the order of at the impact parameter (up to ), with a weight factor .

An evolution equation for the amplitude with the rapidity of the scattering can be established. It is enough to know how the dipole density in the target evolves when rapidity is increased, since all the rapidity dependence is contained in in the factorization (3), and such an equation may easily be worked out with the help of the splitting rate distribution (1). It reads [26]

 ∂n(y,(x0,x1))∂(¯αy)=∫d2x22π|x01|2|x02|2|x12|2[n(y,(x0,x2))+n(y,(x2,x1))−n(y,(x0,x1))], (6)

where . The very same equation holds for . The elementary scattering amplitude only appears in the initial condition at , which is not shown in Eq. (6). In a nutshell, the integral kernel encodes the branching diffusion of the dipoles. The total number of dipoles at a given impact parameter grows exponentially, and their sizes diffuse. The appropriate variable in which diffusion takes place is . (This is due to the collinear singularities in Eq. (1).) This equation is nothing but the BFKL equation. A complete solution to this equation, including the impact-parameter dependence, is known [45].

An important property of the amplitude is that it is boost-invariant. This property is preserved in the BFKL formulation. We could have put the evolution in the projectile instead of the target, or shared it between the projectile and the target: The result for the scattering amplitude would have been the same. In a frame in which the target carries units of rapidity and the projectile , the amplitude reads

 A(y,(x0,x1))=∫d2z02πd2z12πd2z′02πd2z′12πn% projectile(y−y′,(z0,z1)|(x0,x1))×Tel((z0,z1),(z′0,z′1))ntarget(y′,(z′0,z′1)). (7)

is the density of dipoles found in a dipole of initial size after evolution over steps in rapidity. If , one recovers Eq. (3). If , then all the evolution is in the projectile instead.

The amplitude is related to an interaction probability, and thus, it must be bounded: In appropriate normalizations, has to range between 0 and 1. But as stated above, the BFKL equation predicts an exponential rise of with the rapidity for any dipole size, which at large rapidities eventually violates unitarity. Hence the BFKL equation does not provide a complete account of high-energy scattering in QCD.

#### 2.1.3 Unitarity and the Balitsky-Kovchegov equation

It is clear that one important ingredient that has been left out in the derivation of the BFKL equation is the possibility of multiple scatterings between the probe and the target. Several among the dipoles in Eq. (7) may actually interact with the dipoles in the other hadron simultaneously. The only reason why such interactions may not take place is that (see Eq. (2)), and thus the probability for two simultaneous scatterings is of order , which is suppressed. But this argument holds only as long as the dipole number densities are of order . If (which is also the point above which the unitarity of is no longer preserved in the BFKL approach), then it is clear that multiple scatterings should occur.

In order to try and implement these multiple scatterings, we introduce the probability that there be no scattering between a dipole and a given realization of the target with total rapidity , that we shall denote by . Let us start with a system in which the evolution is fully contained in the target. We increase the total rapidity by boosting the projectile (initially at rest) by a small amount . Then there are two cases to distinguish, depending on whether the dipole splits in the rapidity interval . In case it splits into two dipoles and , the probability that the projectile does not interact is just the product of the probabilities that each of these new dipoles do not interact. This is because once created, dipoles are supposed to be independent. In summary:

 S(y+dy,(x0,x1))=⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩S(y,(x0,x1))        with proba 1−¯αdy∫x2dPd(¯αy)(x01→x02,x12)S(y,(x0,x2))×S(y,(x2,x1))        with proba ¯αdydPd(¯αy)(x01→x02,x12) (8)

Taking the average over the realizations of the target and the limit , we get

 ∂∂y⟨S(y,(x0,x1))⟩=¯α∫d2x22πx201x202x221[⟨S(y,(x0,x2))S(y,(x2,x1))⟩−⟨S(y,(x0,x1))⟩] (9)

(See Fig. 4 for a graphical representation.)

We see that this equation is not closed: An evolution equation for the correlator is required. However, we may assume that these correlators factorize

 ⟨S(y,(x0,x2))S(y,(x2,x1))⟩=⟨S(y,(x0,x2))⟩⟨S(y,(x2,x1))⟩. (10)

This assumption is justified if the dipoles scatter off independent targets, for example, off the nucleons of a very large nucleus. Writing , we get the following closed equation for :

 ∂∂yA(y,(x0,x1))=¯α∫d2x22πx201x202x221[A(y,(x0,x2))+A(y,(x2,x1))−A(y,(x0,x1))−A(y,(x0,x2))A(y,(x2,x1))], (11)

which is the Balitsky-Kovchegov (BK) equation [24, 25]. Note that if one neglects the nonlinear term, one gets back the BFKL equation (6) (written for instead of ). A graphical representation of this equation is given in Fig. 5.

It is not difficult to see analytically that the BK equation preserves the unitarity of : When becomes of the order of 1, then the nonlinear term gets comparable to the linear terms in magnitude, and slows down the evolution of with , which otherwise would be exponential.

Let us go back to Eqs. (8),(9) and instead of assuming the factorization of the correlators (10), work out an equation for the two-point correlator . From the same calculation as before, we get

 ∂∂y⟨S02S2′1⟩=¯α∫d2x32πx202x203x232(⟨S03S32S2′1⟩−⟨S02S2′1⟩)+¯α∫d2x32πx212′x213x232′(⟨S2′3S31S02⟩−⟨S02S2′1⟩), (12)

where we have introduced the notation . (See Fig. 6a for the corresponding graphical representation.)

This equation calls for a new equation for the 3-point correlators, and so on. The obtained hierarchy is nothing but the Balitsky hierarchy [18] (see also Ref. [46, 47]) restricted to dipoles. We refer the reader to [48, 49] (see also Ref. [50]) for the detailed relationship of this equation to the B-JIMWLK formalism. Note in particular that the factorized correlators (10) is a solution of the whole hierarchy, and is actually a good approximation to the solution of the full B-JIMWLK equations. This statement was first made after the results of the numerical solution to the JIMWLK equation worked out in Ref. [51].

We may wonder why there are no terms involving one-point functions in the right handside of the previous equation. Actually, such terms would correspond to graphs like the one of Fig. 6b, in which, for example, two dipoles merge. They are expected to occur if saturation is properly taken into account. While the restriction of the Balitsky equation to dipoles does not drastically change the solution for the scattering amplitudes, such terms would instead have a large effect, as we shall discover in the following. In the next section, we will explain why such terms are actually required for physical reasons.

#### 2.1.4 Saturation

The BK equation may be well-suited for the ideal case in which the target is a nucleus made of an infinity of independent nucleons. But it is not quite relevant to describe the scattering of more elementary objects such as two dipoles (or two virtual photons, to be more physical).

Indeed, following Chen and Mueller [52] (see also Ref. [53]), let us consider dipole-dipole scattering in the center-of-mass frame, where the rapidity evolution is equally shared between the projectile and the target (see Fig. 7a). Then at the time of the interaction, the targets are dipoles that stem from the branching of a unique primordial dipole. Obviously, the assumption of statistical independence of the targets, which was needed for the factorization (10) to hold, is no longer justified.

So far, we have seen that nonlinear effects which go beyond the factorization formula (7) are necessary to preserve unitarity as soon as . This came out of an analysis of Eq. (7) in the restframe of the target. The rapidity at which the system reaches this number of dipoles and hence at which the BFKL approach breaks down may be found from the form of the typical growth of with , namely . Parametrically,

 yBFKL∼1¯αlog1α2s. (13)

Now we may go to the center-of-mass frame, where Eq. (7) with would describe the amplitude in the absence of nonlinear effects. There, the typical number of dipoles in the projectile and in the target are well below : . We actually see that the evolution of the dipoles in each of theses systems remains linear until . In that rapidity interval, nonlinear effects consist in the simultaneous scatterings of several dipoles from the target and the projectile but the evolution of still obeys the BFKL equation (see Fig. 7a). Now, performing a boost to the projectile restframe, the evolution goes into the target. Formula (3) should then apply for the amplitude . But if the evolution of the target were kept linear, then the amplitude would not be unitarity because the number of dipoles would be larger than . Hence, through some nonlinear mechanism, which was represented by multiple scatterings between linearly evolving objets in the center-of-mass frame, the dipole number density has to be kept effectively lower than in order to preserve unitarity (see Fig. 7b). This is called parton saturation. The precise saturation mechanism has not been formulated in QCD. It could be dipole recombinations due to gluon fusion, multiple scatterings inside the target which slow down the production of new dipoles [54] (as in Fig. 7b), “dipole swing” as was proposed more recently [55, 56], or any other mechanism.

Hence, unitarity of the scattering amplitudes together with boost-invariance seem to require the saturation of the density of partons.

A pedagogical review of saturation and the discussion of the relationship between saturation and unitarity may be found in Ref. [57]. Original papers include Refs. [58, 59].

#### 2.1.5 The Pomeron language

So far, we have presented in detail a -channel picture of hadronic interactions, and it is in this formalism that we will understand most easily the link with reaction-diffusion processes. In the -channel formulation, all the QCD evolution happens in the form of quantum fluctuations of the interacting hadrons. However, a picture maybe more familiar to the reader is a -channel picture, where the rapidity evolution is put in the -channel, while the projectile and target are in their bare states. This picture directly stems from the usual Lorentz-invariant formulation of quantum field theory, while the dipole picture (or the parton model) is derived in the framework of time-ordered perturbation theory.

Classes of Feynam diagrams can be grouped into “Pomerons” (or Reggeized gluons, see Fig. 8), in terms of which scattering processes may be analyzed. (A pedagogical review on how to derive the BFKL equation in such a formalism is available from Ref. [60]).

An effective action containing Pomeron fields and vertices may be constructed. In these terms, the -channel diagrams of Figs. 5,7a may be translated in terms of the diagrams of Fig. 9. The effective action formalism was initially developped in Refs. [61, 62, 63]. More recently, there has been some progress in the definition of the effective action [64], some of it with the help of the correspondence with statistical physics processes [65, 66].

We will not expand on this formulation in the present review, because it is difficult to see the analogy with statistical processes in this framework. A -channel picture is much more natural. However, a full solution of high-energy QCD may require to go back to that kind of calculation and compute accurately the Pomeron vertices. This program was formulated some time ago [67, 68], and there is continuing progress in this direction (see e.g. [69, 70]).

### 2.2 Analogy with reaction-diffusion processes

We are now in position to draw the relationship between high-energy QCD and reaction-diffusion processes. In the first section below, we will show that the BK equation is, in some limit, an equation that also appears in the context of statistical physics. Second, we will exhibit a particular reaction-diffusion model, and show in the final section how this model is related in a more general way to scattering in QCD.

#### 2.2.1 The BK equation and the FKPP equation

Let us first show at the technical level that under some well-controlled approximations, the BK equation (11) may be mapped exactly to a parabolic nonlinear partial differential equation. This observation was first made in Ref. [38].

To simplify, we will look for impact-parameter independent solutions: is supposed to depend on and only, not on . We switch to momentum space through the Fourier transformation

 A(y,k)=∫d2x012πx201eikx01A(y,x01). (14)

This transformation greatly simplifies the BK equation [24, 25]. It now reads

 ∂¯αyA(y,k)=χ(−∂logk2)A(y,k)−A2(y,k). (15)

The first term in the right handside, which is a linear term, is actually an integral kernel, obtained by Fourier transformation of the BFKL kernel (first three terms in the right handside of Eq. (11)). It is most easily expressed in Mellin space: is the set of its eigenfunctions, with the corresponding eigenvalues

 χ(γ)=2ψ(1)−ψ(γ)−ψ(1−γ). (16)

This kernel may be expanded around some real , fixed between 0 and 1. Keeping the terms up to is the well-known diffusive approximation, which is a good approximation for large rapidities. Introducing the notations , and , the BK equation reads, within this approximation

 ∂¯αyA=χ′′02∂2logk2A+(γ0χ′′0−χ′0)∂logk2A+(χ0−γ0χ′0+γ20χ′′02)A−A2. (17)

Through some linear change of variable ,

 ¯αy=tχ0−γ0χ′0+γ20χ′′02logk2= ⎷χ′′02(χ0−γ0χ′0)+γ20χ′′0x+γ0χ′′0−χ′0χ0−γ0χ′0+γ20χ′′02t, (18)

one may get rid of the first-order partial derivative in the right handside. We then find that the new function

 u(t,x)=A(y(t),logk2(t,x))χ0−γ0χ′0+γ20χ′′02 (19)

obeys the equation

 ∂u(t,x)∂t=∂2u(t,x)∂x2+u(t,x)−u2(t,x), (20)

which is the Fisher [36] and Kolmogorov-Petrovsky-Piscounov [37] (FKPP) equation. This equation was first written down as a model for gene propagation in a population in the large population size limit. But it turns out to apply directly or indirectly to many different physical situations, such as reaction-diffusion processes, but also directed percolation, and even mean-field spin glasses [71]. A recent comprehensive review on the known mathematics and the phenomenological implications of the FKPP equation can be found in Ref. [72].

As a side remark, we note that if is chosen such as , then the mapping drastically simplifies. Actually, this choice has a physical meaning, as we will discover in Sec. 4 when we try and solve the BK equation.

Beyond the exact mapping (20) between an approximate form of the BK and the FKPP equations, the full BK equation is said to be in the universality class of the FKPP equation. All equations in this universality class share some common properties, as will be understood below. The exact form of the equation is unessential. As a matter of fact, recently, it has been checked explicitely that the BFKL equation with next-to-leading order contributions to the linear evolution kernel (but keeping the QCD coupling fixed) is also in the same universality class. A mapping to a partial differential equation (which involves higher-order derivatives in the rapidity variable) was exhibited [73]. What defines physically the universality class of the FKPP equation is a branching diffusion process with some saturation mechanism. The details seem unimportant.

In the next section, we shall give a concrete example of a reaction-diffusion process: We will see how the FKPP equation appears as a fluctuationless (or “mean-field”) limit of some stochastic reaction-diffusion process. In Ref. [38], it had not been realized that the analogy of QCD with such processes is in fact much deeper than the formal mapping between the BK equation and the FKPP equation that we have just outlined. But this is actually the case, as we shall shortly argue.

#### 2.2.2 Reaction-diffusion processes: An example

We consider a reaction-diffusion model, which was introduced in Ref. [74]. Particles are evolving in discrete time on a one-dimensional lattice. At each timestep, a particle may jump to the nearest position on the left or on the right with respective probabilities and , and may split into two particles with probability . We also allow that each of the particles on site at time to die with probability .

From these rules, we may guess what a realization of this evolution may look like at large times. The particles branch and diffuse (they undergo a linear evolution) until their number becomes of the order of , at which point the probability that they “die” starts to be sizable, in such a way that their number never exceeds by a large amount, on any site. But if the initial condition is spread on a finite number of lattice sites, the linear branching-diffusion process may always proceed towards larger values of , where there were no particles in the beginning of the evolution. Hence a realization will look like a front connecting an ensemble of lattice sites where a quasi-stationary state in which the number of particles is (up to fluctuations) has been reached, to an ensemble of empty sites (towards ). This front moves with time as the branching diffusion process proceeds. The position of the front may be defined in different ways, leading asymptotically to equivalent determinations, up to a constant. For example, one may define as the rightmost bin in which there are more than particles, or, alternatively, as the total number of particles in the realization whose positions are greater than 0, scaled by . A realization and its time evolution is sketched in Fig. 10.

Between times and , particles out of move to the left and of them move to the right. Furthermore, particles are replaced by their two offspring at , and particles disappear. Hence the total variation in the number of particles on site reads

 n(t+Δt,x)−n(t,x)=−nl(t,x)−nr(t,x)−n−(t,x)+n+(t,x)+nl(t,x+Δx)+nr(t,x−Δx). (21a) The numbers describing a timestep at position x have a multinomial distribution: P({nl,nr,n+,n−})=n!nl!nr!n+!n−!Δn!pnllpnrr×λn+(λn/N)n−(1−pl−pr−λ−λn/N)Δn, (21b)

where , and all quantities in the previous equation are understood at site and time . The evolution of is obviously stochastic. One could write the following equation:

 u(t+Δt,x)=⟨u(t+Δt,x)⟩+√⟨u2(t+Δt,x)⟩−⟨u(t+Δt,x)⟩2ν(t+Δt,x) (22)

where the averages are understood over the time step that takes the system from to . They are conditioned to the value of at time . is a noise, i.e. a random function. The equation was written in such a way that it has zero mean and unit variance. Note that the noise is updated at time in this equation.

One can compute the mean evolution of in one step of time which appears in the right handside of Eq. (22) from Eq. (21). It reads

 ⟨u(t+Δt,x)|{u(t,x)}⟩=u(t,x)+pl[u(t,x+Δx)−u(t,x)]+pr[u(t,x−Δx)−u(t,x)]+λu(t,x)[1−u(t,x)]. (23)

The mean evolution of the variance of that appears in Eq. (22) may also be computed. The precise form of the result is more complicated, but roughly speaking, the variance of after evolution over a unit of time is of the order of for small . This is related to the fact that the noise has a statistical origin: Having particles on the average in a system means that each realization typically consists in particles.

When is infinitely large, one can replace the ’s in Eq. (23) by their averages: This would be a mean field approximation. Obviously, the noise term drops out, and the equation becomes deterministic. Note that if we appropriately take the limits and , setting

 λ=Δt,  pR=pL=Δt(Δx)2, (24)

the obtained mean-field equation is nothing but the FKPP equation (20). For the numerical simulations of this model that we will perform in Sec. 4, we will keep and finite, which is usually more convenient for computer implementation.

Thus we have seen that the evolution of reaction-diffusion systems is governed by a stochastic equation (22) whose continuous limit () and mean-field limit () is a partial differential equation of the form of (exactly actually, in our simple case study) the FKPP equation. We shall now argue that partons in high-energy QCD form such a system.

#### 2.2.3 Universality class of high-energy QCD

Let us come back to the QCD dipole model. We have seen that rapidity evolution of the hadron wavefunctions proceeds through a branching diffusion process of dipoles. Let us denote by the scattering amplitude of the probe dipole off one particular realization of the target at rapidity and at a given impact parameter. This means that we imagine for a while that we may freeze the target in one particular realization after the rapidity evolution , and probe the latter with projectiles of all possible sizes. Of course, this is not doable in an actual experiment, not even in principle. But it is very important for the statistical picture to go through such a “gedanken observable”. The amplitude , which is related to the measurable total cross section, is nothing but the average of over all possible realizations of the fluctuations of the target, namely

 A(y,r)=⟨T(y,r)⟩. (25)

The branching diffusion of the dipoles essentially occurs in the variable. The scattering amplitude is roughly equal to the number of dipoles in a given bin of (logarithmic) dipole size, multiplied by . From unitarity arguments and consistency with boost-invariance, we have seen that the branching diffusion process should (at least) slow down in a given bin as soon as the number of objects in that very bin is of the order of , in such a way that effectively, the number of dipoles in each bin is limited to . A typical realization of is sketched in Fig. 11. As in the case of the reaction-diffusion process, from similar arguments, it necessarily looks like a front. The position of the front, defined to be the value of for which is equal to some fixed number, say , is related to the saturation scale defined in the Introduction: .

We now see that there is a very close analogy between what we are describing for QCD here and the model that we were introducing in the previous section. So in particular, one might be able to formulate interaction processes in QCD with the help of a stochastic nonlinear evolution equation for the “gedanken” amplitude . We already know the mean-field limit that one should get, when is very large: This is the BK equation, as was rigorously proven above. Thus we know the equivalent of the term in Eq. (22). The noise term is not known, but since it is of statistical origin, we know that it must be of the order of the square root of the number of dipoles normalized to , that is to say, of order . We may write an equation of the form

 ∂¯αyT(y,k)=χ(−∂logk2)T(y,k)−T2(y,k)+αs√2T(y,k)ν(y,k), (26)

where is a noise, uncorrelated in rapidity and transverse momentum, with zero mean and unit variance. (The factor of 2 under the square root is essentially arbitrary). This equation is to be compared to the following one:

 ∂tu(t,x)=∂2xu(t,x)+u(t,x)−u2(t,x)+√2u(t,x)Nν(t,x), (27)

which is the so-called “Reggeon field theory” equation when the noise is exactly a normal Gaussian white noise, that is to say, of zero mean and whose non-vanishing cumulant reads

 ⟨ν(t,x)ν(t′,x′)⟩=δ(t−t′)δ(x−x′). (28)

It is a stochastic extension of Eq. (20). If the noise term were of the form

 √2u(t,x)(1−u(t,x))Nν(t,x) (29)

instead, then this equation would be what is usually referred to as the stochastic Fisher-Kolmogorov-Petrovsky-Piscounov equation. The sFKPP equation and the physics that it represents is reviewed in Ref. [75].

Taking averages over events converts this equation into a hierarchy of coupled equations, which has a lot in common in its structure with the Balitsky hierarchy (12). (Actually, there are some extra terms compared to the Balitsky hierarchy, which were first found from the analogy with reaction-diffusion processes, and which precisely represent nonlinear effects inside the wavefunctions. A detailed study may be found in Ref. [76]). We will perform explicit calculations in this spirit within simpler models in Sec. 3 below.

Based on these considerations, we may establish a dictionary between QCD and reaction-diffusion processes. The correspondence is summarized in Tab. 1.

The mechanism for saturation of the parton densities (i.e. of the dipole number density) is not known for sure in QCD. There are also important differences between the reaction-diffusion model introduced above and QCD that lie in the “counting rule” of the particles (provided by the form of in the QCD case, see Eq. 2). But from the general analysis of processes described by equations in the universality class of the stochastic FKPP equation and the underlying evolution mechanisms presented in Sec. 4, we will understand that most of the observables have universal properties in appropriate limits, which do not depend on the details of the mechanism at work. We draw the reader’s attention to Refs. [77, 78], where a precise stochastic equation was searched for in QCD. Some of the problems one may face with the use and the very interpretation of such equations were studied in Ref. [79].

The way in which we view high energy QCD is actually not particularly original: It is nothing but the QCD dipole model, which was implemented numerically in the form of a Monte Carlo event generator by Salam [80, 81, 82] (see also [55] for another more recent implementation). He also devised and implemented a saturation mechanism [54] that went beyond the original dipole model pictured in Fig. 7a, but which is necessary, as was argued before.

Before discussing more deeply the physical content of equations of the form of Eq. (26), we shall first study a model in which spatial dimensions are left out, that we will be able to formulate in different ways.

## 3 Zero-dimensional model

In the previous section, we have understood that scattering at high energy in QCD may be viewed as a branching-diffusion process supplemented by a saturation mechanism. We have exhibited a simple toy model with these characteristics, whose dynamics is represented by an equation of the type (27).

Unfortunately, even that toy model is too difficult to solve analytically. We shall study a still simplified model, where there is no diffusion mechanism: Realizations are completely specified by the number of particles that the system contains at a given time. Of course, in this case, a saturation scale cannot be defined, which limits the relevance of this model for QCD. However, we will be able to formulate this model in many different ways, and to draw parallels with QCD.

We start by defining precisely the model. Then, two approaches to the computation of the moments of the number of particles are presented. The first set of methods relies on field theory (Sec. 3.2). The second method relies on a statistical approach (Sec. 3.3) and will be extended in a phenomenological way to models with a spatial dimension in Sec. 4. We shall then draw the relation to a scattering-like formulation (Sec. 3.4). Finally (Sec. 3.5), some variants of the basic model are reviewed.

### 3.1 Definition

Let us consider a simple model in which the system is characterized by its number of particles at each time . Between times and , each particle has a probability to split in two particles. For each pair of particles, there is a probability that they merge into one. We may summarize these rules in the following form:

 nt+dt=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩nt+1 % proba ntdtnt−1 proba nt(nt−1)dtNnt proba 1−ntdt−nt(nt−1)dtN. (30)

From this, one can easily derive an equation for the time evolution of the probability of having exactly particles in the system at time :

 ∂P∂t(n,t)=(n−1)P(n−1,t)+n(n+1)NP(n+1,t)−(n+n(n−1)N)P(n,t). (31)

This is the master equation for the Markovian process under consideration. The two first terms with a positive sign represent the process of going from one state containing particles to an adjacent one containing or particles respectively, while the last term simply corrects the probability to keep it unitary.

By multiplying both sides of this equation by and summing over , we get an evolution equation for the average number of particles :

 d⟨nt⟩dt=⟨nt⟩−1N⟨nt(nt−1)⟩ (32)

Unfortunately, this equation is not closed, and one would have to establish an equation for , which would involve 3-point correlators of , and so on, ending up with an infinite hierarchy of equations, exactly like in Sec. 2 for QCD (see Eq. (12)).

This illustrates the difficulties one has to face before one can get an analytical expression for , even in such a simple model.

### 3.2 “Field theory” approach

In the next subsections, we will follow different routes to get analytical results on the moments of the number of particles in the system at a given time . The first one will be similar to the -channel picture of QCD (see Sec. 2), since it will consist in computing the time (rapidity in QCD) evolution of realizations of the system. The second one will be closer to the -channel picture of QCD. We will see how “Pomerons” may appear in these simple systems. We will then examine a formulation in terms of a stochastic nonlinear partial differential equation, which is nothing but the sFKPP equation in which the space variable () has been discarded.

#### 3.2.1 Particle Fock states and their weights

Statistical problems were first formulated as field theories by Doi [83] and Peliti [84]. Different authors have used these methods (see Ref. [85] for a review). We shall start by following the presentation given in Ref. [86].

We would like to interpret the master equation (31) as a quasi-Hamiltonian evolution equation of the type of the ones that appear in quantum mechanics. To this aim, we need to introduce the basis of states of fixed number of particles. We define the ladder operators and by their action on these states:

 a|n⟩=n|n−1⟩, a†|n⟩=|n+1⟩ (33)

and which obey the commutation relation

 [a,a†]=1. (34)

The -particle state may be constructed from the vacuum (zero-particle) state by repeated application of the ladder operator:

 |n⟩=(a†)n|0⟩. (35)

The normalization is not standard with respect to what is usually taken in quantum mechanics. In particular, the orthogonal basis is normalized in such a way that . This implies that the completeness relation reads

 ∑n1n!|n⟩⟨n|=1. (36)

We also introduce the state vector of the system at a time as a sum over all possible Fock states weighted by their probabilities:

 |ϕ(t)⟩=∑nP(n,t)|n⟩. (37)

It is straightforward to see that the master equation (31) is then mapped to the Schrödinger-type equation

 ∂∂t|ϕ(t)⟩=−H|ϕ(t)⟩, (38)

where is the “Hamiltonian” operator

 H=(1−a†)a†a−1N(1−a†)a†a2. (39)

The first term represents the splitting of particles, while the second one, proportional to , represents the recombination. We may rewrite as

 H=H0+H1, (40)

where

 H0=a†a (41)

is the “free” Hamiltonian whose eigenstates are the Fock states. We now go to the interaction picture by introducing the time-dependent Hamiltonian

 HI(t)=eH0tH1e−H0t (42)

and the states . The solution of the evolution reads

 |ϕ⟩I=Texp(−∫t0dt′HI(t′))|ϕ0⟩I=|ϕ0⟩I−∫t0dt′HI(t′)|ϕ0⟩I+∫t0dt′∫t′0dt′′HI(t′)HI(t′′)|ϕ0⟩I+⋯ (43)

We may then compute the weights of the successive Fock states by applying this formula. Let us show how it works in detail by computing the state of a single particle evolved from time 0 to time , in the limit in which there are no recombinations. We follow the usual method to deal with such problems in field theory. We insert repeatedly complete basis of eigenstates of into Eq. (43), namely

 |ϕ⟩I=|1⟩−∫t0dt′∑n11n1!|n1⟩⟨n1|HI(t′)|1⟩+⋯ (44)

(We have kept the first two terms in Eq. (43) explicitely). Using the expression for as a function of and , together with the knowledge that the Fock states are eigenstates of , we get

 |ϕ⟩=e−t|1⟩−∑n1e−n1t∫t0dt′en1t′−t′1n1!|n1⟩⟨n1|H1|1⟩+⋯ (45)

Inserting the expression for , one sees that in the infinite- limit, there is only one possible transition, . Performing the integration over and computing in the same manner the higher orders, one finally gets the expansion

 |ϕ⟩=e−t|1⟩+e−t(1−e−t)|2⟩+⋯+e−t(1−e−t)n−1|n⟩+⋯ (46)

from which one can read the probabilities of the successive Fock states. This expansion is similar to the expansion in dipole Fock states introduced in Sec. 2: The -particle states correspond to -dipole states in QCD, and their weights are computed by applying successive splittings to the system, whose rates are given by Eq. (1). (They are just unity in the case of the zero-dimensional model.)

We see that this method is well-suited to compute the probabilities of the lowest-lying Fock-states, and their successive corrections at finite . But in general we are rather interested in averages such as , for which the weights of all Fock states are needed. We will develop a slightly different (but equivalent) formalism below, that will enable us to get these averages in a much more straightforward way.

#### 3.2.2 Pomeron field theory

Let us introduce the generating function of the factorial moments of the distribution of the number of particles

 Z(z,t)=∑n(1+z)nP(n,t). (47)

The evolution equation obeyed by can easily be derived from the master equation (31):

 ∂Z∂t=z(1+z)(∂Z∂z−1N∂2Z∂z2). (48)

We may represent this equation in a second-quantized formalism by introducing the operators

 b†=z, b=∂∂z=¯z (49)

acting on the set of states consisting in the analytic functions of . Then we may write

 ∂Z∂t=−HPZ, (50)

where

 HP=HP0+HP1, with HP0=−b†b, HP1=−b†b†b+1Nb†(1+b†)b2. (51)

A basis for the states is

 |k⟩=zk, ⟨k|=¯zk (52)

which is orthogonal with respect to the scalar product

 ⟨Z1|Z2⟩=∫dzd¯z2iπe−|z|2¯Z1(z,¯z)Z2(z,¯z), (53)

and obeys the normalization condition . We shall call these states “-Pomeron” states, by analogy with high-energy QCD. We may apply exactly the same formalism as before, since the operators , have the same properties as the , .

From the definition of the scalar product, it is not difficult to see that the -th factorial moment of may be obtained by a mere contraction of the state vector , computed by solving the Hamiltonian evolution, with a -Pomeron state. The following identity holds:

 ⟨k|Z⟩=⟨nt!(nt−k)!⟩, (54)

where the average in the right handside goes over the realizations of the system. As for the initial condition, starting the evolution with one particle means taking as an initial condition the superposition of zero- and one-Pomeron states respectively. The zero-Pomeron state does not contribute to the evolution, hence a one-Pomeron state is like a one-particle state.

In order to simplify the systematic computation of these moments, we may use a diagrammatic method and establish Feynman rules. To this aim, we write the contribution of the graphs with -vertices (corresponding to the term of order in the expansion of Eq. (43)), starting with a one-Pomeron state:

 ⟨k|Z⟩⊃(−1)l∫t0dt1∫t10dt2