text:
# Deconfinement mechanism in three dimensions for gauge fields coupled to bosonic matter fields with fundamental charge

(Received )

## Abstract

We propose a mechanism by which electric charges deconfine in an Abelian Higgs model with matter fields belonging to the fundamental representation of the gauge group. Kosterlitz-Thouless like recursion relations for a scale-dependent stiffness parameter and fugacity are given, showing that for a logarithmic potential between point charges in any dimension, there exists a stable fixed point at zero fugacity, with a dimensionality dependent universal jump in the stiffness parameter at the phase transition.

It is well known that a pure Maxwell theory with gauge fields arising from a compact $`U(1)`$ gauge group permanently confines electric charges in three dimensions Polyakov . In such a case, specifically in the absence of matter fields coupled to the gauge fields, the Wilson loop is a good probe for confinement and satisfies the area law, with a linearly confining potential between electric test charges. The periodic character of the gauge fields induces topological defects in the theory. In the case of compact Maxwell theory in three dimensions these defects are magnetic monopoles. This is easily seen as the $`2\pi `$ jumps can be accounted for by writing the field strength as

$$F_{\mu \nu }=f_{\mu \nu }2\pi ϵ_{\mu \nu \lambda }_\lambda d^3yD(xy)n(y),$$
(1)
where $`f_{\mu \nu }`$ is the non-singular part of the field strength, $`n(x)=_iq_i\delta ^3(xx_i)`$ is the magnetic charge density and

$$^2D(xy)=\delta ^3(xy).$$
(2)
Thus, the dual field strength $`F_\mu ^{}=ϵ_{\mu \nu \lambda }F_{\nu \lambda }/2`$ satisfies

$$_\mu F_\mu ^{}(x)=2\pi n(x).$$
(3)
In the absence of monopoles the potential between electric test charges is given simply by the Coulomb interaction in two space dimensions:

$$V(R)\mathrm{ln}R.$$
(4)
The monopoles of the dual theory, on the other hand, interact through a Coulomb potential in three dimensions, i.e., $`V_{\mathrm{mon}}1/r`$. In order to show that in the compact Maxwell theory electric charges never deconfine, we have to study a classical Coulomb gas of monopoles in three dimensions. Polyakov Polyakov  carried out this study by showing that for low fugacity the Coulomb gas maps onto a sine-Gordon theory, which he solved in a saddle point approximation in the presence of external sources corresponding to the test electric charges. The end result is a potential between electric test charges $`V(R)R`$ in two space dimensions. Thus, the monopoles act in such a way as to produce an anti-screening effect in the potential between electric charges. The force between the electric charges, which decreases as $`1/R`$ in two space dimensions for the non-compact Maxwell theory, does not vary with the distance any longer in the case of compact Maxwell theory.

Note that in two space-dimensions, the absence of monopoles in the gauge-field produces a potential between test charges which is given by $`V(R)\mathrm{ln}R`$. In the presence of fluctuating matter fields in three space-time dimensions (two space dimensions), this is changed to $`V(R)1/R`$, due to the presence of an anomalous scaling dimension of the gauge field. The gauge-field propagator in $`d`$ space-time dimensions in the presence of critical matter-field fluctuations is given by $`D(|x|)1/|x|^{d2+\eta _A}`$, where $`\eta _A`$ is the anomalous scaling dimension of the gauge-field, gauge-invariance dictates that $`\eta _A=4d`$ Herbut , and $`d`$ is the space-time dimensionality. Therefore, we see that $`D(|x|)1/|x|^2`$ when $`d(2,4]`$. This is the same as we would get with or without matter-field fluctuations in four-dimensional space-time. Hence, we see that the role of the anomalous scaling of the gauge field is to essentially produce extra space dimensions. It is as if the test charges were living in three space dimensions, but confined to moving in two.

An alternative derivation of Polyakov’s result amounts to showing that a Coulomb gas of monopoles does not undergo any phase transition between a dielectric and metallic phase in three dimensions. The system is always in a “high-temperature” metallic phase. This can be seen by deriving the corresponding recursion relations for a three-dimensional Coulomb gas. The recursion relations for a $`d`$-dimensional Coulomb gas were derived by Kosterlitz Kosterlitz . The result is

$$\frac{dK^1}{dl}=4\pi ^2y^2(2d)K^1,$$
(5)
$$\frac{dy}{dl}=\left[d2\pi ^2f(d)K\right]y,$$
(6)
where $`f(d)=(d2)\mathrm{\Gamma }[(d2)/2]/(4\pi )^{d/2}`$. Here, $`y(l)`$ and $`K(l)`$ are essentially the scale-dependent fugacity and inverse dielectric constant of the $`d`$-dimensional Coulomb gas, respectively. These renormalization group equations are therefore basically nothing but self-consistency equations for scale-dependent electrostatics. For $`d=2`$, the above equations reduce to the celebrated recursion relations for the Kosterlitz-Thouless (KT) phase transition KT . In this case a line of fixed points occurs in the flow diagram. However, for $`d>2`$ there are no zero-fugacity fixed points to the above recursion relations. For $`d>2`$, the second term on the right hand side in Eq. (5) is positive, which means that even at arbitrarily low fugacity, $`K^1(l)>0`$ will increase indefinitely as the (logarithmic) length scale $`l`$ increases. Hence, no matter how large we make the bare value of $`K(l)`$, it will eventually be reduced enough, by Eq. (5), to make the right hand side of Eq. (6) positive. Hence, it is inevitable that $`y(l)`$ will eventually start increasing with $`l`$, thus destroying the zero-fugacity fixed point well known to exist in the case $`d=2`$ KT . In particular, this result means that no phase transition occurs in $`d=3`$ in the ordinary Coulomb gas. Thus, the Coulomb gas of monopoles is always in the plasma phase. This is nothing but a statement which is equivalent to permanent confinement of electric test charges in compact three-dimensional photodynamics Polyakov .

When matter fields are present, the confinement properties of the theory are likely to be changed. There are, however, many subtleties involved when the matter fields are in the so called fundamental representation of compact $`U(1)`$. To see what are the main points, let us consider the lattice abelian Higgs model, whose action is given by

$$S=\beta \underset{x,\mu }{}\mathrm{cos}[\mathrm{\Delta }_\mu \theta (x)qA_\mu (x)]\kappa \underset{x,\mu ,\nu }{}\mathrm{cos}[F_{\mu \nu }(x)],$$
(7)
where $`q`$ is the charge carried by the Higgs field . The case where the Higgs field carries the fundamental charge ($`q=1`$) differs in an essential way from the case $`q>1`$. This can be seen by considering the limiting cases $`\beta \mathrm{}`$ and $`\kappa \mathrm{}`$. Let us consider first the case $`q=1`$. First of all, for all values of $`q`$ the limit $`\kappa \mathrm{}`$ corresponds to the $`3DXY`$ model. This is so because when $`\kappa \mathrm{}`$ all gauge-field fluctuations are supressed except those that are indistinguishable from vortex-loop fluctuations in the matter sector. Hence, for all $`q`$, the model exhibits a phase transition when $`\kappa \mathrm{}`$. The limit $`\beta \mathrm{}`$, on the other hand, is trivial when $`q=1`$, and there is no phase transition associated with it.

The situation for $`q=2`$ is drastically different since in this case the limit $`\beta \mathrm{}`$ leads to a $`Z_2`$ gauge theory, which exhibits a phase transition in the Ising model universality class when the space-time dimensionality $`d=3`$. Thus, when $`q=2`$ it is natural to think that there is a critical line in the phase diagram of the $`q=2`$ theory that interpolates between the two limiting critical regimes, and this can indeed be demonstrated Bhanot ; Sudbo ; Smiseth . The case corresponding to the Higgs field with the fundamental charge does not have two asymptotic critical regimes to be interpolated. For this reason, it is generally thought that there is no phase transition in the $`q=1`$ three-dimensional lattice Abelian Higgs model and that therefore the theory permanently confines electric test charges FradShe , as is the case in the pure compact Maxwell theory Polyakov .

Another obstacle against a deconfinement phase transition in the $`q=1`$ case comes from Elitzur’s theorem Elitzur . This theorem simply states that averages of non-gauge invariant operators are always zero in the absence of gauge fixing. Only gauge-invariant operators can have a nonzero expectation value. In other words, a local gauge symmetry cannot be spontaneously broken. This is important, since it implies that the Higgs mechanism can only occur upon some gauge fixing. In the lattice action a natural gauge fixing is the unitary gauge parametrization. In the unitary gauge there is a residual global gauge symmetry left, which can in turn be broken. In the $`q=2`$ case the residual global symmetry corresponds to the $`Z_2`$ group. Thus, in this case the Higgs phase can be distinguished from the confinement phase FradShe . In the $`q=1`$ case, however, the residual symmetry is just the identity group and it is therefore trivial. The Higgs phase cannot be distinguished from the confinement phase. Again, there is only one phase and it seems that there is no way out from permanent confinement for fields carrying the fundamental charge.

In the case of continuous global symmetries the Mermin-Wagner theorem MW  forbids spontaneous symmetry breaking in two dimensions. Elitzur’s theorem is far more restrictive than Mermin-Wagner’s theorem, since it applies to any dimension and to discrete gauge groups. A well known way out from the Mermin-Wagner theorem is the KT transition KT , where a phase transition occurs in the absence of long range order. The KT transition occurs precisely in the case of a global $`U(1)`$ group.

Since the global symmetries do not suffer the severe restriction imposed by Elitzur’s theorem, duality transformations where a locally gauge invariant theory is mapped into a globally invariant theory is a powerful tool. We can look for phase transitions there where the symmetry can be spontaneously broken and, in some cases, even to look for phase transitions without spontaneous symmetry breaking, like the KT phase transition at $`d=2`$. Unfortunately, it seems to be generally the case that the phase transition in the dual model corresponds to a nontrivial residual symmetry in the original model. Thus, as far as the phase transition is concerned, in those cases it can easily be established in the original model as well. In such cases, the dual model is still a powerful tool to establish the universality class. It would be interesting if a KT-like phase transition could occur in the dual three-dimensional theory, elevating the three-dimensional $`q=1`$ Abelian Higgs model to a system sustaining a confinement-deconfinement transition. Recently it was pointed out in Ref. KNS  that such a transition indeed occurs in this case. A deconfinement phase transition is driven by a KT-like phase transition in the monopole plasma. Note that in general, the topological defects of theory are more complicated objects. The monopoles are generally connected by magnetic vortex lines, and there are also closed vortex loops. However, at the critical point the vortex line looses tension and we have once more an effective description in terms of a gas of magnetic monopoles. The difference to the usual Coulomb gas is that the interaction between the monopoles are no longer $`1/r`$ as in the pure compact Maxwell theory. The matter fields induce an anomalous scaling behavior and the monopole-monopole interaction becomes $`\mathrm{ln}r`$ in three dimensions KNS ; KNS1 . This behavior strongly suggests that a KT-like transition may also occur for this anomalous Coulomb gas in three dimensions.

For low fugacity the anomalous Coulomb gas can be brought in the form of a sine-Gordon theory with a $`1/|p|^3`$ free propagator KNS ; KNS1 :

$$S=\frac{1}{8\pi ^2K}d^3x[\phi (^2)^{3/2}\phi 2z\mathrm{cos}\phi ],$$
(8)
where $`K=1/g^2`$, with $`g`$ being the gauge coupling. From the above sine-Gordon theory we could in principle raise the following objection to a KT-like behavior at $`d=3`$. It could happen that a $`\phi (^2)\phi `$ term is generated by fluctuation effects. If the theory behaves in such a way, the generated term would clearly dominate at large distances and the $`\phi (^2)^{3/2}\phi `$ would become irrelevant in the renormalization group (RG) sense. The resulting effective action would just correspond to an ordinary sine-Gordon theory in three dimensions and, therefore, no phase transition occurs in this case. The $`\mathrm{ln}r`$ interaction is screened into a $`1/r`$ potential. This argument, if correct, would spoil the deconfinement phase transition for the $`q=1`$ case. Let us show that this is not the case by using two different arguments.

The first argument relies on simple power counting. It turns out that $`d=3`$ is the lower critical dimension of the problem since the field $`\phi `$ is dimensionless. This means that renormalization proceeds in precisely the same way as in the $`d=2`$ well known counterpart of the present problem.

The second argument relies on exact scaling and duality properties of the theory. The sine-Gordon action (8) is a result of a duality transformation of the critical effective action:

$$S_{\mathrm{eff}}d^3xF_{\mu \nu }\frac{1}{\sqrt{^2}}F_{\mu \nu }.$$
(9)
The above corresponds to an exact scaling behavior of the gauge field propagator in the three-dimensional Abelian Higgs model Herbut . Next, we proceed by reductio ad absurdum to prove that the effective sine-Gordon action (8) gives the dominant contribution at large distances. This is seen as follows. The effective action (9) clearly corresponds to the dominant large distance behavior, i.e., the usual Maxwell term $`F_{\mu \nu }^2`$ is obviously irrelevant in the infrared against the anomalous contribution given in Eq. (9). Here, it is crucial that the sign of coefficient of the anomalous term is positive, to ensure that the usual Maxwell term is indeed irrelevant. A negative sign in front of the anomalous term would make the usual Maxwell term relevant. Next, assume that a $`\phi (^2)\phi `$ term is generated by fluctuation effects in Eq. (8) and that this term has a positive sign. This clearly implies that the $`\phi (^2)^{3/2}\phi `$ term becomes irrelevant in the infrared. Therefore, the corresponding effective sine-Gordon action is of the usual type, thus dualizing back to an ordinary Maxwell theory. This contradicts the exact result that the dominant effective critical theory is given by Eq. (9).

Next, to further substantiate our scenario, we derive the corresponding recursion relations for a gas of point-charges in three dimensions interacting through a logarithmic pair-potential, i.e., a three-dimensional logarithmic plasma or anomalous Coulomb gas, and for which Eq. (8) is a field-theoretical description KNS . As in the case of Eqs. (5) and (6), we will find it useful to consider the problem in $`d`$ dimensions. Furthermore, we will consider a more general propagator of the form $`1/|p|^\sigma `$ for the anomalous sine-Gordon model. Thus, we consider a bare potential given by $`U_0(r)=4\pi ^2KV(r)`$, where

$$V(r)=\frac{\mathrm{\Gamma }\left(\frac{d\sigma }{2}\right)}{2^\sigma \pi ^{d/2}\mathrm{\Gamma }(\sigma /2)}[(\mathrm{\Lambda }r)^{\sigma d}1],$$
(10)
with $`\mathrm{\Lambda }`$ being an ultraviolet cutoff. For the particular case $`d=\sigma `$, corresponding to the lower critical dimension in this generalized problem, we have

$$V(r)|_{d=\sigma }=\frac{2^{1\sigma }\pi ^{\sigma /2}}{\mathrm{\Gamma }(\sigma /2)}\mathrm{ln}(\mathrm{\Lambda }r).$$
(11)
The bare electric field is given by $`E_0(r)=4\pi ^2KA(d,\sigma )r^{\sigma d1}/r_0^{\sigma d}`$, where $`r_01/\mathrm{\Lambda }`$ and $`A(d,\sigma )=(d\sigma )\mathrm{\Gamma }[(d\sigma )/2]/[2^\sigma \pi ^{d/2}\mathrm{\Gamma }(\sigma /2)]`$. The bare electric field is renormalized by the other dipoles which are treated as a dielectric medium. The renormalized electric field is then given by

$$E(r)=\frac{4\pi ^2KA(d,\sigma )r^{\sigma d1}}{\epsilon (r)},$$
(12)
where $`\epsilon (r)`$ is the scale-dependent dielectric constant of the medium. We can write this in the form $`\epsilon (r)=1+S_d\chi (r)`$, where $`S_d=2\pi ^{d/2}/\mathrm{\Gamma }(d/2)`$, and the electric susceptibility

$$\chi (r)=N_d_{r_0}^rn(s,\theta )\alpha (s)s^{d1}\mathrm{sin}^{d2}\theta dsd\theta ,$$
(13)
with $`N_d=2\pi ^{(d1)/2}/\mathrm{\Gamma }[(d1)/2]`$. In Eq. (13), $`n(r,\theta )`$ is the density of pairs and $`\alpha (r)`$ is the polarizability of a dipole. We have for small separation of the pairs

$$\alpha (r)=\frac{4\pi ^2Kr^2}{d}+𝒪(r^4).$$
(14)
For $`n(r,\theta )`$ we have at lowest order in the bare fugacity $`y_0`$

$$n(r)=\frac{y_0^2}{r_0^{2d}}e^{U(r)},$$
(15)
where $`U(r)`$ is the effective potential obtained by integrating the renormalized electric field:

$$U(r)=U(r_0)+4\pi ^2KA(d,\sigma )_{r_0}^r𝑑s\frac{s^{\sigma d1}}{\epsilon (s)}.$$
(16)
The effective stiffness $`K_{\mathrm{eff}}(l)`$ is related to the dielectric constant $`\epsilon (r)`$ through

$$\frac{1}{K_{\mathrm{eff}}(l)}=\frac{\epsilon (r_0\mathrm{exp}l)}{K}e^{(\sigma d)l}.$$
(17)
Let us define $`u(l)=U(r_0\mathrm{exp}l)`$, out of which we obtain

$$u(l)=u(0)+4\pi ^2A(d,\sigma )_0^l𝑑vK_{\mathrm{eff}}(v).$$
(18)
Thus,

$$\frac{du}{dl}=4\pi ^2A(d,\sigma )K_{\mathrm{eff}}(l).$$
(19)
We next define the square of the effective fugacity as follows

$$y^2(l)=\frac{2S_d^2}{dr_0^{\sigma 2}}y_0^2e^{(2d\sigma +2)lu(l)}.$$
(20)
From Eqs. (17), (19), and (20) we finally obtain

$$\frac{dK_{\mathrm{eff}}^1}{dl}=4\pi ^2y^2(\sigma d)K_{\mathrm{eff}}^1,$$
(21)
$$\frac{dy}{dl}=\left[d\eta _y2\pi ^2A(d,\sigma )K_{\mathrm{eff}}\right]y,$$
(22)
where the anomalous dimension of the fugacity is given by $`\eta _y=(\sigma 2)/2`$. When $`\sigma =2`$, we recover the recursion relations (5) and (6), which were originally derived by a completely different method Kosterlitz . Such a theory does not exhibit a phase-transition for $`d>2`$. The case relevant for the three-dimensional Abelian Higgs model is, on the other hand, $`\sigma =d=3`$. In this case Eqs. (21) and (22) are very similar to the usual KT recursion relations, except for the presence of the anomalous scaling dimension of the fugacity, $`\eta _y`$ which is nonzero in our case and given by $`\eta _y=1/2`$.

By integration of the above recursion relations Eqs. (21), (22), for the case $`d=\sigma =3`$, we may compute explicitly the screened effective potential $`u(l)`$ on the ordered side, the result is Note

$`u(l)u(0)={\displaystyle \frac{1}{\omega _{}\omega _+}}\left[{\displaystyle \frac{5}{2}}\omega _+l+\mathrm{ln}\left({\displaystyle \frac{\omega _+e^{2\theta }+\omega _{}}{\omega _+e^{2u}+\omega _{}}}\right)\right]`$ (23)
with $`\omega _\pm =1\pm \omega `$, $`u=(5/2)\omega l+\theta `$, and $`\omega `$ and $`\theta `$ are integration constants determined from initial conditions on the flow equations. Asymptotically, we have $`u(l)l\mathrm{ln}(r/r_0)`$, which shows that the effective potential is also logarithmic after screening effects are taken into account. Hence, we conclude that the statement, alluded to above, that the transition is destroyed by screening a bare $`\mathrm{ln}(r)`$-potential into a $`1/r`$-potential, is not correct. This can also be seen from a simple Debye-Hückel theory for a $`\mathrm{ln}(r)`$-potential in $`d`$ dimensions. Indeed, due to Gauss’s theorem in $`d`$-dimensions, the field equation in the corresponding Debye-Hückel theory is given by Dixit  $`(𝐄r^{2d})=S_d[q\delta ^d(𝐫)+\rho (𝐫)]`$, where $`S_d=2\pi ^{d/2}/\mathrm{\Gamma }(d/2)`$ and $`\rho (𝐫)`$ is the variation of the charge density in a plasma with density $`n_0`$ Ichimaru . The electric field is $`𝐄=U`$, where $`U`$ is the screened effective potential. The high temperature limit corresponds to the Debye-Hückel approximation, where the differential equation for $`U(r)`$ can be solved exactly Dixit . The result is

$$U(r)=\frac{2q}{d}K_0[(r/\lambda _D)^{d/2}]q\mathrm{ln}\lambda _D,$$
(24)
where $`K_0`$ is a modified Bessel function of the second kind and the inverse of the Debye screening length $`\lambda _D`$ is given by

$`\lambda _D^1=\left({\displaystyle \frac{16q^2n_0\pi ^{d/2}}{d^2T\mathrm{\Gamma }(d/2)}}\right)^{1/d},`$ (25)
where $`T`$ is the temperature. For $`r/\lambda _D1`$, Eq. (24) has the expansion $`U(r)=2q(\mathrm{ln}2\gamma )/dq\mathrm{ln}r+𝒪(r^2/\lambda _D^2)`$.

The recursion relations Eqs. (21) and (22) thus predict the possibility of a topological phase-transition in the system of point charges interacting with the potential Eq. (10) from a “low-temperature” dielectric phase to a “high-temperature” metallic phase, with a universal jump in the stiffness of the theory given by $`K_{\mathrm{eff}}^{}=2/5`$. This should be contrasted with the universal jump $`K_{\mathrm{eff}}^{}=2/\pi `$ that is found in the two-dimensional case. A specific realization of such a KT-like scenario has recently been suggested in the context of discussing the physics of strongly correlated fermion systems at zero temperature in two spatial dimensions KNS .

In summary, we have proposed a mechanism in three space-time dimensions by which a deconfinement transition occurs for $`U(1)`$ compact gauge fields coupled to bosonic fundamental matter fields. The proposed mechanism relies on a KT-like phase transition in three space-time dimensions. Therefore, no symmetry breaking is involved and Elitzur’s theorem is not violated. Instead, a topological phase transition occurs.


format: text

text:
# Geometric scaling in exclusive processes

## I Introduction

The striking phenomenon of geometric scaling discovered in the HERA data Stasto:2000er  has triggered a wide interest. The inclusive deep-inelastic scattering (DIS) cross section $`\sigma (x,Q)`$ measured as a function of the photon virtuality $`Q`$ and the Bjorken variable $`x`$, was shown to depend on the single combined variable $`Q^2(x/x_0)^{\lambda _0}`$ where $`\lambda _00.3`$ and $`x_0310^4`$. Empirically, this scaling turns out to be valid for $`x10^2`$ and for all available values of $`Q^2`$ (in practice, for $`Q^2=0`$ to $`400\text{GeV}^2`$). Geometric scaling has also been found more recently in deep-inelastic scattering off nuclei Freund:2002ux .

A similar scaling behaviour was discussed a long time ago in the framework of elastic proton-proton scattering Buras:1974km ; Amaldi:1980kd . At that time, an impact parameter analysis revealed that the rise of the cross section as a power of the center-of-mass energy was due to an effective growth of the size of the interaction region, which was almost black in its center. Geometric scaling in $`\gamma ^{}p`$ scattering was discovered in the context of the discussion of saturation effects, and is often considered as an evidence for the fact that the proton looks already quite black, or “saturated”, at HERA.

The scaling variable turns out to be the dimensionless ratio of the photon virtuality $`Q`$ and the so-called saturation scale, which has the following dependence upon the Bjorken variable:

$$Q_s(x)1\text{GeV}\left(\frac{x}{x_0}\right)^{\lambda _0/2}.$$
(1)
The latter is believed to characterize the mean transverse momentum of partons with longitudinal momentum fraction $`x`$. The growth of this typical momentum scale with energy $`W=Q/\sqrt{x}`$ is understood within several different formulations of dense parton evolution Balitsky:1996ub ; Kovchegov:1999yj ; McLerran:1994ni ; McLerran:1994ka ; McLerran:1994vd ; Iancu:2001md ; Iancu:2000hn ; Ferreiro:2001qy  which technically lead to non-linear integro-differential equations. Geometric scaling was seen in numerical simulations of these models Lublinsky:2001bc ; Golec-Biernat:2001if . It was shown to be consistent with QCD evolution in the range in which it is verified in the data Kwiecinski:2002ep . In Ref.Iancu:2002tr ; Mueller:2002zm ; Triantafyllopoulos:2002nz , it was derived from perturbative QCD for $`Q`$ above the saturation scale as a feature of the Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation Kuraev:1977fs ; Balitsky:1978ic  with appropriate boundary conditions. In Ref.Munier:2002gf , the symmetry of the quantity $`Q/Q_s\sigma (Q/Q_s)`$ under the interchange of $`Q`$ and $`Q_s`$ already noticed in Ref.Stasto:2000er  was interpreted as a possible evidence for a proton looking like a collection of independent dipoles of typical size $`1/Q_s`$. The generalized vector dominance model combined with the color dipole model as formulated in Kuroda:2002ry  also leads to geometric scaling. More recently, geometric scaling for all $`Q^2`$ was deduced from the critical behaviour of the correlation function of Wilson lines which appears at small $`x`$ in a near-light-cone Hamiltonian formalism Pirner:2002us . This approach looks fundamentally different from the other ones, in the sense that it does not rely on the parton model.

Most of the latter discussions were only concerned with inclusive quantities. The proton was assumed an homogeneous disk, and its size large with respect to the size of the probe, so that boundary effects could be neglected. The reason for these assumptions is the technical complication of the saturation equations when full dependence upon transverse coordinates is taken into account.

However, one knows that the saturation scale itself depends on the impact parameter. Indeed, it is related to the density of partons inside the proton, which obviously decreases smoothly when one moves from the center to the periphery. The saturation scale in the ideal case described before was only a kind of effective one. Actually, it is theoretically not clear what this scale is in a real proton, because all momenta from $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ to $`Q`$ are a priori expected to contribute to cross sections.

The importance of impact parameter studies for the phenomenology of saturation effects at HERA was explained in Ref.Munier:2001nr . A theoretical discussion of the impact parameter dependence of the saturation scale $`Q_s(x,b)`$ was provided in Ferreiro:2002kv . In the latter paper, some arguments were put forward in favor of a “local” geometric scaling, which would manifest itself in the fact that cross sections for exclusive processes only depend on the ratio of scales $`Q/Q_s(x,b)`$. Let us also stress that geometric scaling for exclusive processes has recently been predicted qualitatively in the context of the discussion of generalized parton distributions Freund:2002ff .

The aim of this paper is to show that “local” geometric scaling is testable experimentally. In the first section, we explain why this new scaling is expected. We then propose a method to investigate it at HERA. Finally, we discuss the qualitative results obtained so far by applying our method to the data.

## II Saturation and geometric scaling

In this section, we show how local scaling appears within different theoretical approaches. We consider the idealized process of scattering of a dipole of size $`1/Q`$ off a dipole of size $`1/\mathrm{\Lambda }`$. $`\mathrm{\Lambda }`$ and $`Q`$ are fixed momentum scales, which are taken in the perturbative region in the following discussion, and the ordering $`Q>\mathrm{\Lambda }`$ is assumed. The small dipole is called “projectile” and the large one “target”. In the end, $`\mathrm{\Lambda }`$ will be formally extrapolated to $`\mathrm{\Lambda }_{\mathrm{QCD}}`$. We denote by $`y\mathrm{log}(1/x)`$ the maximum rapidity available at c.o.m. energy $`W=Q/\sqrt{x}`$.

### II.1 Global and local scaling above the saturation scale from linear evolution

We explore here the perturbative picture of high energy scattering. We consider dipole-dipole scattering at fixed impact parameter $`b`$ as given by the BFKL equation. We assume that $`b`$ is larger than the sizes of the initial-state dipoles: $`b1/Q,1/\mathrm{\Lambda }`$. The Mellin representation of the solution of the BFKL equation then reads Navelet:1998tx ; Munier:1998nj

$$𝒜_\mathrm{d}(y,Q,b)=\frac{1}{2}\frac{1}{Q^2\mathrm{\Lambda }^2b^4}\frac{d\gamma }{2i\pi }(12\gamma )d(\gamma )\left(16b^2Q\mathrm{\Lambda }\right)^{2\gamma }e^{\overline{\alpha }_sy\chi (\gamma )},$$
(2)
where the integration goes over a line in the complex plane parallel to the imaginary axis and intersecting the real axis between 0 and 1. In our conventions, the $`S`$ and $`T`$ matrices are related by $`S=1T`$, so that generally speaking the amplitudes are essentially real at high energy. This is true in particular for the BFKL solution (2). The following standard notations have been used: $`\overline{\alpha }_s=\alpha _sN_c/\pi `$ and $`\chi (\gamma )=2\psi (1)\psi (\gamma )\psi (1\gamma )`$. The function $`d(\gamma )=\alpha _s^2/(16\gamma ^2(1\gamma )^2)`$ is the Mellin transform of the elementary dipole-dipole elastic amplitude with respect to the ratio of their sizes, with appropriate normalization.

The forward amplitude $`𝒯_\mathrm{d}`$ is obtained by integrating (2) over the impact parameter $`b`$. As formula (2) holds only in a large $`b`$ approximation, we must set the lower bound of the integral to a number of the order of the size of the largest interacting dipole. We choose $`1/(4\mathrm{\Lambda })`$ and check the consistency a posteriori. The forward amplitude then reads

$$𝒯_\mathrm{d}(y,Q)=2\pi _{1/(4\mathrm{\Lambda })}^{\mathrm{}}𝑑bb𝒜_\mathrm{d}(y,Q,b)=\frac{8\pi }{Q^2}\frac{d\gamma }{2i\pi }d(\gamma )\left(\frac{Q}{\mathrm{\Lambda }}\right)^{2\gamma }e^{\overline{\alpha }_sy\chi (\gamma )}.$$
(3)
One checks the normalizations by computing the total cross section, using the optical theorem which reads $`\sigma =2𝒯_\mathrm{d}`$ in our conventions. Setting $`y=0`$, one recovers the well-known expression of the elementary (lowest order) dipole-dipole cross section Mueller:1994jq .

Let us recall how one obtains geometric scaling from this equation. One is interested in the large $`y`$ region, where a saddle point gives the dominant contribution to the integral in (3). Expanding the function $`\chi (\gamma )`$ to second order around $`\gamma =1/2`$, the amplitude (3) reads

$$𝒯_\mathrm{d}(y,Q)=\frac{4\pi \alpha _s^2}{\sqrt{\pi a^2y}}\frac{1^{}}{\mathrm{\Lambda }^2}\mathrm{exp}\left\{\omega y\frac{1}{2}\mathrm{log}(Q^2/\mathrm{\Lambda }^2)\frac{1}{4a^2y}\mathrm{log}^2(Q^2/\mathrm{\Lambda }^2)\right\},$$
(4)
with the definitions $`\omega 4\overline{\alpha }_s\mathrm{log}2`$ (the BFKL intercept) and $`a^214\zeta (3)\overline{\alpha }_s`$. At fixed $`y`$, the amplitude becomes large (of the order of the area of the target) when the argument of the exponential is positive, which occurs below the scale $`Q=Q_s(y)`$ satisfying the relation

$$Q_s^2(y)=\mathrm{\Lambda }^2e^{\lambda y},\text{with}\lambda a^2\left(\sqrt{1+4\omega /a^2}1\right).$$
(5)
Note that this formula can be compared with the form of the saturation scale previously quoted (1) provided one makes the identifications $`\lambda \lambda _0`$ and $`\mathrm{\Lambda }1\text{GeV}x_0^{\lambda _0/2}`$. The equation $`Q^2=Q_s^2(y)`$ defines the so-called critical line below which, per definition, the amplitude (3) is large and non-linear corrections prove to be necessary to ensure probability conservation. The interesting feature of Eq.(4) is that it can be rewritten as a function of the ratio of the variables $`Q`$ and $`Q_s`$, up to some approximation. Putting (4) and (5) together, one arrives at

$$𝒯_\mathrm{d}(y,Q)=\frac{4\pi \alpha _s^2}{\sqrt{\pi a^2y}}\frac{1}{\mathrm{\Lambda }^2}\mathrm{exp}\left\{\frac{1}{a^2y}\mathrm{log}^2\frac{Q}{Q_s(y)}\right\}\left(\frac{Q}{Q_s(y)}\right)^{\sqrt{1+4\omega /a^2}}.$$
(6)
The factor which still contains the rapidity alone is slowly varying, hence it can be approximated by its value for a fixed typical rapidity $`y_0`$. Moreover, if $`Q`$ is sufficiently close to $`Q_s`$, one can trade the exponential for $`1`$. This approximation applies when $`|\mathrm{log}Q/Q_s|a^2y\sqrt{1+4\omega /a^2}`$. Introducing the notation $`\gamma (\sqrt{1+4\omega /a^2})/2`$, the amplitude rewrites

$$𝒯_\mathrm{d}(Q,Q_s(y))=\frac{4\pi \alpha _s^2}{\sqrt{\pi a^2y_0}}\frac{1}{\mathrm{\Lambda }^2}\left(\frac{Q_s^2(y)}{Q^2}\right)^\gamma ,$$
(7)
which exhibits geometric scaling.

To make contact with HERA phenomenology, one has to set $`\mathrm{\Lambda }=\mathrm{\Lambda }_{\mathrm{QCD}}`$ and to average the projectile size $`r`$ with weight given by the squared photon wave function. This is the well-known dipole model Nikolaev:1991ja ; Mueller:1994rr  which provides an appealing and useful picture of DIS in the leading-logarithmic approximation of QCD. The virtual photon-proton cross section reads

$$\sigma _{\gamma ^{}p}(y,Q)=d^2r\psi _\gamma ^{}^{}(r,Q)\psi _\gamma ^{}(r,Q)2𝒯_\mathrm{d}(y,1/r).$$
(8)
$`\psi _\gamma ^{}`$ is the photon wave function on a $`q\overline{q}`$ colour dipole. Besides the transverse size $`r`$ and virtuality of the photon $`Q`$, this wave function also depends on the longitudinal momentum fraction $`z`$ of the antiquark. The notation “$``$” stands for an integration over $`z`$. For negligible quark masses, $`\psi _\gamma ^{}^{}(r,Q)\psi _\gamma ^{}(r,Q)`$ is a function of $`rQ`$ only which vanishes rapidly at large $`rQ`$ Stasto:2000er . Thus the scaling property of the dipole-proton amplitude (7) is transmitted to the photon-proton cross section itself which then exhibits the same behaviour as the amplitude $`𝒯_\mathrm{d}`$ in Eq.(7).

The formula just established would give a reasonable description of the data in a limited kinematical range, provided one considers $`\omega `$ and $`a^2`$ as effective parameters. Indeed, the values given by the leading order BFKL equation are known to be incompatible with the data on one hand, and to receive large corrections from next-to-leading order terms on the other hand. The parameters have to be adjusted such that $`\lambda \lambda _00.3`$ and $`\gamma 1`$ Golec-Biernat:1999qd ; Munier:2002gf . Then formulae (1), (5) show that the choice $`x_0=310^4`$ implies $`\mathrm{\Lambda }_{\mathrm{QCD}}=300\text{MeV}`$, which is reasonable. To obtain more reliable expressions, neglected effects such as diffusion into the infrared have to be taken into account. As shown in Ref.Mueller:2002zm , the latter have sizable effects both on the parametric form of the saturation scale (5) and on the expression of the amplitude (7), but do not spoil the overall picture.

The previous analysis can be taken over to the impact parameter dependent amplitude (2) in a straightforward way. The steepest-descent method gives (see Kovner:2001bh )

$$𝒜_\mathrm{d}(y,Q,b)=\frac{128\pi \alpha _s^2}{(\pi a^2y)^{3/2}}\mathrm{log}(16b^2Q\mathrm{\Lambda })\mathrm{exp}\left\{\omega y\mathrm{log}(16b^2Q\mathrm{\Lambda })\frac{1}{a^2y}\mathrm{log}^2(16b^2Q\mathrm{\Lambda })\right\}.$$
(9)
This amplitude is again of order $`1`$ for $`Q=Q_s`$, the saturation scale $`Q_s`$ being given by

$$Q_s^2(y,b)=\frac{1}{256\mathrm{\Lambda }^2b^4}\mathrm{exp}(\lambda y).$$
(10)
Within the approximations already made in the inclusive case, the amplitude can then be expressed as a function of $`Q`$ and $`Q_s`$ only:

$$𝒜_\mathrm{d}(Q,Q_s(y,b))=\frac{64\alpha _s^2}{\sqrt{\pi a^2y_0}}(2\gamma 1)\left(\frac{Q_s^2(y,b)}{Q^2}\right)^\gamma ,$$
(11)
where the saturation scale now depends on the impact parameter $`b`$. Thus we see that the amplitude for each impact parameter is a function of the ratio of the inverse size of the projectile and of the local saturation scale.

However, the kinematical range in which formula (11) may apply is quite restricted. Let us set $`\mathrm{\Lambda }=\mathrm{\Lambda }_{\mathrm{QCD}}`$. We have required large impact parameters (in practice, $`b>0.2\text{fm}`$) and $`Q_s`$ close to $`Q`$, which has to sit in the perturbative regime ($`Q\mathrm{\Lambda }_{\mathrm{QCD}}`$). From formula (10), one sees that the latter conditions require a large rapidity evolution, during which the saturation regime might already be reached for smaller impact parameters (i.e. $`Q<Q_s(b^{})`$ for some $`b^{}<b`$). It is not possible to quantify this effect without knowing the dynamics for all impact parameters.

We shall now refer to another, maybe safer approach, which will however lead to comparable conclusions on “local” geometric scaling.

### II.2 The impact parameter profile as an initial condition

In Ref.Ferreiro:2002kv , the argument that the quantum evolution was quasi-local in impact parameter was developed on the basis of an analysis of the non-linear evolution equations. As a result, the dipole-proton scattering amplitude at fixed impact parameter can be written in a factorized form, because the latter feature means that the impact parameter dependence remains unchanged throughout the rapidity evolution. The non-perturbative profile of the proton is introduced as an initial condition. The validity of such a factorization has been confirmed very recently in the context of $`\gamma ^{}\gamma ^{}`$ reactions Bondarenko:2003ym . The amplitude then reads

$$𝒜_\mathrm{d}(y,Q,b)=S(b)𝒯_\mathrm{d}(y,Q),$$
(12)
where $`𝒯_\mathrm{d}(y,Q)`$ is the forward amplitude (3). Above the saturation scale, $`𝒯_\mathrm{d}(y,Q)`$ exhibits geometric scaling as seen on Eq.(7), thus the amplitude $`𝒜_\mathrm{d}`$ also scales. The profile function

$$S(b)=\frac{2}{\pi }m_\pi ^2e^{2m_\pi b}$$
(13)
was postulated (we chose the normalization such that $`d^2bS(b)=1`$). It involves the long-distance scale $`1/m_\pi `$, which is put by hand and interpreted as an initial condition.

This leads once again to a scaling of the form (11), with a local saturation scale given by

$$Q_s^2(y,b)=e^{2m_\pi b/\lambda }Q_s^2(y)=e^{2m_\pi b/\lambda }\mathrm{\Lambda }_{\mathrm{QCD}}^2e^{\lambda y}.$$
(14)
We see that this result is qualitatively close to Eq.(10), except for the $`b`$-dependence which is exponential rather than powerlike. This stems from the fact that non-perturbative physics with short range pion fields were introduced here, whereas previously only long-range Coulomb-like fields were considered.

At this point, it may be useful to evaluate the saturation scale at HERA from these formulae. For values of the rapidity $`y`$ of the order of those measured at HERA in the small-$`x`$ region ($`xx_0310^4`$), one sees from formula (14) that $`Q_s(b=0)1\text{GeV}`$. $`Q_s`$ reaches the value of $`\mathrm{\Lambda }_{\mathrm{QCD}}`$ at the impact parameter $`b0.6\text{fm}`$. This means that for a probe of size $`1/Q<1\text{GeV}^1`$, the proton looks grey everywhere and extends up to a radius of about $`0.6\text{fm}`$.

## III A way to check scaling at HERA

As geometric scaling at fixed impact parameter seems to be a definite prediction of saturation models, it is worth testing it against the data. A former analysis Munier:2001nr  has already shown the relevance of impact parameter dependent analysis in the discussion of saturation. There, the $`S`$-matrix element for dipole-proton scattering at fixed impact parameter was extracted from diffractive electro-production of vector mesons at HERA. This quantity can be interpreted as a transparency coefficient, and thus quantifies the “blackness” of the proton as seen by the projectile. Such kind of analysis is always model-dependent, since to get the correct normalization of the $`S`$-matrix element, one has to rely on an ad hoc model for the final-state vector meson. Although it was shown in Ref.Munier:2001nr  that requiring that the initial state be a longitudinal photon limits the model-dependence, an estimated uncertainty of 20% was still recognized. However, dependence of the cross section on a scaling variable can be tested in a model-independent way.

The scattering amplitude at fixed impact parameter can be extracted from differential cross sections for high energy quasi-elastic processes by a Fourier transform with respect to the momentum transfer. By quasi-elastic we mean that the initial and final states have the same number of particles. In practice, the process we are thinking of will be diffractive electroproduction of vector mesons.

Let us consider such a quasi-elastic process. Its amplitude $`𝒜`$ is related to the differential cross section through the formula

$$\frac{d\sigma }{dt}=\frac{1}{16\pi }|𝒜(t)|^2,$$
(15)
where $`t`$ is the usual Mandelstam variable and all other dependences are implied. At high energies, it is well-known that the scattering amplitude $`𝒜`$ is essentially real. Formula (15) shows that $`𝒜(t)`$ is proportional to the square root of the differential cross section $`d\sigma /dt`$. A further Fourier transform with respect to the two-dimensional momentum transfer then gives the scattering amplitude for a fixed impact parameter

$$𝒜(b)=\sqrt{\pi }_0^+\mathrm{}d(t)J_0(b\sqrt{|t|})\sqrt{\frac{d\sigma }{dt}},$$
(16)
where $`J_0`$ is a Bessel function. For formula (16) to be valid, it is important that the quasi-elastic final state be well identified as a single particle or resonance. A more complex final state like $`\gamma ^{}p\pi ^+\pi ^{}p`$ (which is what is actually observed in the experiments) can proceed through several different channels (with intermediate states $`\rho `$, $`\omega `$, $`f_0`$…). Its study would require the careful treatment of phases and interference effects, see for example Refs.Hagler:2002nh ; Hagler:2002nf . We shall assume that we have a single-channel process, and we use formula (16) to extract the $`b`$-dependent amplitude from the available data.

We shall now be more specific, specializing to the case of electroproduction of vector mesons $`\gamma ^{}pVp`$, where $`V`$ will in practice be a $`\rho `$ meson. The picture is the following. At high energy, the photon splits into a dipole of fixed size $`r`$ which scatters off the proton before recombining into a meson. Due to its high relative velocity, the dipole does not change size during the scattering. Accordingly, the scattering amplitude $`𝒜`$ is the convolution of the photon wave function $`\psi _\gamma ^{}`$, the meson wave-function $`\psi _V`$ and the dipole amplitude $`𝒜_\mathrm{d}`$ in the same way as in Eq.(8)

$$𝒜(y,Q,b)=d^2r\psi _\gamma ^{}^{}(r,Q)\psi _V(r)𝒜_\mathrm{d}(y,1/r,b),$$
(17)
the forward amplitude $`2𝒯_\mathrm{d}`$ appearing in (8) being replaced by the impact parameter dependent one $`𝒜_\mathrm{d}`$. Due to the fact that the product of the wave functions $`\psi _\gamma ^{}^{}\psi _V`$ is peaked around a typical dipole size depending only on the external transverse momentum scales, the scaling of the dipole amplitude $`𝒜_\mathrm{d}`$ (see Eqs.(11)) is transmitted to the photon-proton amplitude $`𝒜`$. As it is well-known from both empirical tests Nemchik:1994fp ; Nemchik:1997cw  and experimental Marage:2001sy  studies, the relevant distance scale at the photon-meson vertex is a combination of the photon virtuality $`Q`$ and of the meson mass $`M`$, namely $`1/\sqrt{Q^2+M^2}`$ instead of $`1/Q`$ as in the case of inclusive scattering. We then obtain from Eqs.(10), (14) and (11) that the amplitude $`𝒜`$ in Eq.(17) should be a function of the scaling variable $`\tau =(Q^2+M^2)(x/x_0)^\lambda `$ only, once the impact parameter $`b`$ is fixed.

To check this statement, we apply the transformation (16) to the HERA published data on electroproduction of $`\rho ^0`$ vector mesons. In the presently published experimental analysis Breitweg:1998nh ; Adloff:1999kg , the total cross section as well as the logarithmic $`t`$-slope $`B`$ are quoted. The data can then be represented by the parametrization

$$\frac{d\sigma }{dt}(Q,W)=B(Q,W)\sigma (Q,W)e^{B(Q,W)|t|},$$
(18)
up to $`|t|`$ of the order of $`|t_{\mathrm{max}}|=1\text{GeV}^2`$, which corresponds to the range in which data have been collected. The scattering amplitude at each impact parameter reads

$$\begin{array}{cc}\hfill 𝒜(y,Q,b)& =\sqrt{\pi }\left(_0^{|t_{\mathrm{max}}|}d(t)J_0(b\sqrt{|t|})\sqrt{B\sigma }e^{B|t|/2}+_{|t_{\mathrm{max}}|}^+\mathrm{}d(t)J_0(b\sqrt{|t|})\sqrt{\frac{d\sigma }{dt}}\right)\hfill \\ & \sqrt{\pi }\left(_0^+\mathrm{}d(t)J_0(b\sqrt{|t|})\sqrt{B\sigma }e^{B|t|/2}\right)=2\sqrt{\frac{\pi \sigma (Q,W)}{B(Q,W)}}e^{b^2/2B(Q,W)}.\hfill \end{array}$$
(19)
We take the total cross section from Breitweg:1998nh ; Adloff:1999kg . The kinematical range of these measurements is shown on Fig.1. The experimental points range from $`Q^2=0.47`$ to $`27\text{GeV}^2`$ and $`W=23.4`$ to $`150\text{GeV}`$. Due to the lack of data for the slope $`B`$ for all values of $`Q^2`$ and $`W`$, we shall assume

$$B(Q,W)=0.6\left(\frac{14}{(Q^2+M^2)^{0.26}}+1\right)+4\alpha ^{}\mathrm{log}(W/75),$$
(20)
all quantities appearing in the previous formula being expressed in powers of 1 GeV. The parametrization (20) at $`W=75\text{GeV}`$ is taken from Caldwell:2001ky . A logarithmic energy dependence has been added according to the Donnachie-Landshoff parametrization for the Pomeron Donnachie:1984xq . We set the Regge slope $`\alpha ^{}`$ to the value $`0.25`$ which is measured in hadron-hadron cross sections. The relevance of this choice and its implications will be discussed in the next section.

The results of our analysis are displayed on Fig.2 for three different impact parameters, $`b=0.3,1.0\text{and}1.5\text{fm}`$. Needless to say, our plot is only illustrative of the method and, being theoreticians, we cannot provide a quantitative analysis. The error bars shown take into account only the error on the total cross section ($`\mathrm{\Delta }𝒜/𝒜=(\mathrm{\Delta }\sigma /\sigma )/2`$) and the contour gives a rough idea of the total uncertainty, evaluated as if they were uncorrelated and small by the formula

$$\frac{\mathrm{\Delta }𝒜}{𝒜}=\frac{1}{2}\frac{\mathrm{\Delta }\sigma }{\sigma }+\frac{1}{2}\frac{\mathrm{\Delta }B}{B}\left(1+\frac{b^2}{B}\right).$$
(21)
We assumed an error on the slope $`B`$ of 5% which seems realistic regarding the quality of recently measured data (see the $`J/\mathrm{\Psi }`$ analysis of Chekanov:2002xi ). We see from (21) that the uncertainty on $`B`$ affects mostly the large impact parameter region. For small $`b`$, an additional error would come from the fact that the data does not extend to arbitrarily large $`|t|`$, hence the exponential form (19) is not constrained for small $`b`$. This is the reason why we performed our analysis only down to impact parameters as large as 0.3 fm (see Munier:2001nr  for a more complete discussion).

The results on Fig.2 are consistent with “local” geometric scaling within uncertainties, for each of the displayed impact parameters. Although the kinematical range of the data is not very large, scaling is a non-trivial feature since data points of different $`Q`$ and $`W`$ overlap, as seen on Fig.1. We stress that from the technical point of view, the fact that the scaling holds for various values of $`b`$ comes from a subtle interplay between the $`(Q,W)`$ dependence of $`\sigma `$ and $`B`$ in formula (19). Finally, let us note that the data points denoted by a cross correspond to very low values of $`Q^2`$ ($`0.47\text{GeV}^2`$), and are kinematically quite separated from the other data. As they are measured with a different method, we do not make a point of the fact that they do not seem to lie on the same curve as the larger $`Q^2`$ data.

## IV Discussion and conclusion

In this paper, we have demonstrated that “local” geometric scaling, predicted by the saturation models for deep-inelastic scattering, is compatible with the HERA data. However, a consistent experimental analysis with careful treatment of uncertainties as well as an extension of the kinematical range is needed to draw a strong conclusion. In particular, an extension of the range in $`W`$ would enable us to compare measurements sitting at very different $`Q^2`$ but at the same scaling variable $`\tau `$. This would constitute a definite test of this local scaling.

Nevertheless, let us note the surprizing fact that scaling seems to work well for large impact parameters, very far from the center of the proton, and even outside the “grey” zone ($`b=1.5\text{fm}`$), where the concept of saturation scale makes little sense. This feature depends of course a lot on the form of the measured $`B`$, and particularly on its energy dependence, which was introduced in an ad hoc way in Eq.(20). Interestingly enough, if the energy dependence were less steep for large momentum transfers $`\sqrt{|t|}`$, scaling would also be slightly better verified at small impact parameters. Such a feature is likely to happen since it seems that the harder the momentum scales, the weaker the energy-dependence of the $`t`$-slope. This rule of thumb is most seen in large-$`t`$ proton-dissociative photoproduction of $`\rho `$ mesons Chekanov:2002rm , although of course it is not a quasi-elastic process and our analysis would not apply. Moreover, recent preliminary experimental results indicate that $`\alpha ^{}`$ is smaller than the Donnachie-Landshoff value, and decreases with the scale $`Q^2+M^2`$. If scaling were confirmed also for large impact parameters, it might be an indication of some more general symmetry of high energy dynamics.

One would in principle be able to plot all data for different $`b`$ as a function of a unique scaling variable $`Q^2/Q_s^2(y,b)`$, see formula (11). However, we checked that the parametrizations for the saturation scale proposed so far in Eqs.(5), (10) do not lead to a superposition of all points on a single curve. This means that the available approaches to the computation of the impact parameter dependence of the saturation scale are not good enough yet.

The theoretical status of the impact parameter dependence of the saturation scale is still unclear, mainly due to technical difficulties in solving the full saturation equations. Better insights into these issues would help a lot the understanding of high energy scattering within QCD. Important issues like unitarity problems require a good understanding of the evolution in impact parameter Kovner:2001bh ; Iancu:2002tr  (see Ref.Shoshi:2002in  for a model which predicts such an evolution). On the experimental side, diffractive vector meson production in DIS is being investigated intensively at HERA Levy:2003qz , and more data will be available soon. We believe that it would be worth to plot the forthcoming data in the way we propose in this paper.

###### Acknowledgements.

This work is funded by the European Commission IHP program under contract HPRN-CT-2000-00130. We thank H.-J. Pirner for a critical reading of the manuscript and useful suggestions. We are also grateful to A. Freund, E. Iancu and B. Pire whose remarks helped us to improve the manuscript.


text:
# Quantum dynamics in high codimension tilings: from quasiperiodicity to disorder

## Abstract

We analyze the spreading of wavepackets in two-dimensional quasiperiodic and random tilings as a function of their codimension, i.e. of their topological complexity. In the quasiperiodic case, we show that the diffusion exponent that characterizes the propagation decreases when the codimension increases and goes to 1/2 in the high codimension limit. By constrast, the exponent for the random tilings is independent of their codimension and also equals 1/2. This shows that, in high codimension, the quasiperiodicity is irrelevant and that the topological disorder leads in every case, to a diffusive regime, at least in the time scale investigated here.

It is now well established that quasiperiodic order has a strong influence on the quantum dynamics of wavepackets. Indeed, the nature of the eigenstates in quasiperiodic systems, which are neither spatially extended (as in periodic systems) nor localized (as in disordered systems) but critical, is often responsible for a sub-ballistic motion. Although most of the studies about this anomalous diffusion concern one-dimensional systems such as the Fibonacci or the Harper chain, there has also been a great interest for the, more physical, higher-dimensional ones Passaro\_Octo ; Zhong ; Yuan ; Triozon\_Rauzy ; Vidal\_ICQ7 . However, the parameters that determine the characteristics of the long time dynamics, such as the diffusion exponent $`\beta `$, remains misunderstood.

In this paper, we investigate the quantum dynamics of wavepackets in two-dimensional quasiperiodic tilings built with De Bruijn grid method deBruijn1 ; deBruijn2 . This construction allows us to change easily the codimension of the structures that determines the number of different tiles and thus, in a sense, its complexity. The codimension reveals itself as a fundamental quantity for the dynamics since the diffusion exponent decreases when the codimension increases and converges towards $`1/2`$ in the high codimension limit. We also analyze the influence of phason flips that consists in turning the quasiperiodic tiling into a random one without changing the stœchiometry of the tiles (see Fig. 1). We show that contrary to what was claimed in Ref. Passaro\_Octo  for the octagonal tiling, the phason disorder do not fasten the spreading but, always slow it down towards a diffusive regime.

The De Bruijn grid algorithm deBruijn1 ; deBruijn2  relies on two steps. First, we define $`D`$ families of $`N_l`$ equally spaced parallel lines, rotated one from the other by an angle $`2\pi /D`$. These $`D\times N_l`$ lines form a grid and separate different connected regions of the plane called cells of the decomposition. Second, we dualize this grid by associating a tiling vertex to each cell and connecting any vertices associated to two cells separated by a line.

This algorithm is simpler to implement in two dimensions than the standard Cut-and-Project KD ; KKL ; Elser\_random\_tiling  method where tilings are generated from high-dimensional spaces. In this latter approach, the codimension of the tiling is defined by the difference between the dimension $`D`$ of the initial lattice in which the points are selected, and the dimension $`d`$ of the tiling. For two-dimensional tilings, the Cut-and-Project method requires the analyzis of a $`(D2)`$-dimensional (perpendicular) space whose complexity increases very quickly with $`D`$ so that, in practice, this method is not tractable beyond $`D=5`$. By contrast, in the grid method, a tiling of codimension $`(D2)`$ is straightforwardly obtained by choosing $`D`$ families of lines. Note that for $`D=2`$, the grid method builds the usual square lattice, and, for $`D=3`$, it builds the well-known (periodic) dice lattice that displays interesting properties when embedded in a magnetic field Vidal\_Cages\_big . Nonperiodic (quasiperiodic) tilings are thus obtained for $`D4`$. The number of lines $`N_l`$ determines the size of the tiling since the number of sites is $`N_sN_l^2D(D1)/2`$. However, with this method, different tilings can be obtained depending on the relative position of each line families. Here, we have focussed on a generic family of tilings and checked that our results were weakly sensitive to this dephasing choice. In addition, to avoid spurious effects due to fixed boundary conditions Destainville  ($`N_l`$ is finite), we have isolated a circular central cluster of about $`10^6`$ sites, which is quasiperiodic.

The dynamics is given by a standard tight-binding Hamiltonian:

$$H=\underset{i,j}{}t_{ij}|ij|\text{,}$$
(1)
where $`|i`$ denotes a localized orbital on site $`i`$ and where the hopping term $`t_{ij}`$ equals to 1 if $`i`$ and $`j`$ are nearest neighbors and $`0`$ otherwise. Of course, the time evolution of any wavepacket is directly obtained by diagonalizing $`H`$, but unfortunately, exact diagonalizations are restricted to rather small system size ($`10000`$ sites). Since we are interested in the long time behavior of the dynamics that requires large systems, we have used an approximate method, the Second Order Differencing Scheme Numerique , which consists, for a given initial state $`|\psi (0)`$, in writing:

$$|\psi (t+\mathrm{\Delta }t)=|\psi (t\mathrm{\Delta }t)2i\mathrm{\Delta }tH|\psi (t)\text{.}$$
(2)
Practically, a time step $`\mathrm{\Delta }t=0.05`$ is sufficient to get a good accuracy and all our results have been obtained with this value. Note that despite the very low order of this development for the evolution operator $`e^{iHt}`$, this algorithm is actually an very efficient tool (see for example Vidal\_ICQ8 ).

The time evolution of a wavepacket can be characterized by several observables. Here, we focus on the mean square spreading defined by:

$$\mathrm{\Delta }R^2(t)=\psi (t)|\widehat{R}^2|\psi (t)[\psi (t)|\widehat{R}|\psi (t)]^2\text{,}$$
(3)
where $`\widehat{R}`$ is the position operator, the origin being taken at the center of the cluster. The diffusion exponent $`\beta `$ is defined by the long time behavior of $`\mathrm{\Delta }R^2(t)t^{2\beta }`$. A priori, nothing prevents from getting a different asymptotic regime but, actually, $`\mathrm{\Delta }R^2`$ often seems to be rather well described by a power-law at large time. The case $`\beta =1`$ obtained for periodic potential defines a ballistic propagation whereas $`0<\beta <1`$, observed in many quasiperiodic systems, corresponds to a sub-ballistic spreading. Note that the diffusive motion ($`\beta =1/2`$) is obtained in several situations among which the three-dimensional Anderson model in the metallic part of the spectrum. A crucial issue is to understand the parameters which determine the value of this exponent.

As already noticed in the octagonal tiling Passaro\_Octo  for wavepackets initially localized on a single site, $`\beta `$ strongly depends on this site. In the Labyrinth tiling Yuan  and in the generalized Rauzy tilings Triozon\_Rauzy  where energy-filtered wavepackets have been studied, $`\beta `$ has also been shown to be energy-dependent. In the present work, we do not claim to give a precise characterization of the diffusion exponent for each tiling but we wish to investigate the behavior of $`\beta `$ with respect to the codimension. Therefore, we have chosen to consider random phase states initially spread over a finite portion of the cluster (typically $`1\%`$). This leads to an exponent which is “averaged” over the total density of states. We have considered several configurations of the phases and several initial radius, and checked that although $`\beta `$ is different from one state to another, it fluctuates of less than $`5\%`$. In addition, since we have finite size systems with open boundary conditions, the propagation must be stopped when the wavepacket reach the boundaries. As a criterion for this maximum time $`t_{\text{max.}}`$, we have chosen the time for which the presence probability integrated over all the border sites reaches $`1/N_s`$. Typically, this leads to $`t_{\text{max.}}300`$. As usual, we can never ensure that the real asymptotic regime (if it exists !) is reached but this time scale certainly gives some hints concerning the long time behavior.

The mean square spreading $`\mathrm{\Delta }R^2(t)`$ is displayed in Fig. 2 for quasiperiodic tilings with $`D[3,6]`$. The dice lattice $`D=3`$ is given as a reference for a ballistic motion ($`\beta =1`$). For low codimension tilings, we observe a superdiffusive regime with an exponent that decreases when the codimension increases.

More precisely, we get $`\beta =0.77(1)`$ for the octagonal tiling $`(D=4)`$, $`\beta =0.71(1)`$ for the Penrose tiling $`(D=5)`$, and $`\beta =0.65(1)`$ for the dodecagonal tiling $`(D=6)`$. Nevertheless, we emphasize that even for these codimensions, the value of the diffusion exponent is not constant over a wide time range and actually slightly decreases when the time increases. This leads us to conclude that either the asymptotic regime is not reached, or the long time behavior is not “strictly” power-law like.

For higher codimension, the situation is more complicated: no power-law regime can be clearly distinguished since $`\beta `$ strongly depends on the time range and decreases when the time increases.

To determine whether it converges towards a asymtotic value (at large time) or not, we have focussed on the behavior of $`\mathrm{\Delta }R^2(t)/t`$ which gives us an indication on the super- or sub-diffusive character of the propagation. This is an important issue since in such two-dimensional quasiperiodic systems, no sub-diffusive propagation has ever been reported. As it can be seen in Fig. 3, although no power-law behavior can be extracted in the time range investigated, a strong indication of a diffusive asymptotic regime ($`\underset{t\mathrm{}}{lim}\mathrm{\Delta }R^2(t)/t=\text{Cte}`$) is provided.

These results show that the codimension is an important parameter for the wavepacket propagation. Indeed, noting that for the two-dimensional generalized Rauzy tiling Vidal\_ICQ8  $`(D=3)`$ one obtains an exponent $`\beta =0.95(1)`$, it is clear that the higher the codimension, the lower $`\beta `$. This can be easily understood in terms of the complexity of the structure for the following reasons. An important characteristic of quasiperiodic tilings is the so-called repetitivity property initially studied in two dimensions by Conway. In the class of tilings considered here, it states that for any local environment of linear typical size $`L`$, a similar environment can be found at a distance $`\xi L^z`$ with $`z1`$. In the simplest case $`D=4,5`$ one has $`z=1`$. Moreover, when the codimension increases, the number of possible edge orientations at each node is lower or equal to $`2D`$ and for a given $`L`$, $`\xi (L)`$ also grows. Consequently, the specific effects of the quasiperiodic order are expected to decrease when the codimension increases. Note that in periodic system, $`\xi `$ is a constant and in disordered systems, it grows exponentially with $`L`$.

Nevertheless, these tilings are quasiperiodic and this special order is responsible, at least in low codimension, for a superdiffusive regime ($`\beta >1/2`$). For higher codimension the quasiperiodicity becomes irrelevant and the exponent $`\beta `$ seems to converge towards $`1/2`$ which is the value expected for a weakly disordered system in two dimensions (for example a periodic lattice with random impurities or random magnetic fluxes Kawarabayashi , or a disordered quasiperiodic Fibonacci lattice Roche\_Fibo\_desordre ) where such a regime is expected when $`\sqrt{\mathrm{\Delta }R^2}`$ is lower than the localization length. Unfortunately, since the asymptotic regime is not reached with the system size studied here, we cannot estimate the codimension dependence of the diffusion constant. However, the results shown in Fig. 3 indicate that this constant does not seem to simply decrease when the codimension increases.

If this analysis is correct, destroying quasiperiodicity in these tilings should also lead to a diffusive regime. To check this assumption and to corroborate our previous statement, we have studied the quantum dynamics in random tilings obtained from the quasiperiodic tilings by making elementary phason flips (about $`10^{11}`$) as shown in Fig. 1 for $`D=4`$. The random flip Markov chain has been studied into detail in the cases $`D=3`$ and $`D=4`$ (see Destain02 ), and the number of flips needed to get a random tiling uncorrelated to the initial quasiperiodic tiling is indeed of order $`N_s^2`$. Furthermore, this number has been computed numerically in higher codimension tilings and is, in fact, always smaller than $`N_s^2/2`$ WDMB . In the present case, this leads to an upper bound of $`5.10^{11}`$ flips. However, we have checked for $`D=4`$, that the spreading in tilings with $`10^{11}`$ and $`5.10^{11}`$ flips have essentially the same diffusion exponent.

As shown in Fig. 4, although no power-law regime can be clearly extracted, the propagation is very similar to those displayed in Fig. 3 for high codimension quasiperiodic tilings. For the octagonal tiling ($`D=4`$), this result refutes the conclusion of Passaro et al. claiming that the diffusion exponent increases when a phason disorder is introduced. This contradiction is due to the short time scale investigated in Ref. Passaro\_Octo . Indeed, we also observe that for $`t20`$, the propagation is faster in the $`D=4`$ random tiling than in the $`D=4`$ quasiperiodic tiling. However, at larger times, the motion in the random tiling seems to become diffusive whereas, in the octagonal tiling, it remains superdiffusive with $`\beta =0.77(1)`$ ZZZ . We have also observed a similar crossover for higher codimension $`(D=5,6)`$ (see Fig. 5).

A possible explanation of these two regimes is that, at short times, the local quasiperiodic order is likely responsible for coherent quantum interferences (absent in a random tiling) that slow down the dynamics. In this regime, the topological disorder thus favors, in a sense, the propagation. At large times, the influence of the quasiperiodicity remains the same but a weak localization process becomes dominant for the random tiling and leads to a diffusive regime. Of course, for larger codimension, the situation is different since the asymptotic regime is likely the same both in the quasiperiodic and in the random case. Thus, one expects $`\mu =\mathrm{\Delta }R_{\text{quasi.}}^2/\mathrm{\Delta }R_{\text{random}}^2`$ to become constant at large time and it is actually the case (see Fig. 5 for $`D=10`$).

This study raise interesting questions concerning the nature of the spectrum in quasiperiodic and random tilings. A simple way to characterize the spectral measure of a state $`|\psi `$ is to compute the averaged autocorrelation function

$$C(t)=\frac{1}{t}_0^t𝑑t^{}P(t^{})t^\alpha \text{,}$$
(4)
where $`P(t)=|\psi (0)|\psi (t)|^2`$. If $`\alpha =1`$, the spectral measure is absolutely continuous, if $`\alpha =0`$, it is pure point, and if $`0<\alpha <1`$ it is singular continuous. As previously, for a random phase state, $`\alpha `$ gives an average information over the total density of states. Unfortunately, in the tilings studied here, there exist strictly localized states known as confined states due to specific local environment (see Ref. Kohmoto\_Sutherland ; Arai2  for the case $`D=5`$). These states eventually leads to $`\alpha =0`$ provided the initial random phase state has a nonzero overlap with them. The analyzis of the spectral measure with respect to the codimension thus requires to get rid of these confined states. Further investigations are needed to clarify this point which is certainly of great interest in the understanding of the electronic properties of quasiperiodic systems.

###### Acknowledgements.

We would like to thank M. Baake and M. Duneau for fruitful remarks about the repetitivity property of quasiperiodic tilings.


text:
# Percolation, Morphogenesis, and Burgers Dynamics in Blood Vessels Formation

## Abstract

Experiments of in vitro formation of blood vessels show that cells randomly spread on a gel matrix autonomously organize to form a connected vascular network. We propose a simple model which reproduces many features of the biological system. We show that both the model and the real system exhibit a fractal behavior at small scales, due to the process of migration and dynamical aggregation, followed at large scale by a random percolation behavior due to the coalescence of aggregates. The results are in good agreement with the analysis performed on the experimental data.

The problem of morphogenesis and self-organization in biological systems is both an actual and a long standing one yates88 ; jensen00 . In living beings, complex structures arise from the aggregation of separate constituents. Such structures, although irregular, often show some sort of hidden order, like self-similarity and scaling laws mandelbrot88 ; stanley2 . Scaling laws are known to emerge not only in the physics of phase transition stanley  and in percolation sa94 , but also from several kinds of aggregation dynamics, and are often considered as fingerprints of the process which led to the formation of the structure itself vicsek92 ; vicsek01 . A typical example of natural structure characterized by non trivial scaling laws is the vascular network wbe . In recent years many experimental investigations have been performed on the mechanism of blood vessel formation carmeliet00 . Cells are cultured on a gel matrix and their migration and aggregation observed through videomicroscopy. Tracking of individual trajectories shows marked persistence in the direction, with a small random component superimposed bio . The motion is directed towards zone of higher concentration of cells, suggesting that chemotactic factors play a role. Cells migrate over distances which are an order of magnitude larger than their radius and aggregate when they get in touch with one of their neighbors. In a time of the order of 10 hours they form a continuous multicellular network which can be described as a collection of nodes connected by chords. The mean chord length is approximately independent on the initial cell density $`\overline{n}`$, with an average value $`\overline{\mathrm{}}=200\pm 20\mu \mathrm{m}`$ for $`\overline{n}`$ between 100 to 200 cells/mm2. Observing the formed structures, one finds that, by varying the initial cell density, there is a percolative transition. Below a critical value $`n_\mathrm{c}100`$ cells/mm2 groups of disconnected structures are formed (Fig. 1a). For $`\overline{n}n_\mathrm{c}`$ a single connected network is instead visible (Fig. 1b,c). For higher values of $`\overline{n}`$ one observes a sort of “swiss cheese” pattern (Fig. 1d). In the biological system the percolating property is of physiological relevance, since it is directly linked with the functionality of blood vessels.

In this paper we propose a theoretical model which turns to be in good agreement with experimental observations. It allows to well reproduce both the observed percolative transition and the typical scale of observed vascular networks. With respect to standard percolative models, here an important role is played by migration and dynamical aggregation of particles. The model appears thus as a possibly new and physically interesting representative in the class of percolation models. We characterize the phase transition from the point of view of scaling laws, both in the model and in the real system, computing critical indices and the fractal dimension of the percolating cluster through extensive numerical simulations and analysis of the experimental data. The model describes the cell population as a continuous density field $`n`$ and velocity $`𝐯`$; it also assumes the presence of a concentration field $`c`$ of soluble factor. Cells are modeled as a fluid accelerated by gradients of the soluble factor. The latter is supposed to be released by cells, diffuse, and degrade in finite time, in agreement with experimental observations. These assumptions give rise to the following equations:

$`{\displaystyle \frac{n}{t}}+\left(n𝐯\right)=0`$ (1a)
$`{\displaystyle \frac{𝐯}{t}}+𝐯𝐯=\mu c`$ (1b)
$`{\displaystyle \frac{c}{t}}=D^2c+\alpha n\tau ^1c`$ (1c)
where $`D,\alpha ,\tau `$, are respectively the diffusion coefficient, the rate of release and the characteristic degradation time of soluble mediators, and $`\mu `$ measures the strength of cell response. Initial conditions are given in the form of a set of randomly distributed bell-shaped bumps in the density field, representing single cells initially spread in the system, with zero velocities and zero concentration of the soluble factor. In Eqs. (1) a multidimensional Riemann (1b) equation is coupled to the diffusion equation (1c), while (1a) is just a continuity equation expressing conservation of the total cell number. Multivaluedness of solutions to (1b) can be prevented by adding a viscous term $`^2𝐯`$ in the r.h.s. One than gets the analog of Burgers’ equation, which has been used to describe growth of molecular interfaces kpz86  and the emergence of network-like patterns bik02  in the large-scale distribution of masses of the Universe z70 ; vdfn94 . Instead of viscosity, we introduce a phenomenological, density dependent pressure term $`\varphi (n)`$, in the r.h.s. of (1b), with $`\varphi (n)`$ zero for low densities, and rapidly increasing for densities above a threshold $`\sigma ^2`$, where $`\sigma =30\mu \mathrm{m}`$ is the mean cell radius, in order to describe the fact that cells do not interpenetrate.

The model provides a natural lenght scale $`r_0=\sqrt{D\tau }`$, which is the effective range of the interaction mediated by the soluble factor. This suggests that, starting from the above mentioned initial conditions, Eqs. (1) should develop network patterns characterized by the same scale. Direct measurements of $`D`$ and $`\tau `$ give $`D10^7`$ cm2/s, $`\tau `$ 1h, from which one gets $`r_0`$ 200 $`\mu `$m, a value which is in good agreement with the mean chord length $`\overline{\mathrm{}}`$ experimentally observed.

We performed numerical simulations of model (1) on square boxes of size $`L=1`$, 2, 4, $`8\mathrm{mm}`$, with periodic boundary conditions, using a finite volume method leveque . We used the experimental values for $`D`$ and $`\tau `$, while the unknown parameters $`\mu ,\alpha `$ were set to one (which amounts to an appropriate rescaling of time and the concentration field $`c`$). The initial conditions were given by throwing the same number of cells as in the biological experiments in random positions inside the box, with zero velocities and zero concentration of the soluble factor, with a single cell given initially by a Gaussian bump of width $`\sigma `$ and unitary weight in the integrated cell density field $`n`$. Starting from these condition, Eqs. (1) were numerically integrated. The simulation was stopped when the vascular network was formed, or a stationary state was reached. We observed a remarkable correspondence between simulations and experiments both in the the dynamic evolution details (not shown) and in the reproduction of the percolating patterns (Fig. 2).

In order to study the percolative transition we computed between 100 to 200 realizations of the system for each box size and cell density, with different initial positions of the cells. The simulation gives as a result a continuous density field $`n`$. In order to study percolation properties we partitioned the box in little square boxes of size $`2^6`$ mm, forming the sites of a square lattice, and set each site as filled if the mean cell density inside the corresponding little box was larger than the threshold $`\sigma ^2`$, empty otherwise. We then identified clusters of nearest neighbor filled sites. In Fig. 2 the results of four simulations, for different values of the density, are shown. Black areas represent regions filled with cells, that is regions where density exceeds the threshold, white areas represent the underlying substrate. As in the case of experimental data, a percolative transition is observed at some critical cell density $`n_\mathrm{c}`$. Furthermore the model reproduces quite well the typical structure of the vascular network, with chords of length $`\overline{\mathrm{}}`$ 200 $`\mu `$m for a wide range of cell densities, in good agreement with the expected value $`r_0`$.

For each set of realizations, corresponding to a value of the box size $`L`$ and of the mean cell density $`\overline{n}`$, we identified clusters of nearest neighbor filled sites, and measured the following quantities: 1) the percolation probability $`\mathrm{\Pi }`$, defined as the fraction of realizations in which a percolating cluster appears; 2) the infinite cluster density $`P`$, defined as the density of sites belonging to the infinite (greatest) cluster; 3) the mean cluster size $`S=(1/N)_sn_ss^2`$, where $`N`$ is the total number of sites and $`n_s`$ is the number of clusters of size $`s`$, excluding from the sum the infinite cluster.

In percolation models, in presence of a second order transition at some mean cell density $`n_\mathrm{c}`$, these quantities show a characteristic singular behavior near the percolation point. Namely, the percolation probability is zero below the transition and one above it; the infinite cluster density, which plays the role of order parameter, vanishes as $`|\overline{n}n_\mathrm{c}|^\beta `$ from above; and the mean cluster size has a singularity $`|\overline{n}n_\mathrm{c}|^\gamma `$ on both sides of the transition. For finite size systems, the same quantities obey near the transition to the following relations:

$`\mathrm{\Pi }(\overline{n},L)`$ $``$ $`\widehat{\mathrm{\Pi }}[(\overline{n}n_\mathrm{c})L^{1/\nu }]`$ (2a)
$`P(\overline{n},L)`$ $``$ $`L^{\beta /\nu }\widehat{P}[(\overline{n}n_\mathrm{c})L^{1/\nu }]`$ (2b)
$`S(\overline{n},L)`$ $``$ $`L^{\gamma /\nu }\widehat{S}[(\overline{n}n_\mathrm{c})L^{1/\nu }]`$ (2c)
Using Eq. (2a), one can find the critical density $`n_\mathrm{c}`$ as the point where the curves $`\mathrm{\Pi }(\overline{n},L)`$ for different sizes $`L`$ cross. The estimated value is $`n_\mathrm{c}=94\pm 1`$ cells/mm2. Then, the critical index $`\nu `$ is found by plotting $`\mathrm{\Pi }(\overline{n},L)`$ as a function of $`(\overline{n}n_\mathrm{c})L^{1/\nu }`$, where $`n_\mathrm{c}`$ is held fixed to the value previously computed, and looking for the value of $`\nu `$ that gives the best data collapse. The error on the computed value of the index is estimated looking for the range of values that gives an acceptable data collapse. Holding the value of $`n_\mathrm{c}`$ and of $`\nu `$ fixed, the other two relations Eq. (2b) and (2c) are used to estimate the value of the indices $`\beta `$ and $`\gamma `$, by looking for the data collapse of $`L^{\beta /\nu }P(\overline{n},L)`$ and $`L^{\gamma /\nu }S(\overline{n},L)`$ as a function of $`(\overline{n}n_\mathrm{c})L^{1/\nu }`$. We observed that while $`n_\mathrm{c}`$ is sensitive to the choice of the threshold in a neighborhood of the natural value $`\sigma ^2`$, the critical exponents do not depend on it. In Fig. 3 the data collapses for $`\mathrm{\Pi }`$ and $`S`$ are shown. The estimated values of the critical indices are reported in the second column of Tab. 1.

Within the errors, the values obtained are the same that one expects for the simple random percolation model, that is when the sites of the lattice are uncorrelated, and each one occupied with the same probability. As critical indices are linked to the large scale structure of the percolating cluster, this means that, on such a large scale, the structure of the vascular network is mainly determined by the initial random positioning of the cells, and is not altered by the dynamical process of migration and aggregation. On the other hand, on smaller distances clusters are quite different from the ones of random percolation. For example the mean fraction of occupied sites, which is a quantity sensible to the small scale structure of clusters, is about 0.2 at the percolation point, much lower than the corresponding figure for random percolation, that is 0.59. A quantity that can give us informations about the structure of the percolating cluster at different scales is the density $`\rho (r)`$ as a function of the radius. This is defined as the mean density of sites belonging to the percolating cluster, inclosed in a circle of radius $`r`$ centered at one site belonging to the cluster, and averaged over different centering sites and different realizations of the system. For a fractal object with fractal dimension $`D`$, this should scale as $`\rho (r)r^{Dd}`$. For the percolating cluster of random percolation at the critical point, one expects a fractal dimension $`D=1.896`$. We computed $`\rho (r)`$ for our model at the percolation point $`\overline{n}=n_\mathrm{c}`$ (see Fig. 4), and found two different regimes depending on the scale. For $`r>r_\mathrm{c}`$, with $`r_\mathrm{c}=0.77\pm 0.08`$ mm, a behavior compatible with random percolation is found, with fractal dimension $`D=1.87\pm 0.03`$. On the other hand, for $`r<r_\mathrm{c}`$ a different exponent is found, $`D=1.50\pm 0.02`$. This exponent may be the signature of the dynamic growth process driven for $`r<r_\mathrm{c}`$ by the rapidly oscillating components of the concentration field.

We have repeated the measures described above on a set of experimental data, consisting of 28 digital photographs (with a resolution of 10242 pixels) of mature structures obtained on a matrix of size $`L=2`$ mm, starting from initial densities of 50, 75, 100, $`\mathrm{}`$, 200 cells/mm2. The photographs are the results of two experimental sessions, each performed in duplicate, so that 4 photographs are available for each density value. The data show the presence of a percolative transition at $`n_\mathrm{c}125`$ cells/mm2, and are compatible with the values $`\gamma /\nu =1.83`$ obtained from numerical simulations. The available data are not enough to obtain $`\nu `$ and $`\beta /\nu `$ with a reasonable precision. Estimated values of the critical indices are reported in the third column of Tab. 1. Finally, we measured the density of the percolating cluster as a function of the radius at the percolation point, and found a behavior surprisingly similar to the one found in the numerical model (see Fig. 4): for $`r>r_\mathrm{c}`$ the fractal dimension of the percolating cluster is $`D=1.85\pm 0.10`$ which is compatible with $`D=1.87\pm 0.03`$ from the numerical simulations and with the value $`D=1.896`$ of random percolation. For $`r<r_\mathrm{c}`$ instead $`D=1.48\pm 0.05`$.

In conclusion, we have introduced a model which describes the formation of the vascular network with a typical scale lenght $`r_0`$, and is in good agreement with experimental data. The percolative phase transition appears to be second order, at least to the accuracy allowed by the limitations in size which are intrinsic to our numerical experiments. Comparison of critical indices suggests that the phase transition falls in the universality class of random percolation, even in presence of migration and dynamical aggregation. This is confirmed by the fractal dimension of the percolating cluster, on scales larger than $`r_\mathrm{c}0.8`$ mm, that is in agreement with the one expected for random percolation. On the other hand we found that, if observed on scales smaller than $`r_\mathrm{c}`$, the vascular network shows a different scaling behavior, with a fractal dimension $`D1.5`$, that may be the signature of the dynamical process that led to the formation of the network.

Thus, the model appears to be quite successful in describing in vitro experiments, where all the parameters are under control and one can easily tune the cell density. It would also be interesting to understand the relevance of this picture to vascular networks in living beings. It is known carmeliet00  that during embriogenesis blood vessels form through local (in situ) differentiation of a fraction of endothelial cells, that afterwards assemble in a vascular labyrinth, a process similar to the one described in this paper. However, in this case it is more difficult to have complete control over all the relevant parameters. In particular, one would expect that the cell density is fixed through some feedback mechanism, in order to maximize efficiency wbe . Indeed, some experiments gblb95  on in vivo capillaries show compact space-filling structures, with fractal dimension $`D=2`$, which would indicate that cell density is tuned by Nature to a value larger than the percolation threshold. More experiments, focusing on the relevance of the percolation mechanism, and the nature of cell density self-regulation, are needed to shed light on these issues.

Acknowledgments – The authors acknowledge useful discussions with N. Bellomo, I. Kolokolov, P. Netti and M. Vergassola. This work was partially supported by the European TMR Network-Fractals (Contract No. FMRXCT980183), RTN HPRN-CT-2000005, CNR 00.00141.ST74, MURST-PRIN-2000, MURST-Cofin-2001, INFM-PRA (HOP), and MIUR-FIRB-2002.


text:
# Comments on “Improved limit on quantum-spacetime modifications of Lorentz symmetry from observations of gamma-ray blazars”

## I Introduction

In a recent paper Crab  we presented a new constraint on Lorentz violation of order $`E/M_{\mathrm{Planck}}`$ in the electron dispersion relation, based on the observation of synchrotron radiation from the Crab nebula. This constraint improved by a factor of one billion the previous best constraint, which was based on the observation of photon absorption from distant blazars SG ; Short ; Long  on the infrared background radiation. (The latter constraint limits the relevant Lorentz violating parameters to be of order unity in Planck units.) Both of these analyses were criticized in a recent paper by Amelino-Camelia Ref. GAC . The purpose of this note is to respond to these comments, clarifying our strategy and correcting some points that have been misconstrued in GAC .

## II The synchrotron constraint

### II.1 Theoretical Framework of Effective Field Theory

Ref. GAC  observes that the synchrotron constraint of Ref. Crab  relies on assumptions about the dynamics of electrons and electromagnetic fields rather than just the kinematical dispersion relations like some other constraints. It is stated there that if dynamical assumptions are made then “we are back to our starting point, we are actually proposing a full quantum-gravity theory, with all the uncertainties and risks of inconsistencies that plague quantum-gravity research.” It is true that constraints derived using kinematics and dynamics are qualitatively different than constraints derived from kinematics alone. However, we do not agree that our assumption about low energy dynamics is equivalent to proposing a full quantum gravity theory.

The essential assumption in Ref. Crab  about the full quantum gravity theory is that at energies low compared to the Planck energy electrons and the electromagnetic field can be described by an effective field theory preserving gauge invariance and rotation invariance. The only aspect of the interactions in this effective field theory that we use is that the modifications to the usual Lorentz-invariant QED equations of motion are suppressed by an inverse power of the Planck mass.

Let us elaborate a bit on how we view the significance of the effective field theory assumption. At low energies, physics is well described by the standard model and general relativity, which are both believed to be effective field theories. Since the particle energies even in astrophysical observations are far below the Planck energy it is reasonable to simply extend the low energy framework that we know is accurate. Ref. GAC  comments that the dynamics we adopt “appear to be consistent with a classical and continuous spacetime, while most authors would expect deformed kinematics at the Planck scale to be the result of non-classical (discrete, noncommutative,…) aspects of spacetime structure, which should have equally dramatic (but presently unknown) consequences for dynamics.” Exotic possibilities at the Planck scale do not preclude a low energy effective field theory description however. Some examples of this are lattice field theories, effective field theories of the low energy degrees of freedom of condensed matter systems 111Ref. GAC  states that the superluminal dispersion relations considered in Refs. Short ; Long  are “conceptually disfavored” because ordinary media always lead to subluminal velocities. However, there exist condensed matter systems, like Bose-Einstein condensates (see e.g. Liberati ) or superfluid 3He-A (see e.g. Volo ), where Lorentz invariance emerges at low energies for quasiparticles while superluminal dispersion characterizes the high energy propagation., and noncommutative field theory Carroll:2001ws . While we feel that the effective field theory assumption is quite reasonable, it is certainly conceivable that quantum gravity does not satisfy it. The constraints derived in Ref. Crab  apply only modulo the effective field theory assumption.

The starting point for an analysis of possible quantum gravity effects (such as Lorentz symmetry violation) using effective field theory is to add to ordinary field theory the operators that realize the effects. One can then characterize the observational consequences. This approach has been extensively pursued in the case of renormalizable field theory Kostelecky ; Kbook .

In the synchrotron case, the effective field theory being modified is QED. Since there are suggestions from quantum gravity that Lorentz symmetry might be violated with a suppression by $`M_{\mathrm{Planck}}^1`$, we introduce in the low energy effective theory dimension five operators that violate Lorentz symmetry. We assume this is the only standard symmetry that is broken, so in particular we assume both rotation and gauge invariance are preserved. 222 One could of course imagine that more symmetries are broken, but it would take a conspiracy for these to produce the usual symmetric physics when one broken symmetry alone would not. Hence we regard it as a fruitful strategy to begin by constraining (or looking for) minimal deviations from standard physics . The dispersion relations for free electrons and photons determine most of the terms in the effective field theory. In addition to momentum dependence, there can be polarization and chirality dependence of the dispersion. These possibilities correspond to different choices of the effective field theory. It was recently shown in Myers-Pospelov  that there are only three additional terms quadratic in the fields that contribute to the dispersion relation $`E^2(p)`$ at cubic order in the momentum. In particular, the deformation parameters for left and right circular polarized photons must be negatives of each other, while the parameters for each electron chirality are independent. For the synchrotron constraint, the polarization dependence is irrelevant in the interesting region of parameter space. We assumed in our synchrotron calculation that the Lorentz violating parameters for electrons are chirality independent. It may be possible to lift this restriction and still get powerful constraints.

The interaction terms in the effective field theory arise from minimal coupling (replacement of $``$ by $`+ieA`$) and non-minimal coupling in dimension five operators of the form $`\overline{\psi }\mathrm{\Gamma }^{\alpha \beta }\psi F_{\alpha \beta }`$, where $`\mathrm{\Gamma }^{\alpha \beta }`$ is built from gamma matrices and possibly a factor of the preferred timelike direction $`u^\alpha `$. All terms other than the usual minimal coupling are suppressed by a factor of $`1/M_{\mathrm{Planck}}`$ hence can be neglected for our purposes, since the unsuppressed Lorentz invariant interaction terms dominate.

### II.2 Objections to derivation of synchrotron constraint

There are three main theoretical objections in Ref. GAC  to our derivation of synchrotron radiation in the presence of Lorentz violation. We respond to these objections below. First, however, we wish to address the criticism that the observation of Crab synchrotron emission is “at best a promising conjecture”. We believe it is actually much more than a conjecture. The interpretation of the two observed humps in the emission spectrum from Crab and other supernova remnants as due to synchrotron and inverse Compton (IC) is the working hypothesis in the literature (see e.g. AA96 Koyama ). This hypothesis is able to explain the observed fluxes with success at least up to 20 TeV, with a value of the magnetic field that is measured with consistent results by several methods. Uncertainties remain concerning the possible presence of more than one population of accelerated electrons and the specific mechanism responsible for the highest energy IC photons (e.g. synchrotron-self-Compton (SSC) plus IC on infrared and CMBR photons or hadronic contributions AA96 ; Tanimori ), but the synchrotron nature of the first hump is widely accepted. Thus we think that our adoption of the “standard model” for the Crab emission is well-justified. We now turn to the theoretical objections in Ref. GAC  about the derivation of the synchrotron constraint.

#### II.2.1 Validity of heuristic formula for the cutoff frequency of synchrotron radiation

Ref. GAC  states that we assume the Lorentz invariant expression (Eqn. (4) in Crab ) for the synchrotron cutoff frequency $`\omega _c(E)`$ is valid in the Lorentz violating case. This is not an assumption—it is straightforward to verify explicitly and was done. As we state in our paper, this is a purely kinematical result. It involves the radius of curvature $`R(E)`$ of the electron trajectory, the angular width $`\delta (E)`$ of the synchrotron beam, and the group velocities of the electron and light, all of which are purely kinematical quantities (although the energy dependence of the first two depends on the dynamics).

#### II.2.2 Value of electron path radius of curvature $`R(E)`$

Ref. GAC  comments that we assume $`R(E)`$ is equal to its Lorentz invariant value, and that support for this assumption is given by proposing a new dynamics. We are not proposing a new, arbitrary dynamics, but merely implementing the dynamics given by the effective field theory approach discussed above. Since the modification of the electron field equation is strongly suppressed, $`R(E)`$ remains almost unchanged. More explicitly, gauge invariance determines the leading order interaction of the electrons and electromagnetic field via minimal coupling as described above. We show by calculation using this minimal coupling that while $`R(E)`$ is not equivalent to its Lorentz invariant value, it varies by some small relative amount which can be neglected in computing the cutoff synchrotron frequency.

#### II.2.3 Value of the opening angle $`\delta (E)`$ of the synchrotron beam

Ref. GAC  states that the relation between the opening angle $`\delta (E)`$ and $`E`$ is assumed to be unchanged. In fact it is not assumed, but rather argued (albeit briefly) that the scaling $`\delta (E)\gamma ^1(E)`$ follows from the effective field theory. We expand on this argument here. For a fixed source term in the electromagnetic field equation, the effective field equation has a solution for the vector potential of the form $`A=A_{LI}+A_{dev}`$, where $`A_{LI}`$ is the field that would be produced by the same source using the standard Maxwell equations. The deviation $`A_{dev}`$ is a consequence of the Lorentz violation and contains a suppression factor of $`1/M_{\mathrm{Planck}}`$. Dimensional analysis indicates that it will thus be suppressed by a factor $`\omega /M_{\mathrm{Planck}}`$ where $`\omega `$ is the highest frequency in the problem. For the synchrotron emission from the Crab nebula, this frequency is roughly 1 GeV, leading to a suppression of $`A_{dev}`$ by a factor of $`10^{19}`$. This is not competitive with the Lorentz invariant term, hence can be neglected. Thus the angular distribution of the radiation from a given source will be to a very good approximation the same as it is in the Lorentz invariant case.

In Ref. GAC  evidence is given for large deviations of $`\delta (E)`$ by viewing the synchrotron process as an off-shell threshold phenomenon, and observing that the angular distribution in such a process can be very sensitive to Lorentz violation. However, it is a long way from such an observation for individual off-shell processes to a calculation of the classical, coherent effect of radiation from accelerating charges in a slowly varying magnetic field. It is clear from our field theoretic analysis that in the end the angle sensitivity discussed in Ref. GAC  has no impact on the opening angle $`\delta (E)`$. Unlike in threshold phenomena, there is a Lorentz invariant zeroth order contribution that always dominates the synchrotron emission.

## III $`\gamma `$-ray absorption with the infrared background

We turn now to the weaker constraints that can be derived from the absorption of high energy gamma rays from blazars on the cosmic infrared background. Ref. GAC  argues that the type of constraint derived in SG ; Short ; Long  is conditional on unverified assumptions about the source spectrum and IR background. We agree that there is some uncertainty here, although not as much as it is made out to be in GAC . Our reasoning was not very explicit in Long  however, so we elaborate here briefly on our viewpoint, taking the opportunity to strengthen the case somewhat.

Our starting point was the analysis carried out in reference SD . There it was shown that using the most accurate model available for the infrared background the reconstructed spectrum of Mkn 501 shows no sign of anomalous pile up. Moreover, it was shown recently Kono  that the SSC model accounts remarkably well for the intrinsic spectra of the blazars Mkn 501 and 421 (the latter in two different states of emission) consistently in both the synchrotron and IC regions, using the same IR background. Hence the statement in Ref. SG  that there is “no indication of Lorentz invariance breaking” up to 20 TeV is well justified.

The constraint in question corresponds to the statement that the soft photon threshold for absorption of a 20 TeV gamma ray should not be shifted upwards beyond 25 meV (50 $`\mu m`$), 333In Ref. GAC  the author reported an error in our paper Long  saying that 25 meV photons correspond to a wavelength of 8 $`\mu m`$ rather than 50 $`\mu m`$. This statement is not correct. the usual Lorentz invariant threshold ($`\omega _{th}=m^2/k`$) for a 10 TeV photon. To justify this statement in Ref. Long  we said that since there is no evidence for anomalies up to 10 TeV, it is unlikely that the threshold for 25 meV can be raised by more than a factor of order unity and remain consistent with the data. This was criticized in Ref. GAC  on the grounds that one cannot observationally confirm absorption of any given soft photon energies, since it is only the effect on the spectrum of hard photons that is observed. We believe our justification of the constraint was inadequate, but the constraint is nevertheless justified by the absence of anomalies in the reconstructed source spectrum out to 20 TeV rather than just out to 10 TeV. (This factor of two makes a significant difference because it appears cubed in the constraint.)

The grounds for such a constraint follow from the shape of the IR background spectrum reported in SD  If the absorption threshold for 20 TeV gamma rays were shifted up by a factor of two from 12.5 meV to 25 meV, that would eliminate all the absorption from the far infrared hump of the spectrum. This would lead to a sharp downturn in the reconstructed source spectrum above 10 TeV, which would be inconsistent with the SSC source model. (To be more precise about this effect it would be necessary to reconstruct the source spectrum allowing for Lorentz violation in the absorption on the IR background.)

## Acknowledgements

We are grateful to Giovanni Amelino-Camelia for stimulating criticism. This work was supported in part by the NSF under grant No. PHY98-00967.