text:
# Subdiffusion-Limited Reactions

## 1 Introduction

Diffusion-limited reactions in low-dimensional geometries have been studied intensely because they exhibit “anomalous kinetics,” that is, behavior different from that predicted by the laws of mass action in well-stirred systems . Among the simplest and most extensively studied are single species diffusion-limited coagulation ($`A+AA`$ or $`A+AA`$ and annihilation ($`A+A0`$. These reactions, which show anomalous behavior in one dimension, are of particular theoretical interest because they lend themselves to exact solution in one dimension. Exact one-dimensional solutions can be obtained for diffusion in continuous  as well as in discrete  systems. A particularly elegant solution of the coagulation problem is provided by the method of intervals. This method focuses on the diffusive evolution of empty intervals, that is, of intervals that contain no particles . The empty interval equation turns out to be linear and hence exactly solvable. Originally developed for continuous systems, the method has been extended to discrete lattices . The method of intervals can not be directly adapted to the annihilation reaction, but a new formalism, the method of odd/even intervals which keeps track of the parity (even or odd) of the number of particles in an interval, has recently been developed . Again, the odd/even interval equation turns out to be linear and hence exactly solvable.

The anomalies of the $`A+A`$ problem in one dimension are typically displayed in two ways: one is through the time dependence of the reactant concentration $`c(t)`$, which for the $`A+AA`$ and the $`A+A0`$ reactions in infinite systems decays asymptotically as $`t^{1/2}`$ instead of the law of mass action decay $`t^1`$. The other is through the interparticle distribution function $`p(x,t)`$, which is the (conditional) probability density for finding the nearest particle at a distance $`x`$ on one side of a given particle. This function scales as $`x/t^{1/2}`$, in typical diffusive fashion. In one dimension a gap develops around each particle that leads to a more ordered spatial distribution than the exponential distribution implicit in well-stirred systems and “explains” the relative slowing down of the reaction.

In a parallel development, the problem of “anomalous diffusion” has also attracted a great deal of attention . The universally accepted characterization of anomalous (as in “not ordinary”) diffusion is through the mean squared displacement of a process $`x(t)`$ for large $`t`$,

$$x^2(t)\frac{2K_\alpha }{\mathrm{\Gamma }(1+\alpha )}t^\alpha .$$
(1)
Ordinary diffusion ($`\alpha =1`$, $`K_1D`$) follows Gaussian statistics and Fick’s second law for the probability density for finding the process at $`x`$ at time $`t`$,

$$\frac{}{t}P(x,t)=D\frac{^2}{x^2}P(x,t),$$
(2)
leading to linear growth of $`x^2(t)`$ with time. Anomalous diffusion is characterized by a nonlinear dependence. If $`0<\alpha <1`$ the process is subdiffusive or dispersive; if $`\alpha >1`$ it is superdiffusive. Anomalous diffusion is associated with many physical systems and is not due to any single universal cause, but it is certainly ubiquitous. Nor is anomalous diffusion modeled in a universal way; among the more successful approaches to the subdiffusive problem, which is the case we consider in this work, have been continuous time random walks with non-Poissonian waiting time distributions , and fractional dynamics approaches in which the diffusion equation (2) is replaced by the generalized diffusion equation

$$\frac{}{t}P(x,t)=K_\alpha {}_{0}{}^{}D_{t}^{1\alpha }\frac{^2}{x^2}P(x,t)$$
(3)
where $`{}_{0}{}^{}D_{t}^{1\alpha }`$ is the Riemann-Liouville operator,

$${}_{0}{}^{}D_{t}^{1\alpha }P(x,t)=\frac{1}{\mathrm{\Gamma }(\alpha )}\frac{}{t}_0^t𝑑\tau \frac{P(x,\tau )}{(t\tau )^{1\alpha }},$$
(4)
and $`K_\alpha `$ is the generalized diffusion coefficient that appears in Eq. (1). Some limitations of this description and some connections between the generalized diffusion equation and continuous random walk formulations have been discussed recently .

In this paper we consider a combination of these two phenomena, namely, the kinetics of $`A+A`$ reactions of particles that move subdiffusively in one dimension . Some aspects of this problem have also been considered using the waiting time distribution approach . Those solutions require approximations relating the reactant concentration to the distinct number of sites visited by a particle , or the waiting time distributions for single particles to the waiting time distributions for relative motion . Here we adapt the fractional dynamics approach to the problem and take advantage of the fact that the resulting generalized diffusion equations can be solved in closed form . We consider both coagulation and annihilation reactions.

In Section 2 we develop the fractional diffusion equations from which one obtains the concentration of reactant as a function of time. For the coagulation reaction we generalize the empty interval method  and, for the annihilation reaction, the odd/even interval method . The solutions are presented in Section 3. Sections 2 and 3 deal with kinetics on continuum one-dimensional systems; in Section 4 we extend the analysis to reactions on a lattice. The differences between the continuum and discrete results are most pronounced at early times (long-time differences are due to finite size effects that get pushed to longer times if one deals with larger systems). We conclude with a recapitulation in Section 5.

## 2 Evolution equations for subdiffusive particles

### 2.1 Empty interval method for coagulation reaction

Consider first the coagulation reaction $`A+AA`$ when the particles move by ordinary diffusion. The probability distribution function for the position $`x`$ of any one $`A`$ particle in the absence of reaction obeys the diffusion equation (2). The coagulation problem can be formulated in terms of the probability $`E(x,t)`$ that an interval of length $`x`$ is empty of particles at time $`t`$. This “empty interval” function obeys the diffusion equation

$$\frac{}{t}E(x,t)=2D\frac{^2}{x^2}E(x,t).$$
(5)
The derivation of this equation is straightforward and recognizes that an empty interval is shortened or lengthened by movement of particles in and out at either end according to the dynamics described by Eq. (2).

It is instructive to recall this derivation as given in Refs. . The reasoning is more transparent if the diffusion processes that lead to the change $`E(x,t)/t`$ are described on a lattice with lattice spacing $`\mathrm{\Delta }x`$ and sites labeled by integers. In this lattice $`E_n(t)`$ gives the probability that sites $`1`$ through $`n`$, for example, are empty at time $`t`$. Since the event $`E_{n+1}`$ contains the event $`E_n`$, the probability that sites $`1`$ through $`n`$ are empty but that there is a particle at site $`n+1`$ is $`E_{n+1}E_n`$. The particles move randomly to the nearest lattice site with a hopping rate $`2D/(\mathrm{\Delta }x)^2`$, the rate $`D/(\mathrm{\Delta }x)^2`$ to the right being equal to the rate to the left. If the particle at $`n+1`$ moves to site $`n`$ in a short time interval $`\mathrm{\Delta }t`$, $`E_n`$ will decrease by $`[D/(\mathrm{\Delta }x)^2](E_nE_{n+1})\mathrm{\Delta }t`$. On the other hand, $`E_n`$ may increase by the departure of a particle at site $`n`$ to site $`n+1`$; the associated increase is $`[D/(\mathrm{\Delta }x)^2](E_{n1}E_n)\mathrm{\Delta }t`$. The same entry and exit processes can take place at the other end of the interval. Combining these four processes then gives for the change in $`E_n`$

$$\frac{\mathrm{\Delta }E_n}{\mathrm{\Delta }t}=\frac{2D}{(\mathrm{\Delta }x)^2}(E_{n+1}2E_n+E_{n1}).$$
(6)
The coagulation reaction fixes the boundary condition. Coagulation involves adjacent occupied sites and either of the two particles hopping onto the other, thereby disappearing. Two adjacent sites can both be occupied, or one can be occupied and the other empty (probability $`2(E_1E_2)`$ because this can occur in two ways) or both can be empty (probability $`E_2`$). Since this covers all possibilities, the probability that two adjacent sites are occupied is therefore $`12E_1+E_2`$. The change in $`E_1`$ is thus

$$\frac{\mathrm{\Delta }E_1}{\mathrm{\Delta }t}=\frac{2D}{(\mathrm{\Delta }x)^2}(12E_1+E_2),$$
(7)
and for this equation to fit the pattern of Eq. (6) it is necessary to require the boundary condition

$$E_0=0.$$
(8)
Letting $`x=n\mathrm{\Delta }x`$ and both $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }t0`$ then leads to the continuous equation (5) with the boundary condition

$$E(0,t)=0.$$
(9)
The other boundary condition requires that the population of particles be nonvanishing:

$$E(\mathrm{},t)=0.$$
(10)
The empty interval dynamics is thus essentially the same as that of individual particles in the absence of reaction, but with different boundary conditions and a diffusion coefficient $`2D`$ that reflects the fact that the relative motion of two diffusive particles involves twice the diffusion coefficient of each particle alone.

The connection between the empty interval function and the observables of interest is obtained as follows . Back in the discrete formulation, the probability that a site is occupied (i. e. not empty) is $`1E_1`$, and the concentration of particles is therefore

$$c(t)=\frac{(1E_1)}{\mathrm{\Delta }x}c(t)=\frac{E(x,t)}{x}|_{x=0},$$
(11)
where we have exhibited the continuum limit. The interparticle distribution function $`p_n(t)`$ is the probability that the nearest neighbor to one side of a given particle is $`n`$ lattice spacings away. This probability is related to $`E_n`$ by

$$E_n=c\mathrm{\Delta }x\underset{k=n+1}{\overset{\mathrm{}}{}}\underset{m=k}{\overset{\mathrm{}}{}}p_m.$$
(12)
The first sum insures that the next neighbor is at least $`k`$ sites away and the second that $`k`$ is at least $`n+1`$. The normalization takes into account that the average distance between the particles is the reciprocal of the concentration,

$$n\mathrm{\Delta }x=\underset{n=1}{\overset{\mathrm{}}{}}np_n\mathrm{\Delta }x=\frac{1}{c}.$$
(13)
Equation (12) can be inverted:

$$cp_n=\frac{(E_{n+1}2E_n+E_{n1})}{\mathrm{\Delta }x}c(t)p(x,t)=\frac{^2E(x,t)}{x^2}$$
(14)
where again we have exhibited the continuum limit.

The construction of the kinetic equation for $`E(x,t)`$ for subdiffusive particles proceeds along arguments analogous to those used in the diffusive case. Again, one follows the motion of the particles in and out at the ends of the empty interval according to the dynamics (3). However, for particles that move subdiffusively, the concept of a hopping rate (steps per unit time) is not defined because the number of steps performed up to time $`t`$ by a subdiffusive particle goes as $`t^\alpha `$ with $`\alpha <1`$ . Therefore, the derivation of the empty interval equation must be adjusted accordingly. Proceeding as above in the discrete formulation, we see that the probability that sites $`1`$ through $`n`$ are empty but that there is a particle at site $`n+1`$ is still $`E_{n+1}E_n`$. However, the rate of decrease of $`E_n`$ due to a particle that moves from site $`n+1`$ to site $`n`$ is now described by the generalized diffusion process implicit in Eq. (3), that is, $`[K_\alpha {}_{0}{}^{}D_{t}^{1\alpha }/(\mathrm{\Delta }x)^2](E_nE_{n+1})\mathrm{\Delta }t`$. Collecting arrival and departure processes at both ends of the interval then leads to the equation

$$\frac{\mathrm{\Delta }E_n}{\mathrm{\Delta }t}=\frac{2K_\alpha {}_{0}{}^{}D_{t}^{1\alpha }}{(\mathrm{\Delta }x)^2}(E_{n+1}2E_n+E_{n1}).$$
(15)
In the continuum limit this then becomes

$$\frac{}{t}E(x,t)=2K_\alpha {}_{0}{}^{}D_{t}^{1\alpha }\frac{^2}{x^2}E(x,t).$$
(16)
The boundary conditions are exactly as before, Eqs. (9) and (10), as are the relations (11) and (14) to the observables.

### 2.2 Odd/even interval method for annihilation reaction

The empty interval method can not be applied to the annihilation reaction $`A+A0`$ because annihilation leads to a discontinuous growth of empty intervals. However, recently a new method of odd/even intervals has been introduced that leads to exact solution in the diffusion-limited case. It is based on the construction of an equation for $`r(x,t)`$, the probability that an arbitrary interval of length $`x`$ contains an even number of particles at time $`t`$ . This probability does not change if two particles inside the interval react since this process does not change the even/odd parity of the interval. Because $`r(x,t)`$ changes only by the movement of particles in or out of the ends of the interval, arguments similar to those that lead to Eq. (5) lead to the same equation for $`r(x,t)`$. This procedure can again be directly extended to the subdiffusive problem, and $`r(x,t)`$ satisfies the same fractional diffusion equation as $`E(x,t)`$:

$$\frac{}{t}r(x,t)=2K_\alpha {}_{0}{}^{}D_{t}^{1\alpha }\frac{^2}{x^2}r(x,t).$$
(17)
One boundary condition is the same as for $`E(x,t)`$ and is obtained via the same arguments used earlier:

$$r(0,t)=1.$$
(18)
The other boundary condition in general differs from that of the empty interval method and depends on the initial condition. In particular, it is determined by the parity of total initial number of particles. This parity never changes. A random initial placement leads to an equal probability that the system forever contain an even or and odd number of particles, so that

$$r(\mathrm{},t)=\frac{1}{2}.$$
(19)
The concentration of particles is related to $`r(x,t)`$ precisely as in Eq. (11):

$$c(t)=\frac{r(x,t)}{x}|_{x=0}.$$
(20)
## 3 Solution of the fractional diffusion equations and subdiffusive reaction kinetics

### 3.1 Coagulation reaction

The solution of Eq. (5) with boundary conditions (9) and (10) can be written as

$$E(x,t)=_0^{\mathrm{}}𝑑y\left[W(xy,t)W(xy,t)\right]E(y,0)+_0^{t(x/\sqrt{2K_\alpha })^{2/\alpha }}\omega _{\alpha /2}(\eta )𝑑\eta ,$$
(21)
where

$$W(x,t)=\frac{1}{\sqrt{8K_\alpha t^\alpha }}H_{11}^{10}\left[\frac{|x|}{\sqrt{2K_\alpha t^\alpha }}|\begin{array}{c}(1\frac{\alpha }{2},\frac{\alpha }{2})\hfill \\ (0,1)\hfill \end{array}\right],$$
(22)
$$\omega _\beta (x)=\frac{1}{\beta x^2}H_{11}^{10}\left[\frac{1}{x}|\begin{array}{c}(1,1)\hfill \\ (1/\beta ,1/\beta )\hfill \end{array}\right],$$
(23)
and $`H`$ is the Fox H-function . Taking into account that in Laplace transform space (indicated by a tilde over the function and the letter $`u`$ as variable)

$$\stackrel{~}{W}(x,u)=\frac{u^{\alpha /21}}{\sqrt{8K_\alpha }}\mathrm{exp}\left[\frac{u^{\alpha /2}}{\sqrt{2K_\alpha }}|x|\right]$$
(24)
and $`\stackrel{~}{\omega }_\beta (u)=\mathrm{exp}(u^\beta )`$, the solution $`E(x,t)`$ in Laplace space adopts a much simpler form:

$$\stackrel{~}{E}(x,u)=\frac{s}{2u}_0^{\mathrm{}}𝑑y\left(e^{|xy|s}e^{|x+y|s}\right)E(y,0)+\frac{1}{u}\mathrm{exp}\left(xs\right)$$
(25)
where

$$s\frac{u^{\alpha /2}}{\sqrt{2K_\alpha }}.$$
(26)
From Eqs. (11) and (25) one finds

$$\stackrel{~}{c}(u)=\frac{s}{u}\left[1s\widehat{E}(s,0)\right],$$
(27)
where $`\widehat{E}(v,t)=_0^{\mathrm{}}𝑑xe^{vx}E(x,t)`$ is the the spatial Laplace transform of $`E(x,t)`$. Equivalently, from Eq. (14) one obtains

$$\stackrel{~}{c}(u)=\frac{\lambda }{u}\left[1\widehat{p}(s,0)\right]$$
(28)
where $`\lambda c(0)`$ and $`\widehat{p}(u,0)`$ is the spatial Laplace transform of the initial interparticle distribution function $`p(x,0)`$.

A commonly considered initial interparticle distribution is the random (Poisson) distribution of average concentration $`\lambda `$, for which $`E(x,0)=e^{\lambda x}`$ and $`p_0(x)=\lambda e^{\lambda x}`$. For this initial distribution either Eq. (27) or Eq. (28) can be used to obtain

$$\stackrel{~}{c}(u)=\frac{\lambda }{u+\lambda \sqrt{2K_\alpha }u^{1\alpha /2}}$$
(29)
and $`c(t)`$ is given in closed form in terms of the Mittag-Leffler function  of parameter $`\alpha /2`$:

$$c(t)=\lambda E_{\alpha /2}\left(\lambda \sqrt{2K_\alpha }t^{\alpha /2}\right)$$
(30)
The Mittag-Leffler function can be calculated by the series expansion

$$E_\alpha (z)=\underset{k=0}{\overset{\mathrm{}}{}}\frac{z^k}{\mathrm{\Gamma }(k\alpha +1)}.$$
(31)
When $`\alpha =1`$ one recovers the usual result for diffusion-limited coagulation  since the Mittag-Leffler function of parameter $`1/2`$ is $`E_{1/2}(x)=\mathrm{exp}(x^2)\text{erfc}(x)`$.

Another initial particle distribution of interest is an ordered or periodic arrangement, $`p(x,0)=\delta (x1/\lambda )`$. In this case Eq. (28) leads to

$$\stackrel{~}{c}(u)=\frac{\lambda }{u}\left[1\mathrm{exp}\left(\frac{u^{\alpha /2}}{\lambda \sqrt{2K_\alpha }}\right)\right]$$
(32)
and $`c(t)`$ is given in terms of the Fox H-function $`H_{1,1}^{1,0}`$ :

$$\frac{c(t)}{\lambda }=1H_{11}^{10}\left[\frac{1}{\lambda \sqrt{2K_\alpha t^\alpha }}|\begin{array}{c}(1,\frac{\alpha }{2})\hfill \\ (0,1)\hfill \end{array}\right].$$
(33)
For normal ($`\alpha =1`$) diffusive particles one recovers the usual result for diffusion-limited coagulation  since $`H_{11}^{10}\left[z|\begin{array}{c}(1,\frac{1}{2})\hfill \\ (0,1)\hfill \end{array}\right]=\text{erfc}\left(\frac{z}{2}\right)`$.

It is useful to exhibit explicitly the long-time and short-time behavior of the particle concentration. The long-time behavior can be extracted via Tauberian theorems from the small-$`u`$ behavior independent of initial condition and therefore in particular common to both of the initial conditions considered above,

$$\stackrel{~}{c}(u)\frac{u^{\alpha /21}}{\sqrt{2K_\alpha }},$$
(34)
from which it follows that

$$c(t)\frac{t^{\alpha /2}}{\sqrt{2K_\alpha }\mathrm{\Gamma }\left(1\frac{\alpha }{2}\right)}.$$
(35)
For diffusive particles ($`\alpha =1`$) one recovers the well-known result $`c(t)t^{1/2}`$. The short-time behavior of the particle concentration is, not surprisingly, markedly different for the two initial conditions. For the random initial condition an expansion of Eq. (30) using (31) gives

$$\frac{c(t)}{\lambda }=1\frac{\lambda \sqrt{2K_\alpha }}{\mathrm{\Gamma }(1+\frac{\alpha }{2})}t^{\alpha /2}+\mathrm{}$$
(36)
Notice that this implies a infinite initial reaction rate because $`dc(t)/dt\mathrm{}`$ for $`t0`$, a reflection of the fact that a random distribution of particles includes a large probability of proximate particles. On the other hand, using the asymptotic expression of the Fox H-function $`H_{1,1}^{1,0}`$ for large argument  in the solution for the periodic initial distribution of particles, Eq. (33), one finds

$$\frac{c(t)}{\lambda }=1\frac{1}{\sqrt{\pi (2\alpha )}}\left(\frac{\alpha }{2}\right)^{1/(\alpha 2)}z^{1/(\alpha 2)}\mathrm{exp}\left[\frac{2\alpha }{2}\left(\frac{\alpha }{2}\right)^{\alpha /(2\alpha )}z^{2/(2\alpha )}\right]+\mathrm{}$$
(37)
where $`z=[\lambda (2K_\alpha t^\alpha )^{1/2}]^1`$. In this case $`dc(t)/dt0`$ for $`t0`$, since initially there are no proximate particles.

The curves in Fig. 1 shows the time evolution of the concentration $`c(t)`$ for the coagulation reaction and the two initial conditions considered here. The differences at early times are evident, as is the coagulation of the curves into a single asymptotic form at long times.

The time-asymptotic form of the empty interval function can be used to obtain the asymptotic behavior of the interparticle distribution function via Eq. (14). From Eq. (25) one deduces that

$$\frac{^2}{x^2}\stackrel{~}{E}(x,u)=\frac{s^3}{2u}_0^{\mathrm{}}𝑑y\left(e^{|xy|s}e^{|x+y|s}\right)E(y,0)+\frac{s^2}{u}e^{xs},$$
(38)
which for $`u0`$ reduces to

$$\frac{^2}{x^2}\stackrel{~}{E}(x,u)\frac{s^2}{u}e^{xs}.$$
(39)
From Eq. (14) and the fact that the inverse transform of (39) can be expressed in terms of a Fox H-function, one finds

$$c(t)p(x,t)\frac{1}{x^2}H_{11}^{10}\left[\frac{x}{\sqrt{2K_\alpha }t^{\alpha /2}}|\begin{array}{c}(1,\frac{\alpha }{2})\hfill \\ (2,1)\hfill \end{array}\right],$$
(40)
which in combination with Eq. (35) immediately leads to

$$p(x,t)\frac{\mathrm{\Gamma }(1\frac{\alpha }{2})}{\sqrt{2K_\alpha }t^{\alpha /2}}H_{11}^{10}\left[\frac{x}{\sqrt{2K_\alpha }t^{\alpha /2}}|\begin{array}{c}(1\alpha ,\frac{\alpha }{2})\hfill \\ (0,1)\hfill \end{array}\right]$$
(41)
for $`t\mathrm{}`$. This result is independent of the initial distribution of particles. We can also write this in the scaled form

$$p_\alpha (z)\mathrm{\Gamma }^2\left(1\frac{\alpha }{2}\right)H_{11}^{10}\left[\mathrm{\Gamma }\left(1\frac{\alpha }{2}\right)z|\begin{array}{c}(1\alpha ,\frac{\alpha }{2})\hfill \\ (0,1)\hfill \end{array}\right]$$
(42)
where $`zc(t)x`$ is the scaled interparticle distance and $`p_\alpha (z)dzp(x,t)dx`$. This stationary form is shown in Fig. 2 for several values of the diffusion exponent $`\alpha `$.

The interparticle distribution function conveys the interesting “anomalies” of the problem most clearly. For a random distribution of particles on a line this distribution is exponential. In particular, the most probable interparticle gaps are the smallest. For diffusion-limited reactions on a line it is well known that the scaled distribution deviates in two ways from the exponential behavior. First, a gap develops around each particle, and the distribution vanishes near the origin (see $`\alpha =1`$ curve in the figure), indicating an “effective repulsion” of particles. Second, the probability of large gaps decays much more rapidly than exponentially: the decay goes as a power of $`\mathrm{exp}(z^2/2)`$. In the subdiffusive case decreasing $`\alpha `$ leads to the diminution of the gap around each particle, that is, to a weakening of the effective repulsion and to a behavior that appears closer to that of a random distribution in the short-interparticle-distance behavior. This is evident in the progression of the curves with decreasing $`\alpha `$ shown in the figure. Furthermore, the probability of large gaps decays as a power of $`\mathrm{exp}(x^{2/(2\alpha )})`$, thus neatly interpolating between the purely random exponential decay state as $`\alpha 0`$ (since then $`p(x,t)c(t)e^{c(t)x}`$) and the more ordered state corresponding to diffusive particles at $`\alpha =1`$.

### 3.2 Annihilation reaction

It is straightforward to deduce that the solution $`r(x,t)`$ of Eq. (17) with initial condition $`r(x,0)=\frac{1}{2}+E(x,0)/2`$ and boundary conditions (18) and (19) is simply $`r(x,t)=\frac{1}{2}+E(x,t)/2`$ where $`E(x,t)`$ is the solution found in the previous section. We can therefore here rewrite Eq. (27) for the Laplace transform of the concentration of particle undergoing the annihilation reaction as

$$\stackrel{~}{c}(u)=\frac{s}{u}\left[1s\widehat{r}(s,0)\right]$$
(43)
where $`\widehat{r}(v,t)`$ is the the spatial Laplace transform of $`r(x,t)`$.

For a random initial distribution of mean concentration $`\lambda `$ the initial distribution is $`r(x,0)=\frac{1}{2}+\frac{1}{2}e^{2\lambda x}`$ (note that this goes to $`1/2`$ as $`x\mathrm{}`$ since there is an equal probability that we start with en even or odd number of particles, and it goes to $`0`$ as $`x1`$ since zero particles is an even number). Hence

$$\stackrel{~}{c}(u)=\frac{\lambda }{u+2\lambda \sqrt{2K_\alpha }u^{1\alpha /2}},$$
(44)
which can be inverted analytically to yield

$$c(t)=\lambda E_{\alpha /2}\left(2\lambda \sqrt{2K_\alpha }t^{\alpha /2}\right).$$
(45)
Comparing this result with Eq. (30) for the coagulation reaction we see that the rate of change of the concentration for annihilation is exactly twice that for coagulation at all times, a result that has been noted for ordinary diffusion .

For a periodic initial distribution with separation $`1/\lambda `$ between particles we have

$$r(x,0)=\{\begin{array}{cc}2j+1\lambda x\hfill & 2j/\lambda x(2j+1)/\lambda \hfill \\ 2j1+\lambda x\hfill & (2j+1)/\lambda x(2j+2)/\lambda \hfill \end{array}$$
(46)
This is a periodic function of period $`2/\lambda `$, so that

$$\widehat{r}(u,0)=\frac{_0^{2/\lambda }𝑑xe^{ux}.r(x,0)}{1e^{2u/\lambda }}$$
(47)
Carrying out this integral and substituting into Eq. (43) immediately leads to

$$\stackrel{~}{c}(u)=\frac{\lambda }{u}\mathrm{tanh}\left(\frac{s}{2\lambda }\right)=\frac{\lambda }{u}\left[\frac{1\mathrm{exp}\left(\frac{u^{\alpha /2}}{\lambda \sqrt{2K_\alpha }}\right)}{1+\mathrm{exp}\left(\frac{u^{\alpha /2}}{\lambda \sqrt{2K_\alpha }}\right)}\right]$$
(48)
where we have written the second equality to make comparison with Eq. (32) more apparent. We have not been able to invert this expression analytically to find $`c(t)`$ in closed form.

Figure 3 shows the temporal evolution of $`c(t)`$ for the two initial conditions when $`\alpha =1/2`$. The evolution for a periodic initial distribution has been obtained by numerically inverting Eq. (48).

The long-time behavior is again seen to be independent of initial condition and can be obtained directly from Eq. (43). Since $`\widehat{r}(s,0)(2s)^1`$ when $`s0`$, it follows that as $`u0`$,

$$\stackrel{~}{c}(u)\frac{u^{\alpha /21}}{2\sqrt{2K_\alpha }}$$
(49)
which in turn implies that

$$c(t)\frac{t^{\alpha /2}}{2\sqrt{2K_\alpha }\mathrm{\Gamma }\left(1\frac{\alpha }{2}\right)},$$
(50)
when $`t\mathrm{}`$. Comparison with Eq. (35) shows that at long times the concentration for the annihilation reaction is exactly half that of the coagulation reaction. This is in agreement with the general result obtained above comparing the rates of change of the concentrations. The short-time behaviors for the two initial conditions are different, as already evident in Fig. 3. For the random initial condition an expansion of Eq. (45) leads to

$$\frac{c(t)}{\lambda }=1\frac{2\lambda \sqrt{2K_\alpha }}{\mathrm{\Gamma }(1+\frac{\alpha }{2})}t^{\alpha /2}+\mathrm{}$$
(51)
\[compare with Eq. (36)\]. For the periodic initial condition expansion of the denominator of Eq. (48) for large $`u`$ and comparison with Eq. (32) leads to

$$\frac{c(t)}{\lambda }=1\frac{2}{\sqrt{\pi (2\alpha )}}\left(\frac{\alpha }{2}\right)^{1/(\alpha 2)}z^{1/(\alpha 2)}\mathrm{exp}\left[\frac{2\alpha }{2}\left(\frac{\alpha }{2}\right)^{\alpha /(2\alpha )}z^{2/(2\alpha )}\right]+\mathrm{}$$
(52)
where again $`z=[\lambda (2K_\alpha t^\alpha )^{1/2}]^1`$ \[compare with Eq. (37)\].

## 4 Reactions in a lattice

It is sometimes convenient not to carry out the continuous limit in the formulation of the reaction kinetics for coagulation or annihilation reactions. Not only may the actual physical system be discrete, but simulations usually involve discrete lattices, and finite reactant size effects (i.e., small distance scale effects) are more appropriately dealt with through discrete formulations . Indeed, a discrete viewpoint was the starting point of the continuum equations considered in the previous section. Furthermore, the discrete formulation allows consideration of different reaction rules such as $`A\mathrm{}AAAA`$ .

The hierarchy of differential-difference equations on a lattice of $`L`$ sites is

$`I_0(t)`$ $`=`$ $`1`$
$`{\displaystyle \frac{dI_n}{dt}}`$ $`=`$ $`a{}_{0}{}^{}D_{t}^{1\alpha }\left[I_{n+1}(t)2I_n(t)+I_{n1}(t)\right]1nL1`$ (53)
for $`n=1,2,\mathrm{}`$, where $`a2K_\alpha /(\mathrm{\Delta }x)^2`$, $`I_n(t)E(x_n=n\mathrm{\Delta }x,t)`$ for the coagulation reaction, and $`I_n(t)r(x_n=n\mathrm{\Delta }x,t)`$ for the annihilation reaction. The other boundary condition, at $`n=L`$, depends on which reaction is under consideration and will be stated below. The discrete form

$$c(t)=\frac{1I_1(t)}{\mathrm{\Delta }x}$$
(54)
was already introduced in Eq. (11). From here on, we take $`\mathrm{\Delta }x=1`$.

### 4.1 Coagulation in a segment

The problem of the coagulation dynamics in a discrete segment has been solved for the normal diffusive case . Here we present the solution for the coagulation process of subdiffusive particles on a chain of $`L`$ sites with periodic boundary conditions, which is described by the hierarchy (53) together with the additional boundary condition

$$I_L(t)=0$$
(55)
We solve this set of equations by means of the ansatz

$$I_n(t)=\underset{\omega }{}b_n(\omega )E_\alpha (a\omega t^\alpha )$$
(56)
which is closely related to the method of separation of variables for the subdiffusive fractional differential equation . Again, $`E_\alpha `$ is the Mittag-Leffler function of parameter $`\alpha `$. Taking into account that

$$\frac{d}{dt}E_\alpha (\beta t^\alpha )=\beta {}_{0}{}^{}D_{t}^{1\alpha }E_\alpha (\beta t^\alpha ),$$
(57)
the resulting eigenvalue problem can be solved in the same way as for normal diffusion . Notice that $`E_\alpha (x)=e^x`$ for $`\alpha =1`$ so that the ansatz (56) takes the well-known form for normal diffusive particles (for which $`\alpha =1`$). The full solution is

$$I_n(t)=\left(1\frac{n}{L}\right)+\underset{m=0}{\overset{L1}{}}A_m\mathrm{sin}\left(m\frac{\pi }{L}(nL)\right)E_\alpha (a\omega _mt)$$
(58)
with

$$\omega _m=2[1\mathrm{cos}(m\pi /L)].$$
(59)
The coefficients $`A_m`$ are determined from the initial conditions.

For an initially random distribution of particles of concentration $`\lambda `$ we have $`I_n(0)=(1\lambda )^n`$ for $`n=0,1,\mathrm{},L1`$, so that

$$A_m(t)=\frac{2}{m\pi }\frac{(1\lambda )^Lm^2\pi ^2+(1)^mL^2\mathrm{ln}^2(1\lambda )}{m^2\pi ^2+L^2\mathrm{ln}^2(1\lambda )}$$
(60)
For an initially periodic distribution of particles separated by $`1/\lambda `$ lattice sites we have $`I_n(0)=1\lambda n`$ for $`\lambda n1`$ and $`I_n(0)=0`$ for $`\lambda n1`$. In this case

$$A_m(t)=(1)^{m+1}\frac{2\lambda L\mathrm{sin}\left(m\pi /\lambda L\right)}{m^2\pi ^2}.$$
(61)
Figure 1 shows the decay of the concentration $`c(t)`$ given by Eq. (54) with (58) when the initial distribution of subdiffusive particles on the lattice is random and when it is periodic. At longer times both solutions approach those of the continuum equations and would do so more clearly if we were to exhibit the solutions to our difference equations for larger lattices. The marked deviation of the lattice solution from the continuous asymptotic $`t^{\alpha /2}`$ decay occurs when the concentration on the lattice approaches its asymptotic (minimum) value of one particle on the entire lattice, $`c(t\mathrm{})=1/L`$. This time can be estimated using the continuum solution to occur at time $`t_\times `$ such that $`c(t_\times )=1/L`$. In the figure this separation of solutions clearly does occur when $`c/\lambda =1/\lambda L=0.1`$, i.e., when $`\mathrm{log}_{10}(c/\lambda )=1.0`$.

While the long-time behavior of the discrete and continuous solutions coincide (independently of initial condition) until the reactant concentration is extremely low, the short-time behavior of the two solutions is expected to differ, especially when the initial distribution of particles is random. The difference has been established for particles undergoing normal diffusion . To determine the short-time behavior on the lattice we assume that the lattice is initially full, so that $`I_n(0)=0`$ for $`n1`$. (The reasoning follows along the same lines for an arbitrary concentration of particles, but the equations are more cumbersome.) Taking the time Laplace transform of Eq. (53) gives

$$\stackrel{~}{I}_{n+1}(u)(2+\frac{u^\alpha }{a})\stackrel{~}{I}_n(u)+\stackrel{~}{I}_{n1}(u)=0,n1.$$
(62)
This set of difference equations must be solved with the boundary conditions $`\stackrel{~}{I}_0(u)=1/u`$ and $`\stackrel{~}{I}_L(u)=0`$. Since we are interested in the short-time (large-$`u`$) behavior, we can set $`L\mathrm{}`$ without appreciably affecting the outcome. With these boundary conditions we essentially follow the procedure of Ref.  to find

$$\stackrel{~}{I}_n(u)=\frac{1}{2^nu}\left[2+\frac{u^\alpha }{a}\left(\frac{u^{2\alpha }}{a^2}+4\frac{u^\alpha }{a}\right)^{1/2}\right]^n.$$
(63)
In particular, $`\stackrel{~}{I}_1(u)=a/u^{1+\alpha }+\mathrm{}`$ for $`u\mathrm{}`$, so that with the help of appropriate Tauberian theorems we find that $`I_1(t)=at^\alpha /\mathrm{\Gamma }(1+\alpha )+\mathrm{}`$ for $`t0`$. Using Eq. (54), we calculate the short-time behavior of the reactant concentration

$$\frac{c(t)}{\lambda }=1\frac{2K_\alpha \lambda ^2}{\mathrm{\Gamma }(1+\alpha )}t^\alpha +\mathrm{}.$$
(64)
Note that with an initially full lattice there is no distinction between an initially random and initially periodic distribution of particles. For an arbitrary random initial concentration one still finds that $`c(t)/\lambda 1t^\alpha `$. Comparison with Eq. (36) shows that at short times, as in the case of normal diffusion , the concentration on a lattice decays more slowly than in the continuum. This is not easily visible on the scale of Fig. 1.

### 4.2 Annihilation in a segment

Now consider the annihilation reaction $`A+A0`$ and suppose that initially there are $`N`$ particles on a segment of length $`L`$. The hierarchy (53) is now augmented with the second boundary condition

$$I_L(t)=\{\begin{array}{cc}1,\hfill & N=\text{even}\hfill \\ 0,\hfill & N=\text{odd}\hfill \end{array}$$
(65)
where we note that now $`I_n(t)r(x_n=n\mathrm{\Delta }x,t)`$. This condition is of course a result of the fact that the reaction does not change the parity of the number of surviving particles. We proceed as in the coagulation problem and write the solution as

$$I_n(t)=I_n^{}+\underset{m=0}{\overset{L1}{}}A_m\mathrm{sin}\left(m\frac{\pi }{L}(nL)\right)E_\alpha (a\omega _mt)$$
(66)
where

$$I_n^{}=I_n(\mathrm{})=\{\begin{array}{cc}1,\hfill & N=\text{even},\hfill \\ 1\frac{n}{L},\hfill & N=\text{odd},\hfill \end{array}$$
(67)
$`\omega _m`$ is given in Eq. (59), and the $`A_m`$ are coefficients determined from the initial conditions.

For a random initial condition of particle concentration $`\lambda `$ we have $`I_n(0)=\frac{1}{2}+\frac{1}{2}(12\lambda )^n`$ for $`n=0,1,\mathrm{},L1`$, so that

$$A_m(t)=\frac{1}{m\pi }\frac{\left[(1)^N(12\lambda )^L\right]m^2\pi ^2+\left[1(1)^{m+N}\right]L^2\mathrm{ln}^2(12\lambda )}{m^2\pi ^2+L^2\mathrm{ln}^2(12\lambda )}.$$
(68)
For a periodic initial distribution with particle separation $`1/\lambda `$ we have

$$I_n(0)=\{\begin{array}{cc}2j+1\lambda n\hfill & 2j/\lambda n(2j+1)/\lambda \hfill \\ 2j1+\lambda n\hfill & (2j+1)/\lambda n(2j+2)/\lambda \hfill \end{array}$$
(69)
which leads to

$$A_m(t)=(1)^N\left[1(1)^N(1)^m\right]^2\frac{\lambda L}{m^2\pi ^2}\mathrm{tan}\left(\frac{m\pi }{2\lambda L}\right).$$
(70)
Figure 3 shows the evolution of the concentration as given by Eq. (54) with (66) for a random initial distribution (circles) and a periodic initial distribution (squares). A number of points are noteworthy. First, the solution of course decays to a finite constant (one particle remaining) when there is initially an odd number of particles in the system, and to zero when the initial number of particles is even. The deviation of the solutions on a finite segment from those in an infinite continuum again set in around the time when $`c=1/L`$ (strictly speaking around $`2/L`$ for the even case). Note that in the continuum solution one does not distinguish between the even and odd particle number cases: the boundary condition (19) is an average of the two, and so the continuum curves should be compared to the average of the discrete ones for the two boundary conditions.

## 5 Conclusions

We have considered the coagulation dynamics $`A+AA`$ and the annihilation dynamics $`A+A0`$ for particles moving subdiffusively in one dimension. This scenario combines the “anomalous kinetics” and “anomalous diffusion” problems, each of which leads to interesting dynamics separately and to even more interesting dynamics in combination. The fractional diffusion equation plays a central role in our analysis and allows the exact calculation of the density $`c(t)`$ within this formulation. We have calculated $`c(t)`$ explicitly for all times for the coagulation reaction in a one-dimensional continuum for two initial distributions, a random (Poisson) distribution and a periodic distribution. This calculation is based on the empty interval method. Using the odd/even interval formulation, we have also obtained an explicit solution for the annihilation reaction with a random initial distribution. For a periodic initial distribution we are only able to calculate the Laplace transform of $`c(t)`$ analytically and must perform the inversion numerically.

For the coagulation reaction the empty interval method also leads to an explicit solution for the interparticle distribution function $`p(x,t)`$. At long times we find a universal expression for this distribution (i.e., one independent of initial distribution), in terms of a single scaled variable. Anomalous diffusion is characterized by the exponent $`\alpha `$ introduced in Eq. (1), ordinary diffusion corresponding to $`\alpha =1`$. Deviations from ordinary diffusion lead to a curious interplay. On the one hand, with decreasing $`\alpha `$ (and hence increasingly subdiffusive motion) the decay of the particle density towards extinction becomes increasingly slower and in this sense increasingly different from law of mass action behavior. On the other hand, for the coagulation reaction the spatial distribution of initially randomly distributed reactants remains more Poissonian for all time as $`\alpha `$ decreases; indeed, as $`\alpha `$ deviates from unity the relatively empty regions around each particle that are tantamount to (and indeed explain) anomalous kinetics in the usual diffusion-limited case become more populated.

We have also considered the problem of subdiffusion-limited reaction on discrete lattices. We have presented a hierarchy of differential-difference equations for the empty interval probability $`I_n(t)E(x_n,t)`$ and the even/odd interval probability $`I_n(t)r(x_n,t)`$, and have expressed the solutions $`I_n(t)`$ as a superposition of subdiffusive modes whose decay is governed by the Mittag-Leffler function in the same way that the exponential function governs the decay of ordinary diffusive modes. As in ordinary diffusion , the associated concentrations $`c(t)`$ differ from the continuum solutions at early times. At long times the lattice solutions deviate from those of the continuum systems when there is only a small number (one or two) of particles remaining in the system.

## 6 Acknowledgments

This work has been supported in part by the Ministerio de Ciencia y Tecnología (Spain) through Grant No. BFM2001-0718 and by the Engineering Research Program of the Office of Basic Energy Sciences at the U. S. Department of Energy under Grant No. DE-FG03-86ER13606. SBY is also grateful to the DGES (Spain) for a sabbatical grant (No. PR2000-0116) and to the Department of Chemistry and Biochemistry of the University of California, San Diego for its hospitality.


format: text

text:
# A generalization of Abel Inversion to non axisymmetric density distribution

## Acknowledgements

On of the authors (P.T.) wish to acknowledge support from the italian M.U.R.S.T. (Project: ”Metodologie e diagnostiche per materiali e ambiente”). Authors are very grateful to D. Giulietti, from the Dep. of Physics, Univ. of Pisa and to L.A. Gizzi and R. Numico from IFAM-CNR, Pisa, for useful discussions and their encouragement.

## Figures Caption

## .


text:
# High-Mass Stars as Early Signpost of Cluster Formation

## 1. Introduction

High-mass stars, although rare, are a major contributor of energy to the ambient ISM, be it during their birth in the form of jets and outflows, during their UV-radiating lifetime, or during their spectacular death as supernovae. They are commonly detected in rich clusters of stars and to date an example of isolated high-mass star formation has not been observed. Thus, if one wants to understand the process of high-mass star formation, one must understand the physical processes which lead to the formation of entire clusters of stars. Furthermore, the dense clusters may influence the formation of massive stars, and may even be a required precursor for high-mass star formation. Is the formation of high-mass stars a natural implication of the stellar density and the possible interaction and coalescence of low-mass stars in large clusters? Do low-mass stars therefore need to form before a high-mass star can be fused? Or is it that high-mass stars simply form independently in the same dense environment as low-mass stars through the same processes of mass accretion?

In this contribution, we present observations which show the presence of massive reservoirs of seemingly quiescent dense gas in close proximity to highly active cluster-forming cores. We argue that these reservoirs are at an early stage of star formation, where the usual signposts of massive stars and cluster formation – strong mid-IR emission, bright radio continuum sources, and dense clusters of stars in 2 $`\mu `$m imaging – are absent. Instead, the only signatures of star formation comes from deeply embedded high-mass stars, namely masers, powerful outflows, and luminous submm sources. We propose that some of the youngest regions of cluster formation may be found by using these early signposts of high-mass star formation, and describe an ongoing search for such regions.

## 2. High-mass Star Formation in Clusters

From observations of 100 open clusters, Battinelli & Capuzzo-Dolcetta (1991) concluded that the star formation rate in clusters in our galaxy is 225 MMyr-1 kpc-2, whereas the SFR near the sun is 3000–5000 MMyr-1kpc-2. From observations of Perseus, L1641, 50% of L1630, and Mon R2 as part of the Two Micron All Sky Survey, Carpenter (2000) deduces that the largest fraction of stars must form in clusters. In fact, 70% of the young stars identified in Carpenter’s analysis are in clusters with more than 100 members. This then leads one to the conclusion that the 10% still seen as open clusters are simply those that remain bound long after their parental molecular clouds are dissipated. The dispersal of the molecular gas in clusters is very much a function of their most massive members. High-mass stars ($`>5`$ M, corresponding to spectral type B6 or earlier) are always found to be associated with rich clusters with stellar densities of more than 1000 stars per pc3. These clusters are thus affected by the dispersive jets, the intense UV fields and, if the cluster remains bound after the dispersal of the molecular cloud, the ensuing supernovae of the massive stars.

From observation of Herbig Ae/Be stars, Testi, Palla, & Natta (1999) conclude that the lack of clusters around stars later than B7–B8 is an imprint of the stellar formation mode for high-mass stars. The process of high-mass star formation within dense clusters through spherical accretion has been a long-standing problem. At around 8–10 M the radiation pressure of the forming star begins to inhibit further infall of the material accreting (e.g. Yorke & Kruegel 1977, Kahn 1974, Larson & Starrfield 1971). Models including rotating protostellar grainy disks (e.g. Jijina & Adams 1996, Wolfire & Cassinelli 1987), magnetic fields (Balsara et al. 2001, Nakano 1998), and accretion rates growing with time (Behrend & Maeder 2001) can certainly push this limit very much higher, but the theories about the physical backgrounds for the accretion process are still open for much discussion. At stellar densities of around 107pc-3 accretion induced collision may begin to have a significant effect on the formation of high-mass stars up to the point where fusion of lower mass stars may actually form the more massive members (e.g. Bonnell et al. 1998, Stahler et al. 1998). These two scenarios for forming massive stars in clusters – accretion and star coalescence – are the focus of the current debate.

## 3. Early Signposts of High-Mass Stars

High-mass stars are typically formed in warm ($`30`$ K), dense (n(H2) $`10^5`$ cm-3), compact ($`50,000`$ AU) and massive ($`1000`$ M) dust cores, easily traced in CS, CH3CN, HC3N, or NH3 and the submm-continuum with typical FIR luminosities of $`10^310^5`$ L. These cores are usually associated with H2O, CH3OH, H2CO, NH3, and/or OH masers and often, at a later stage when the most massive members are emerging from their natal cocoons, hyper- or ultracompact H ii regions. Often one can find shock-excited H2 emission knots at 2.1 $`\mu `$m and/or high-velocity bipolar outflows from broad wings in SiO or 13CO molecular emission lines. These molecular outflows typically have masses of around 150 (10–5000) M and outflow rates of $`3\times 10^53\times 10^2`$ Myr-1 (Garay & Lizano 1999). Outflows and certain types of masers appear very early on in the formation process, long before luminous IR sources or H ii regions can be detected. Sheperd & Churchwell (1996) have shown that the outflow properties seem to be correlated with the luminosity of the embedded sources driving the outflows. Although new observation by Henning et al. (2000) are in support of this, a recent analysis of a number of massive outflows seems to indicate that this empirical relationship breaks down with increasingly higher mass (cf. contribution by Beuther et al., this issue).

### 3.1. W3 SMS3

W3 Main is a well-studied and very active region of high-mass star formation in the W3 GMC at a distance of 2.3 kpc. The W3 Main region is 10 pc north of the W3(OH) region and consists of two main molecular cores: W3 East & West. Towards W3 East we find a dense stellar cluster associated with the luminous infrared source IRS 5 embedded in the compact molecular core W3 SMS1 (Tieftrunk et al. 1998). Towards this region Megeath et al. (1996) detected the highest density of stars: a cluster of $`80300`$ low-mass stars, corresponding to a stellar density of a few $`10^3`$ pc-3. The bright source IRS 5 itself is associated with a massive bipolar outflow (Hasegawa et al. 1994) and a cluster of H2O masers where Claussen et al. (1994) have detected a number of hypercompact H ii regions with diameters less than 1000 AU. Towards IRS 4, embedded in the compact molecular core W3 SMS2, Tieftrunk et al. (1997) found another similar hypercompact H ii region. The presence of hypercompact sources indicates that the massive stars are just beginning to ionize their surrounding cores (Tieftrunk et al. 1997). Just to the south of IRS 4, is the molecular core W3 West/SMS3 (Ladd et al. 1993). This dense core is seemingly quiescent: it is devoid of mid-IR sources and radio continuum sources, and 2 $`\mu `$m imaging has failed to detect a cluster in SMS3 (Tieftrunk et al. 1998). Since SMS3 has a molecular mass and density higher than the surrounding cluster-forming cores, Tieftrunk, Gaume & Wilson (1998) argued that the SMS3 core may be at an early stage of cluster formation.

### 3.2. NGC 6334 I(N)

Towards the north-eastern edge of the giant molecular cloud NGC 6334 (Fig. 2), at a distance of 1.7 kpc (Neckel 1978), another massive twin core system can be found: NGC 6334 I and I(N) (Fig. 3). As W3 SMS1 & 3, these two cores present a unique duet for studying the evolution of stellar clusters and the interstellar chemistry in young star-forming cores.

NGC 6334 I, with a total mass of $`1000`$ M and a dust temperature of $`100`$ K, is a dense and compact molecular hot core associated with an ultracompact cometary H ii as well as a strong IR continuum region and a group of CH3OH, H2O, OH, and NH3 masers. It is considered to be the youngest of the star-forming regions in the NGC 6334 GMC, exhibiting the richest chemistry of all sources. In fact, our SEST data shows that this source is comparable to Orion-KL and Sgr B2 in the richness of its spectrum. However, in comparison to these sources where the preponderance of line blending is a principal problem for assigning emission lines or detecting new transitions, line widths in NGC 6334 I are 2–4 times more narrow allowing for a significantly higher detection rate of weaker features (Fig. 4). In fact, we detected a plethora of line features down to a noise level of 0.01 K and expect that Sgr B2 and Orion-KL will be replaced by NGC 6334 I as the new standard source in searches for new interstellar molecules and transitions (Tieftrunk & Thorwirth, in prep.).

Only 1.5 north (1 pc at 1.7 kpc) of NGC 6334 I a second core of roughly equal size, comparable mass, and a molecular hydrogen density of $`10^6`$ cm-3 (Kuiper et al. 1995) is apparent (Sandell 2000). Designated NGC 6334 I/(North) \[hereafter I(N)\] it is notable for the lack of mid-IR emission: to date, no thermal infrared emission has been detected from this region at wavelengths $`<100\mu `$m (Gezari 1982). It does not contain an H ii region. However, H2O and CH3OH (class I and II) masers have been detected in the core. Very recently, Megeath & Tieftrunk (1999) presented evidence from SEST and ESO-2.2m data for multiple outflows originating in I(N) (Fig. 5). Using the SiO and 13CO emission in the line wings towards I(N), we deduced an outflow mass of 7 M, transporting $`10^3`$ M yr-1. From the well-known relationship of Shepherd & Churchwell (1996) we thus estimated a luminosity of $`10^3`$ L and a mass of 10 M for the driving source. Thus, the I(N) core already contains at least one high-mass (proto)star and is possibly on the verge of cluster formation (cf. also Sandell 2000). From these observations we conclude that the first detectable evidence for cluster formation are the signposts of outflows and masers of the associated high-mass stars long before a dense 2 $`\mu `$m cluster itself appears. Thus the formation of intermediate to high-mass stars can at least precede the appearance of the associated cluster of low-mass stars. This, however, does not imply that low-mass stars would actually form later as well; they may be too deeply embedded to become visible to the observer at an early stage of the cluster formation.

## 4. IRAS-quiet maser clusters

Surveys of the earliest stages of massive star formation have concentrated on UC H ii regions (e.g. Brand et al. 1984, Wood & Churchwell 1989), some associated with IRAS sources (e.g. Bronfman et al. 1996, Osterloh et al. 1997), and more recently, on luminous IRAS sources which are not associated with H ii regions (e.g. Zhang et al. 2001, Molinari et al. 2000, Sridharan et al. 1999). However, as demonstrated by NGC 6334 I(N) and W3 SMS3, the earliest phases of massive star formation are probably both radio and infrared-quiet – analogous to the Class 0 phase of low-mass star formation – and searches using the IRAS database may miss some of the youngest, most interesting sources. Due to the lack of unbiased submillimeter surveys of the galaxy, only a few truly deeply embedded massive protostars or Class 0 objects are known.

From the ATCA survey by Walsh et al. (1998) we had selected 46 maser clusters that were radio- and IR-quiet, i.e. not associated with UC H ii regions and more than 1 pc (1–14) away from the nearest IRAS point source. With the SEST we detected broad SiO emission line wings towards 15 of these, indicative of outflows. It is clear that these CH3OH maser clusters with SiO outflows must have a driving source. These sources, however, are inconspicuous in the mid-IR and thus must be deeply embedded. NGC 6334 I(N), which has a comparable flux density to NGC 6334 I at 400 $`\mu `$m, is much brighter at 1 mm (Gezari 1982). Thus, these deeply embedded objects appear to be the prominent sources when observed at 1 mm in a system of star-forming cores. With SIMBA, the 37-multichannel bolometer of SEST, we will now observe the 1200 $`\mu `$m flux of the maser sources. Since dust emission at (sub)mm wavelengths is generally optically thin, it provides a much better contrast between high and low density gas than most molecular lines. From this sample of regions, we expect to find not only the youngest examples of high-mass star formation, but some of the youngest regions of cluster formation. With such a sample we can better study the relationship between high- and low-mass stars, as well as construct a better picture of the physical and chemical environment of molecular clouds at the onset of cluster formation.

## References

Balsara, D. S., Ward-Thompson, D., & Crutcher, R. M. 2001, MNRAS, 327, 715

Battinelli, P. & Capuzzo-Dolcetta, R. 1991, MNRAS, 249, 76

Behrend, R. & Maeder, A. 2001, A&A, 373, 190

Bonnell, I. A., Bate, M. R., & Zinnecker, H. 1998, MNRAS, 298, 93

Brand, J., van der Bij, M. D. P., de Vries, C. P., Leene, A., Habing, H. J., Israel, F. P., de Graauw, T., van de Stadt, H. & Wouterloot, J. G. A. 1984, A&A, 139, 181

Bronfman, L., Nyman, L.-A., & May, J. 1996, A&AS, 115, 81

Carpenter, J. M. 2000, AJ, 120, 3139

Claussen, M. J., Gaume, R. A., Johnston, K. J., & Wilson, T. L. 1994, ApJ, 424, L41

Garay, G. & Lizano, S. 1999, PASP, 111, 1049

Gezari, D. Y. 1982, ApJ, 259, L29

Hasegawa, T. I., Mitchell, G. F., Matthews, H. E., & Tacconi, L. 1994, ApJ, 426, 215

Henning, Th., Schreyer, K., Launhardt, R., & Burkert, A. 2000, A&A, 353, 211

Jijina, J. & Adams, F. C. 1996, ApJ, 462, 874

Kahn, F. D. 1974, A&A, 37, 149

Kuiper, T. B. H., Peters, W. L., III, Foster, J. R., Gardner, F. F., & Whiteoak, J. B. 1995, ApJ, 446, 692

Ladd, E. F., Deane, J. R., Sanders, D. B., Wynn-Williams, C. G. 1993, ApJ, 419, L186

Larson, R. B. & Starrfield, S. 1971, A&A, 13, L190

Megeath, S. T. & Tieftrunk, A. R. 1999, ApJ, 526, L113

Megeath, S. T., Herter, T., Beichman, C., Gautier, N., Hester, J. J., Rayner, J., & Shupe, D. 1996, A&A, 307, 775

Molinari, S., Brand, J., Cesaroni, R., & Palla, F. 2000, A&A, 355, 617

Nakano, T. 1998, ApJ, 494, 587

Neckel, T. 1978, A&A, 69, 51

Osterloh, M., Henning, Th., & Launhardt, R. 1997, ApJS, 110, 71

Sandell, G. 2000, A&A, 358, 242

Shepherd, D. S. & Churchwell, E. 1996, ApJ, 472, 225

Stahler, S. W., Palla, F., & Ho, P. T. 1998, in: Protstars and Planets IV, eds. Mannings, V. et al. Univ. of Arizona Press, Tucson, 397

Sridharan, T. K., Menten, K. M., Wyrowski, F., & Schilke, P. 1999, in: Proceedings of Star Formation 1999, eds. Nakamoto, T., Nobeyama Radio Observatory, 183

Testi, L., Palla, F., & Natta, A. 1999, A&A, 342, 515

Tieftrunk, A. R., Gaume, R. A., Claussen, M. J., Wilson, T. L., & Johnston, K. J. 1997, A&A, 318, 931

Tieftrunk, A. R., Gaume, R. A., & Wilson, T. L. 1998, A&A, 340, 232

Tieftrunk, A. R., Megeath, S. T., Wilson, T. L., & Rayner, J. T. 1998, A&A, 336, 991

Walsh, A. J., Burton, M. G., Hyland, A. R., & Robinson, G. 1998, MNRAS, 301, 640

Wolfire, M. G. & Cassinelli, J. P. 1987, ApJ, 319, 85

Wood, D. O. S. & Churchwell, E. 1989, ApJS, 69, 831

Yorke, H. W. & Kruegel, E. 1977, A&A, 54, 183

Zhang, Q., Hunter, T. R., Brand, J., Sridharan, T. K., Molinari, S., Kramer, M. A., & Cesaroni, R. 2001, ApJ, 552, 167


text:
# Heavy Flavor in Photoproduction at HERA

## 1 Introduction

During the first phase of operation (1992–2000), HERA collided electrons and positrons with energy $`26.727.6`$GeV and protons with energy $`820920`$GeV yielding a center-of-mass energy of $`296318`$GeV. In this period, extensive measurements of heavy flavor photoproduction were made by the H1 and ZEUS collaborations . Some recent charm and beauty photoproduction results are discussed in this paper.

In photoproduction processes at HERA, a quasi-real photon with virtuality $`Q^20`$ is emitted by the incoming electron and interacts with the proton. At leading order (LO) in QCD, two types of processes are responsible for the production of heavy quarks: the direct photon processes, where the photon participates as a point-like particle, and the resolved photon processes, where the photon acts as a source of partons. The dominant direct photon process is photon-gluon fusion where the photon fuses with a gluon from the incoming proton. In resolved photon processes, a parton from the photon scatters off a parton from the proton. Charm and beauty quarks present in the parton distributions of the photon, as well as of the proton, lead to processes like $`cgcg`$ and $`bgbg`$, which are called heavy flavor excitation processes. In next-to-leading order (NLO) QCD, only the sum of direct and resolved processes is unambiguously defined.

## 2 Charm photoproduction

Differential cross sections for $`D^\pm `$ and $`D_s^\pm `$ photoproduction in $`p_{}^D`$ and $`\eta ^D`$ were measured in the kinematic range $`Q^2<1`$ GeV2, $`130<W<280`$ GeV, $`3<p_{}^D<12`$ GeV and $`|\eta ^D|<1.5`$ . Here $`W`$ is the $`\gamma p`$ center-of-mass energy, and $`p_{}^D`$ and $`\eta ^D`$ are the $`D`$-meson transverse momentum and pseudorapidity, respectively. The pseudorapidity $`\eta `$ is defined as $`\mathrm{ln}(\mathrm{tan}\frac{\theta }{2})`$, where the polar angle $`\theta `$ is measured with respect to the proton beam direction.

In Fig. 1, NLO calculations  obtained with the MRS(G) and GRV-G HO  parton density parametrisations for the proton and photon, respectively, are compared with the differential cross sections. The thick curves were obtained with the renormalization scale $`\mu _R=m_{}\sqrt{m_c^2+p_{}^2}`$ ($`m_c=1.5`$ GeV) and the factorization scales of the photon and proton structure functions were set to $`\mu _F=2m_{}`$. For the thin curves, a rather extreme value for the pole $`c`$-quark mass, $`m_c=1.2`$ GeV, and a $`\mu _R`$ value of $`0.5m_{}`$ were used. The Peterson fragmentation function  was used for charm fragmentation in this calculation. Following the results of the NLO fits  to ARGUS data , the same values of the Peterson parameter, $`ϵ=0.035`$, were used for both $`D^\pm `$ and $`D_s^\pm `$ cross section calculations. The fractions of $`c`$ quarks hadronizing as $`D^\pm `$ or $`D_s^\pm `$ mesons were extracted from results on charm production in $`e^+e^{}`$ annihilations . The NLO calculations underestimate the measured cross sections. The shapes of the $`p_{}^D`$ distributions are not completely reproduced. For the $`\eta ^D`$ distributions, the NLO predictions are below the data in the central and forward (proton direction) regions.

An experimental separation of the direct and resolved processes was obtained by using the variable $`x_\gamma ^{\mathrm{obs}}`$, which is the fraction of the photon momentum contributing to the production of the two jets with the highest transverse energies within the accepted pseudorapidity range. The charm photoproduction differential cross section was measured as a function of $`x_\gamma ^{\mathrm{obs}}`$ . A comparison of the $`x_\gamma ^{\mathrm{obs}}`$ distribution with LO Monte Carlo (MC) simulations indicated the existence of charm excitation in the photon parton density. Recently, the ZEUS collaboration has performed new measurements of charm dijet photoproduction . The angle between the jet-jet axis and the beam axis in the dijet rest frame has been approximated by the variable $`\mathrm{cos}\theta ^{}`$, which is a function of the pseudorapidities of the two jets:

$$\mathrm{cos}\theta ^{}=\mathrm{tanh}(\frac{\eta ^{jet1}\eta ^{jet2}}{2}).$$

In the previous inclusive dijet analysis , it was shown that in direct photon processes, where the propagator in the leading-order (LO) QCD diagrams is a quark, the differential cross section $`d\sigma /d|\mathrm{cos}\theta ^{}|`$ rises slowly towards values of $`|\mathrm{cos}\theta ^{}|1`$, while in resolved photon processes, where in most cases the propagator is a gluon, it rises steeply with increasing $`|\mathrm{cos}\theta ^{}|`$, consistent with the behavior of a spin-1 propagator.

The differential distribution $`dN/d|\mathrm{cos}\theta ^{}|`$ for events with at least two jets and a $`D^\pm `$ is shown in Fig. 2. The data points are given separately for direct photon (open dots) and for resolved photon (black dots) events. In this analysis, a direct (resolved) photon process was defined by the selection $`x_\gamma ^{\mathrm{obs}}>0.75`$ ($`x_\gamma ^{\mathrm{obs}}<0.75`$). The dashed (full) histogram is the PYTHIA  distribution for the direct (resolved) photon MC events. All the distributions are normalized to the resolved data distribution in the first 4 bins. An enhancement at high values of $`|\mathrm{cos}\theta ^{}|`$ is seen for the resolved photon sample, both in the data and in the MC simulation. The direct photon samples do not show this strong peaking. This observation is consistent with a significant gluon-exchange contribution and, consequently, with a significant contribution of charm excitation processes to charm photoproduction at HERA energies.

## 3 Beauty photoproduction

Beauty photoproduction cross sections were measured by the H1 and ZEUS collaborations using events with at least two jets and a lepton in the final state . The beauty signals were extracted by fitting the $`p_T^{\mathrm{rel}}`$ distributions of the data with Monte Carlo predictions for beauty and the lighter flavor components, where $`p_T^{\mathrm{rel}}`$ is the transverse momentum of the lepton relative to the axis of the associated jet. The cross sections were found to be above NLO QCD expectations.

Recently, the H1 collaboration has performed new measurements of the beauty photoproduction using a central silicon vertex detector . The cross section was extracted from the $`p_T^{\mathrm{rel}}`$ and impact parameter ($`\delta `$) distributions of muons in dijet events. Fig. 3 shows the observed impact parameter distribution in the data together with histograms indicating the contributions from beauty production and from backgrounds. The decomposition was obtained from a likelihood fit using the shapes of the $`\delta `$ distributions of beauty and charm events from Monte Carlo simulations and of fake muons from real data. Using both $`p_T^{\mathrm{rel}}`$ and $`\delta `$ observables, and in combination with the earlier result, the open beauty cross section was determined in the visible range $`Q^2<\mathrm{\hspace{0.17em}1}\mathrm{GeV}^2`$, $`0.1<y<\mathrm{\hspace{0.17em}0.8}`$, $`p_{}^\mu >\mathrm{\hspace{0.17em}2.0}\mathrm{GeV}`$, $`35^{}<\theta ^\mu <\mathrm{\hspace{0.17em}130}^{}`$ : $`\sigma _{vis}(epb\overline{b}X\mu X)=(170\pm 25)\mathrm{pb}`$.

The NLO QCD prediction, using the fixed-order NLO calculation  and a fragmentation parameterization described in , is $`(54\pm 9)\mathrm{pb}`$, where the error was estimated by varying the renormalization and factorization scales, the beauty quark mass value and the fragmentation parameters. The NLO prediction undershoots the measured cross section significantly.

The ZEUS collaboration has recently performed the first measurement of the beauty photoproduction differential cross sections . Events with a muon and at least two jets were selected. Differential cross sections as a function of the muon pseudorapidity, $`\eta ^\mu `$, and transverse momentum, $`p_T^\mu `$, were calculated. The beauty fraction in each bin of the selected variable was extracted with the $`p_T^{\mathrm{rel}}`$ fit.

Fig. 4 shows the measured beauty differential cross sections in comparison with predictions of PYTHIA  and HERWIG . The PYTHIA predictions are in reasonable agreement with the data. The agreement is worse in the most-forward $`\eta ^\mu `$ bin, in which the contribution from the beauty excitation processes is expected to be large. As shown, the beauty excitation is a substantial component of the PYTHIA cross section. The HERWIG predictions, which also include the beauty excitation component, are lower than PYTHIA but still compatible with the data within errors.

## 4 Summary

Charm and beauty photoproduction cross sections have been measured by the H1 and ZEUS collaborations. The perturbative QCD calculations are generally smaller than the measured cross sections, in particular in the forward (proton) direction. In the charm sector, theoretical uncertainties are larger than the experimental ones. The study of dijet photoproduction associated with charm is consistent with a significant contribution of charm excitation processes.


text:
# Topological Vacuum of the Closed Universe as a Gauging Factor

## 1 Introduction

Quantization of physical fields in a gravitational field of the Universe encounters the standard problem of vacuum divergencies. As, however, in quantum geometrodynamics (QGD) the standard perturbation theory with its renormalization procedure is inapplicable, the only known base for solution of the problem is supersymmetry  (its urgency for QGD at present is a commonplace). The vacuum contributions of fermion and boson fields are opposite on a sign, and the topological effect of a closed universe – discretization of the vacuum fluctuations spectra – allows to formulate conditions of cancellation of the vacuum divergencies. The most remarkable thing is that these conditions sharply limit gauge arbitrariness in QGD, in particular, the homogeneous isotropic model turns out completely deprived of this arbitrariness by them. Moreover, if, taking into account the known ambiguity in ordering and parametrization  of the Wheeler – DeWitt (WDW) theory , one imposes on it the mentioned conditions, the Hamiltonian of this theory will gain a nontrivial spectrum of eigenvalues, in essence conterminous with an energy spectrum of matter including its vacuum component. This radically changes the structure of the theory, approximating it to the dynamic one proposed in the work . Below, for distinguishing the two versions of geometrodynamics, we understand by QGD the dynamic theory , and the WDW structure, in which there is no time as a formal parameter of dynamics, we shall name, as it is accepted now, the quantum cosmology (QC WDW).

In the following section the statement of the problem in the simplest, semiclassical, version of a closed homogeneous isotropic universe is adduced. In Sec. 3. we shall formulate the problem of topological vacuum spectrum in the closed isotropic Universe, give the outcomes of its solution and show that the effect of topological vacuum (TV) can not be coordinated with the semiclassical theory. Sections 4. and 5. are devoted to solution of the problem within the framework of quantum theory: QC WDW (Sec. 4.) and QGD (Sec. 5.).

## 2 Semiclassical aspect

In the self-consistent semiclassical theory the $`\left\{\genfrac{}{}{0pt}{}{0}{0}\right\}`$-Einstein equation of a closed homogeneous isotropic universe at the presence of matter has the form:

$$\frac{3}{\kappa }\left(\frac{\dot{a}^2}{𝒩^2a^2}+\frac{1}{a^2}\right)=\frac{𝒩}{2\pi ^2a^5}H_{mat},$$
(1)
where $`\kappa `$ is the gravitational constant, $`𝒩=\sqrt{g_{00}}`$ is the “laps function”, $`a`$ is the scale factor, $`H_{mat}`$ is the quantum average Hamiltonian of material fields, as which we shall consider massless conformal two-component objects: neutral vector, spinor and charged scalar fields:

$$H^{\left(s\right)}=\underset{𝐤}{}H_𝐤^{\left(s\right)},$$

$`s`$ beeng the spin quantum number,

$$𝐤=(\lambda ,j,m),\lambda =1+s,2+s,\mathrm{},\mathrm{},j=s,s+1,\mathrm{},\lambda 1,m=j,\mathrm{},+j.$$

The Hamiltonian of boson fields can be presented as

$$H_𝐤^{\left(0\right)},H_𝐤^{\left(1\right)}=\lambda _𝐤\left(a_𝐤^{}a_𝐤+\overline{a}_𝐤^{}\overline{a}_𝐤+1\right),$$

where the two annihilation operators $`a_𝐤,\overline{a}_𝐤`$ concern two different components of a field (accordingly – the creation operators $`a_𝐤^{},\overline{a}_𝐤^{}`$).

The Hamiltonian of fermions

$$H_𝐤^{\left(1/2\right)}=\lambda _𝐤\left(b_𝐤^{}b_𝐤+\overline{b}_𝐤^{}\overline{b}_𝐤1\right).$$

In a state with definite energy

$$H^{\left(s\right)}=^{\left(s\right)}=\underset{𝐤}{}_𝐤^{\left(s\right)},_𝐤^{\left(s\right)}=\lambda _𝐤^{\left(s\right)}\left(n_{𝐤1}^{\left(s\right)}+n_{𝐤2}^{\left(s\right)}\pm 1\right),$$
(2)
where “$`\pm `$”, correspondingly, for bosons and fermions, $`n_𝐤`$ are the occupation numbers. According to (2), the quantity $`^{\left(s\right)}`$ consists of two physically different parts,

$$^{\left(s\right)}=_{part}^{\left(s\right)}+_{vac}^{\left(s\right)},$$

$$_{part}^{\left(s\right)}=\underset{𝐤}{}\lambda _𝐤^{\left(s\right)}\left(n_{𝐤1}^{\left(s\right)}+n_{𝐤2}^{\left(s\right)}\right)$$

is the particle energy,

$$_{vac}^{\left(s\right)}=\pm \underset{𝐤}{}\lambda _𝐤^{\left(s\right)}=\pm \underset{\lambda =1+s}{\overset{\mathrm{}}{}}\underset{j=s}{\overset{\lambda 1}{}}\underset{m=j}{\overset{j}{}}\lambda =\pm \underset{\lambda =1+s}{\overset{\mathrm{}}{}}\lambda \left(\lambda ^2s^2\right)$$
(3)
is the vacuum component. The equation (1) takes the form:

$$\frac{3}{\kappa }\left(\frac{\dot{a}^2}{𝒩^2a^2}+\frac{1}{a^2}\right)=\frac{𝒩}{2\pi ^2a^5},$$

$`=_s\left(_{part}^{(s)}+_{vac}^{(s)}\right)`$. In particular, in the gauge of “conformal time”

$$𝒩=a$$
(4)
we have

$$\frac{6\pi ^2}{\kappa }\left(\dot{a}^2+a^2\right)=.$$
(5)
As was mentioned in Introduction, the diverging vacuum contributions (3) of various fields can not be renormalized by the standard quantum field methods. However opposition of signs of these contributions from bosons and fermions specifies that the cancellation of divergencies is possible in the framework of the boson-fermion supersymmetry .

## 3 Topological vacuum. Inconsistency of semiclassical approximation

Let us consider an elementary supermultiplet of the two-component boson field and the two-component fermion one. According to (3), the vacuum energy of $`N`$ scalar-spinor supermultiplets is determined by the

$$_{vac}=\underset{a=1}{\overset{2N}{}}\underset{\lambda =1}{\overset{\mathrm{}}{}}\lambda ^3\underset{a=1}{\overset{2N}{}}\underset{\lambda =3/2}{\overset{\mathrm{}}{}}\lambda \left(\lambda ^2\frac{1}{4}\right).$$
(6)
To reduce the sums to common limits let us make in the spinor sum the replacement $`\lambda \lambda +1/2`$, and then in both sums $`\lambda ^3`$ put shifts $`\lambda =k+\alpha _a`$ with $`\alpha _a`$ integers:

$$\underset{a}{}\underset{\lambda =1}{\overset{\mathrm{}}{}}\lambda ^3=\underset{a}{}\underset{k=1}{\overset{\mathrm{}}{}}\left(k+\alpha _a\right)^3.$$

The shifts $`\alpha _a`$ should be chosen so that the divergencies were mutually compensated, and then the residual $`_{vac}`$ will develop out of the sums $`_{k=1\alpha _a}^0`$ if $`\alpha _a>0`$, and $`_{k=1}^{\alpha _a}`$ if $`\alpha _a<0`$. The degrees of the divergencies are determined by the formulas:

$$\underset{k=1}{\overset{K}{}}k^3=\frac{1}{4}K^2(K+1)^2;$$

$$\underset{k=1}{\overset{K}{}}k^2=\frac{1}{6}K\left(K+1\right)\left(2K+1\right);$$

$$\underset{k=1}{\overset{K}{}}k=\frac{1}{2}K\left(K+1\right).$$

So,

$$\frac{1}{6}AK\left(K+1\right)\left(2K+1\right)+\frac{1}{2}BK\left(K+1\right)+CK=0;$$

$$A=3\underset{a=1}{\overset{N}{}}\left(\alpha _a\beta _a\right)3\frac{N}{2},$$

$$B=3\underset{a=1}{\overset{N}{}}\left(\alpha _a^2\beta _a^2\right)\frac{N}{2},$$

$$C=\underset{a=1}{\overset{N}{}}\left(\alpha _a^3\beta _a^3\right),$$

$`\alpha _a`$ are the shifts in the scalar sums, $`\beta _a`$ are those in the spinor sums.

The condition of divergency cancellation is $`A=B=C=0`$, that means the set of equations:

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a\beta _a\right)=\frac{N}{2},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^2\beta _a^2\right)=\frac{N}{6},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^3\beta _a^3\right)=0.$$

The second version of divergency cancellation appears after the replacement $`\lambda \lambda 1/2`$ in the spinor sum (6). In this case the set of equations for $`\alpha _a,\beta _a`$ takes the form

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a\beta _a\right)=\frac{N}{2},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^2\beta _a^2\right)=\frac{N}{6},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^3\beta _a^3\right)=0.$$

For $`N`$ vector-spinor supermultiplets there are 4 different versions in correspondence with the two possible appropriate replacements of $`\lambda `$ in each of the two sums – boson and fermion:

1. $$\underset{a=1}{\overset{N}{}}\left(\alpha _a\beta _a\right)=\frac{N}{2},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^2\beta _a^2\right)=\frac{N}{2},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^3\beta _a^3\right)=0.$$
2. $$\underset{a=1}{\overset{N}{}}\left(\alpha _a\beta _a\right)=3\frac{N}{2},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^2\beta _a^2\right)=\frac{N}{2},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^3\beta _a^3\right)=0.$$
3. $$\underset{a=1}{\overset{N}{}}\left(\alpha _a\beta _a\right)=\frac{N}{2},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^2\beta _a^2\right)=\frac{N}{2},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^3\beta _a^3\right)=0.$$
4. $$\underset{a=1}{\overset{N}{}}\left(\alpha _a\beta _a\right)=\frac{N}{2},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^2\beta _a^2\right)=\frac{N}{2},$$

$$\underset{a=1}{\overset{N}{}}\left(\alpha _a^3\beta _a^3\right)=0.$$

These sets of equations are solvable at $`N`$ multiple 12, and the vector-spinor sets 1. and 2. are solvable at $`N`$ multiple 8 as well. The vacuum spectrum $`_{vac}`$ obtained as a result of solution of these sets of equations for various $`N`$ has proved to be equidistant, extending from $`\mathrm{}`$ up to $`+\mathrm{}`$ with the same pitch equal 6:

$$_{vac}=6n+v,n=0,\pm 1,\pm 2,\mathrm{},$$
(7)
and $`v`$ can take values $`0,\pm 2`$ depending on a set of supermultiples.

Coming back to the semiclassical equation (5), we must note the following. The $`_{part}`$ being a function of the occupation numbers, the $`_{vac}`$ is an independent variable, and the negative half of its rigid discrete spectrum cannot be coordinated with the positive definite left part of the equation. So, the semiclassical approach to the physics of the early Universe is incomplete: the presence of the TV demands appealing to the full quantum theory.

## 4 The Wheeler – DeWitt quantum cosmology.

The WDW equation

$$H\mathrm{\Psi }=0$$
(8)
can be arranged in such a way that the operator $`H`$ divides into tree parts,

$$H=H_{grav}+H_{part}+_{vac},$$

where $`H_{grav}`$ is its gravitational component, $`H_{part}=_i_𝐤H_𝐤^{(i)},H_𝐤^{(i)}=\lambda _𝐤^{(i)}a_𝐤^ia_𝐤^i`$, the index $`i`$ numbers components of material fields, $`a_𝐤^i`$ are the annihilation operators of the appropriate particles. The solution of the equation (8) can be factored: $`\mathrm{\Psi }=\mathrm{\Psi }_{part}\mathrm{\Psi }_{grav}`$, so that

$$H_{part}\mathrm{\Psi }_{part}=_{part}\mathrm{\Psi }_{part},$$
(9)
$$H_{grav}\mathrm{\Psi }_{grav}=(_{grav}_{part}_{vac})\mathrm{\Psi }_{grav}.$$
(10)
Thus, an eigenvalue $`_{part}`$ of the operator $`H_{part}`$ has to be considered as a parameter of the equation (10). As to the spectrum $`_{vac}`$, it should be obtained as a result of solution of the WDW equation responsible for a state of the Universe as a whole. Fortunately, parametrization/gauge and permutation ambiguities of this equation allow to come near the solution of the problem.

At first, the oscillator aspect of the left part of the equation (5) and the spectrum nature of the$`_{vac}`$ (7) specify the definite advantage of the gauge (4). If, in addition, one puts $`_{part}=0`$ and chooses a system of units in which $`6\pi ^2/\kappa =1/2`$, the appropriate WDW equation can be presented as:

$$\frac{1}{2}\left(\frac{^2}{a^2}a^2\right)\mathrm{\Psi }_{vac}=E_{vac}\mathrm{\Psi }_{vac},a0,$$

$$E_{vac}=2m+3/2,m=0,1,2,\mathrm{}.$$

Comparing the spectrum of eigenvalues of this equation with (7), we should, first of all, remove all the values of the $`m`$ that are not multiple 3. Existing $`_{part}0`$ requires presence at the left part of the equation of an element that could absorb this component. And here comes to the aid the second ambiguity of the quantum theory in general and of the WDW equation in particular, known as the ordering problem. Here it manifests itself in that in the general case the kinetic part of the WDW Hamiltonian can be presented as

$$\frac{1}{2}f^1\left(a\right)\frac{}{a}f^2\left(a\right)\frac{}{a}f^1\left(a\right),$$

where $`f(a)`$ is an arbitrary differentiable function of the scale factor. Let us put $`f=a^{l+1}`$. Then the WDW equation takes the form of the equation for a “centrifugal oscillator”

$$\left[\frac{1}{2}\left(\frac{^2}{a^2}a^2\frac{l\left(l+1\right)}{a^2}\right)+\right]\mathrm{\Psi }_{}=0$$

with the spectrum (comp. )

$$=2m+l+\frac{3}{2},l=_{part}.$$
(11)
What do we have in the result? Presence of the TV with its distinctive spectrum (7)

1. fixes the gauge and solves the ordering problem;
2. transforms QC WDW to the theory with a nontrivial spectrum.

Thus, the TV has properties intrinsic to the gravitational vacuum condensate (GVC) which existence the consecutive construction of QGD on the basis of the usual quantum theory  inevitably leads to.

On the other hand, the WDW equation cuts the vacuum spectrum removing, in particular, the negative values $`n`$, but does not contain any absorption mechanism for this part of the energy spectrum of the subsystem “topological vacuum”. This is a superfluous testimony of the stated in  incompleteness of the theory.

The other situation is in QGD.

## 5 Quantum geometrodynamics

In the dynamic Schrödinger equation for the physical part of a wave function in the gauge (4) the Hamiltonian has the form:

$$H_{ph}=\frac{1}{8}\frac{1}{M\left(a\right)}\frac{}{a}M\left(a\right)\frac{}{a}2a^2,$$
(12)
$`M`$ being a probability measure. Let us put

$$M=a^{2\left(l+1\right)},l0.$$
(13)
The probability measure, beeng as well a measure of a path integral (PI), is generated by the certain dependence on $`a`$ of a pitch $`ϵ`$ of splitting up the time interval at a final PI definition through a multiple integral, as the measure is connected with $`ϵ`$ by the relation , which in this case has the form

$$\gamma ^{\frac{1}{2}}M=\mathrm{const},\gamma =a^{4(l+1)},$$

where $`\gamma `$ is defined by

$$ϵ=\chi \gamma ,\chi 0.$$

By substituting the measure (13) into the Hamiltonian (12), we come to the equation for stationary states

$$\left[\frac{1}{8}\frac{^2}{a^2}+\frac{1}{4a}\frac{}{a}2a^2\frac{l\left(l+1\right)}{8a^2}+\right]\mathrm{\Psi }=E\mathrm{\Psi },$$

$`E`$ being a parameter of GVC . The eigenvalue spectrum of the equation is defined by the same expression (11),

$$E=_{part}+_{vac}E=2n+l+\frac{3}{2},n=0,1,2,\mathrm{};l=_{part},$$

but this time there is the GVC parameter $`E`$ allowing to put it in agreement with the TV: presenting (7) in the form $`_{vac}=6(nm+v),n,m=1,2,\mathrm{}`$ we have

$$E=4n6m+6v\frac{3}{2}.$$
(14)
Thus, GVC and TV represent a united subsystem, breaking gauge symmetry concerning the diffeomorphism group.

## 6 Conclusion.

The above-stated analysis shows that presence of material fields with physical topological vacuum breaks the basic equation of the WDW theory – the quantum version of the classical Hamilton constraint, – revealing GVC. This confirms our conclusion made in the work  that the quantum principle in cosmology is incompatible to the principle of gauge invariance concerning the diffeomorphism group and that QC WDW has not any general-theoretical foundation.

Here we have considered the elementary cosmological model, which, however, most likely can serve a good “zero approximation” to the exact model of the early Universe. It definitely allows to expect that the obtained results, though will not be so hard for the real model, nevertheless, in essential will be applicable to it. Or else, it is possible to state, that TV considerably limits the gauge/parametrization arbitrariness and the ordering ambiguity (choice of a probability measure) in QGD.


text:
# The Diffractive Logarithmic Slope and the Saturation Phenomena

## 1 Introduction

The measurements of the derivative quantity $`F_2`$-slope on $`Q^2`$ allowed a renewed interest of testing the matching of hard and soft approaches and provided constraints for the saturation formalisms. The reported turnover on the $`x`$ dependence was readly associated with the transition region between the interplaying domains. The most recent determinations of the slope still deserve a better theoretical description, keeping a turnover pattern at fixed c.m.s. energy $`W`$. When we focus in diffractive DIS, in particular the structure function $`F_2^D`$, the situation is far from clear: initially considered as a predominantly soft process, the experimental results suggests that the diffractive cross section at HERA contains hard and soft components. However, the hard piece is stronger than the previously expected. Fortunately, if we consider the pQCD approach, diffraction stands a more profitable field to study saturation effects than the inclusive case. This comes from the fact that in DDIS the interaction probes larger dipole size configuration (soft content) than in the DIS reaction. Although the quite different pictures considered to interpret the diffractive measurements, almost all of them fit the data set properly. We have proposed a derivative quantity, the diffractive logarithmic slope, which would help to disentangle the underlying dynamics in diffractive DIS, settling its validity range, if such observable is measured. In Ref. one perform studies of the proposed quantity for two sound models, representing the essential features in both Regge and pQCD formalisms. It was found that important deviations in the predictions between the models emerge, considering the available kinematic spectrum in diffractive DIS . Here we report in particular the calculations considering the saturation pQCD model applied to diffraction dissociation and its comparison with the non-satured QCD approach, performed in details in. Such results allow to verify the role played by the saturation effects in these processes and the extension of the standard QCD formalism to kinematical regions forbidden perturbatively.

## 2 The Diffractive Logarithmic Slope and the Saturation Phenomenon

In perturbative QCD, the $`\gamma ^{}p`$ process is described in terms of the photon splitting into a $`q\overline{q}`$ pair, far upstream of the nucleon, which then scatters in the proton. This reaction is mediated by the one gluon exchange which turns out into a multi-gluon one when the saturation region is approached. A remarkable feature is that the mechanism leading to the photon dissociation and the further scattering has a factorizable form, written in terms of a convolution between a photon wavefunction (dipole transverse size $`r`$ and quark longitudinal momentum fraction $`z`$) and a $`q\overline{q}`$ cross section. For transverse and longitudinally polarized photons the $`\gamma ^{}p`$ cross section has the form

$`\sigma _{T,L}(x,Q^2)={\displaystyle d^2𝐫_0^1𝑑\alpha |\mathrm{\Psi }_{T,L}(\alpha ,𝐫)|^2\widehat{\sigma }(x,r^2)},`$ (1)
where $`\mathrm{\Psi }_{T,L}`$ are the photon wavefunctions. Although they are calculated from perturbation theory, the $`q\overline{q}`$ cross section, $`\widehat{\sigma }`$, has strong nonperturbative contributions that should be modeled. We choose the saturation model, which reproduces the experimental results at both inclusive and diffractive electroproduction. The dynamics of saturation is present in the effective dipole cross section in an eikonal way $`\widehat{\sigma }(x,r^2)=\sigma _0\left[\mathrm{\hspace{0.17em}1}\mathrm{exp}\left(\frac{r^2}{4R_0^2(x)}\right)\right]`$, where the $`x`$-dependent saturation scale is $`R_0(x)=\frac{1}{Q_0}\left(\frac{x}{x_0}\right)^{\lambda /2}`$. The normalization $`\sigma _0`$ and the remaining parameters were determined from data with $`x<10^2`$. This model consistently interpolates between the saturation regime and the scaling regime of $`\sigma _T(F_2)`$. The elastic scattering of the $`q\overline{q}`$ pair dominates the diffractive $`\gamma ^{}p`$ process for not too large values of the diffractive mass $`M_X`$. Instead, at large $`M_X`$ the emission of a gluon becomes the leading contribution. Then, the diffractive cross section is, given the slope $`B_D`$,

$`\sigma ^D(x,Q^2)={\displaystyle \frac{1}{B_D}}{\displaystyle \frac{1}{16\pi }}{\displaystyle d^2𝐫_0^1𝑑\alpha |\mathrm{\Psi }_{T,L}(\alpha ,𝐫)|^2\widehat{\sigma }^2(x,r^2)}.`$ (2)
In ref. it has been found a ratio between the diffractive and the inclusive cross sections roughly constant, and that the diffractive cross section is sensitive to the infrared cutoff which is effectively provided by $`2R_0(x)`$ concluding that diffraction probes the transition region.

We performed numerically the $`F_2^D`$ logarithmic slope, using the same parameters values from Ref.. In Fig. (1), one shows the $`x_{IP}`$ behavior for both the approaches with and without saturation, at typical $`\beta `$ values. We analyze in particular the transition region between hard and soft dynamics, settled by the low virtualities $`Q^21.59`$ GeV2. The saturation model produces a transition between positive and negative slope values at low $`\beta =0.04`$, while presents a positive slope for medium and large $`\beta `$. This feature is corroborated by the preliminary ZEUS studies. The pQCD approach without saturation, shows a positive slope for the whole $`Q^2`$ and $`x_{IP}`$ range. In , one observed such a smooth behavior for the pQCD predictions, with a positive slope in almost the whole range, whereas a sound Regge model predicted predominantly negative slope values in this kinematical domain, converging to a flat value for larger $`x_{IP}`$. A clear difference between the predictions of Regge and QCD is the change of signal of the slope with the $`Q^2`$ evolution at medium $`\beta `$. Returning, since the diffractive cross section is strongly sensitive to the infrared cutoff, one of the main differences between these models is the assumption related to the small $`Q^2`$ region. In the pQCD models without saturation, an ad hoc cutoff in the transverse momentum is inserted, as well as the energy dependence of the unintegrated gluon distribution. In the saturation model, instead, the saturation radius $`R_0(x)`$ gives the infrared cutoff (the saturation momentum scale) and determines the energy dependence. If this scale is large (1-2 GeV2), then the resulting process is not soft and can be completely calculated using pQCD methods. Therefore, the saturation model extends the pQCD approach towards lower $`Q^2`$ values. We conclude that, the difference between the behaviors predicted by these two models for the $`x_{IP}`$ spectrum, mainly in the region of small $`\beta `$ and medium $`Q^2`$, is large, which should allow to discriminate the dynamics in future experimental analyzes.