text:
# References

Differential constraints and exact solutions
of nonlinear diffusion equations

Oleg V. Kaptsov and Igor V. Verevkin

Institute of Computing Modeling RAS, Academgorodok,

660036, Krasnoyarsk, Russia

E-Mail: kaptsov@ksc.krasn.ru

Abstract

The differential constraints are applied to obtain explicit solutions of nonlinear diffusion equations. Certain linear determining equations with parameters are used to find such differential constraints. They generalize the determining equations used in the search for classical Lie symmetries.

PACS numbers: 02.30.Jr, 02.30.Ik, 44.05.+e

Mathematics Subject Classification: 34G20, 35K57

1. Introduction.

Differential constraints arisen originally in the theory of partial differential equations of the first order. In particular Jacobi used differential constraints to find the total integral of nonlinear equation

$$F(x_1,\mathrm{},x_n,z,z_{x_1},\mathrm{},z_{x_n})=0,$$

König applied them to the equation of the second order . They required that the corresponding over-determined system was compatible. The general theory of over-determined systems was developed by Delassus, Riquier, Cartan, Ritt, Kuranishi, Spencer and others. One can find references in the book of Pommaret . Now the applications of over-determined systems include such diverse fields as differential geometry, continuum mechanics and nonlinear optics.

The general formulation of the method of differential constraints requires that the original system of partial differential equations

$$F^1=0,\mathrm{},F^m=0$$
$`\left(1\right)`$
be enlarged by appending additional differential equations (differential constraints)

$$h_1=0,\mathrm{},h_p=0,$$
$`\left(2\right)`$
such that the over-determined system (1), (2) satisfies some conditions of compatibility.

One can derive many exact solutions of partial differential equation by means of differential constraints. It was particularly shown in  that some soliton solutions can be found using differential constraints. Olver and Rosenau , Olver , Kaptsov , Levi and Winternitz  show that many reduction methods such as non-classical symmetry groups, partial invariance, separation of variables, the Clarkson–Kruskal direct method can be included into the method of differential constraints. In practice, methods based on the Riquier–Ritt theory of over-determined systems of partial differential equations may be difficult. The problem of finding all differential constraints compatible with certain equations can be more complicated than the investigation of the original equations.

It was recently proposed a new method for finding differential constraints which uses linear determining equations. These equations are more general than the classical determining equations for Lie generators  and depend on some parameters. Given an evolution equation

$$u_t=F(t,x,u,u_1,\mathrm{},u_n),$$
$`\left(3\right)`$
where $`u_k=\frac{^ku}{x^k}`$, then according to  the linear determining equation corresponding to (3) is of the form

$$D_t\left(h\right)=\underset{i=0}{\overset{N}{}}\underset{k=0}{\overset{i}{}}b_{ik}D_x^{ik}\left(F_{u_{Nk}}\right)D_x^{Ni}\left(h\right),b_{ik}R.$$
$`\left(4\right)`$
Here and throughout $`D_t,D_x`$ are the operators of total differentiation with respect to $`t`$ and $`x`$. Equality (4) must hold for all solutions of (3). The function $`h`$ may depends on $`t,x,u,u_1,\mathrm{},u_p`$. The number $`p`$ is called the order of the solution of equation (4). If we have some solution $`h`$ then corresponding differential constraint is

$$h=0.$$
$`\left(5\right)`$
It was also shown in  that equations (4) and (5) constitute the compatible system. Thus we sketch the derivation of some solutions to the evolution equation (4).
(I) Find solutions of the linear determining equations (4).
(II) Fixing the function $`h`$, we obtain differential constrain (5).
(III) Find the general solution of (5) which includes some arbitrary functions $`a_i`$ depending on $`t`$.
(IV)Substitute the general solution into (4). It leads to ordinary differential equations for functions $`a_i`$.
(V) Solve the ordinary differential equations and obtain a solution of the evolution equation (4).

In this paper we start with determination of the solutions of linear determining equations of the second and third orders for the nonlinear diffusion equation

$$u_t=\left(u^ku_x\right)_x+f\left(u\right).$$
$`\left(6\right)`$
These solutions exist only if $`f`$ belongs to the special forms. Then we use the obtained functions $`h`$ to find solutions of equation (6). In final section we derive exact solutions of two-dimension equation

$$u_t=\mathrm{\Delta }\mathrm{ln}\left(u\right).$$

2. Solutions of linear determining equations.

The nonlinear diffusion equation

$$u_t=\left(Q\left(u\right)u_x\right)_x+f\left(u\right).$$
$`\left(7\right)`$
often arises in the description of various physical processes. The group classification of the equation has been carried out in . Some exact solutions of (7) can be found in . In physical applications $`Q`$ is usually taken to be a power function. In this section we consider the equation

$$u_t=\left(u^qu_x\right)_x+f\left(u\right),$$
$`\left(8\right)`$
where $`f`$ is an differentiable function, $`q0`$. If $`q=2`$, $`f=u`$ or $`f=const`$ then the equation (8) can be linearized. We shall not discuss this case here.

The linear determining equation, which corresponds to (8), is

$$D_th=u^qD_x^2h+b_1qu_xu^{q1}D_xh+\left(b_3qu^{q1}u_{xx}+b_2q\left(q1\right)u^{q2}u_x^2+b_4f_u\right)h,$$
$`\left(9\right)`$
where $`b_1,\mathrm{},b_4R.`$ We shall seek solutions to (9) in the form

$$h=u_n+g(t,x,u,\mathrm{},u_{n1}),$$

where $`n2,`$ $`u_k=\frac{^ku}{x^k}.`$ The method for finding solutions is very similar to the standard procedure applied in the group analysis of differential equations  and only one of all possibilities is described here for the sake of brevity.

We set $`n=2.`$ First, let us express all $`t`$-derivatives in (9) using (8). As a result, the left-hand side of (3.4) becomes a polynomial with respect to $`u_3,u_2`$. The polynomial must identically vanish. Collecting similar terms we obtain the following relations for the coefficients of $`u_3`$ and $`u_2^2`$

$$q\left(b_14\right)=0,ug_{u_1u_1}+q\left(b_33\right)=0.$$

Thus $`b_1=4`$ and $`g`$ can be represented as follows

$$g=\frac{\left(3b_3\right)q}{2u}u_1^2+a(u,t,x)u_1+g_1(u,t,x);$$

here $`a`$ and $`g_1`$ must be functions of $`u,t`$ and $`x`$ alone. Collecting the coefficients of $`u_2u_1^2`$ and $`u_2u_1`$, we have the equations

$$2b_2q2b_2b_3^2q+b_3q+4b_36q=0,$$
$`\left(10\right)`$
$$2ua_u+q\left(b_3+1\right)a=0.$$

In follows from the last equation that

$$a=a_1(t,x)u^{\frac{\left(1+b_3\right)}{2}q},$$

where $`a_1`$ is a function of $`t`$ and $`x`$. Next we consider the coefficient $`u_1^3`$ and obtain equation

$$4b_2q4b_2+b_3^2q4b_3q+2b_39q+6=0.$$
$`\left(11\right)`$
From (10) and (11) it follows that $`b_3=1`$ or $`b_3=\frac{q+2}{q}.`$

Assuming $`b_3=1`$, we obtain $`b_2=\frac{3q2}{q1}.`$ The coefficient of $`u_2`$ give us equation

$$u^{q+2}\left(2a_{1_x}+f_u\left(b_41\right)\right)+u^{2q+1}qg_1=0$$

The equation enables us to express

$$g_1=\frac{1}{q}u^{1q}\left(f_u\left(1b_4\right)2a_{1_x}\right)$$

The coefficient of $`u_1^2`$ yields Euler equation

$$u^3\left(1b_4\right)f_{uuu}+u^2\left(2qb_42b_4\right)f_{uu}uq^2f_u+q^2f=0.$$

Consider for simplicity the case $`b_4=1`$. It is easy to see that the last equation has two types of solutions:

$$f=ku+nu^q,q1$$

or

$$f=ku+nu\mathrm{ln}u,q=1$$

where $`k,n`$ are arbitrary constants. Let us focus on $`f=ku+nu^q`$. It follows from above calculations and equation (9) that

$$u^{q+1}\left(qa_{1_t}3u^qqa_{1_{xx}}4u^qa_{1_{xx}}+kq^2a_1\right)u_1+$$

$$+nq^2a_1u_1+2u^{q+2}\left(a_{1_{tx}}u^qa_{1_{xxx}}kqa_{1_x}\right)+2una_{1_x}=0.$$
$`\left(12\right)`$
From (12) we have $`na_1=0`$.

If $`a_1=0`$ then the solution of (9) is $`h=u_2+qu_1^2/u`$. If $`n=0`$ and $`q4/3`$ then one easily computes

$$a_1=\left(rx+s\right)\mathrm{exp}\left(kqt\right),$$

$$h=u_2+q\frac{u_1^2}{u}+\left(\left(rx+s\right)u_1\frac{2}{q}ur\right)u^q\mathrm{exp}\left(kqt\right).$$

In the case $`q=4/3`$, we obtain

$$a_1=\left(rx^2+sx+p\right)\mathrm{exp}\left(\frac{4}{3}kt\right),$$

$$h=u_2\frac{4u_1^2}{3u}+\left(\left(rx^2+sx+p\right)u_1+\frac{3}{2}u\left(2rx+s\right)\right)u^{4/3}\mathrm{exp}\left(\frac{4}{3}kt\right).$$

It can be shown that the found functions $`h`$ lead to invariant solutions of the corresponding equation (8).

We omit here for the sake of brevity intermediate calculations and give the list of solutions to the equation (9):
(1) if $`q=1`$ and $`f=su+ru\mathrm{ln}(u)`$ then

$$h=u_2\frac{u_1^2}{u};$$

(2) if $`q1`$ and $`f=su+ru^q`$ then

$$h=u_2+\frac{qu_1^2}{u};$$

(3) if $`q=2`$ and $`f=su+ru^3`$ then

$$h=u_2\frac{3u_1^2}{2u};$$

(4) if $`q=1`$ and $`f=ru`$ then

$$h=u_2+s\mathrm{exp}\left(rt\right)u^2u_1+r/3;$$

(5) if $`q`$ is an arbitrary constant and $`f=su+ru^{1q}`$ then

$$h=u_2\frac{\left(q1\right)u_1^2}{u},$$

with $`r,sR.`$
We did not include functions $`h`$ that correspond to invariant solutions of the equation (8).

If we will look for solutions to the equation (9), which depend on third derivative, then obtain the following list:
(1) if $`q`$ is an arbitrary constant and $`f=su+ru^{1q}+\frac{n(q+1)}{q^2}u^{q+1}`$ then

$$h=u_3+\frac{3\left(q1\right)}{u}u_1u_2+\left(q^23q+2\right)\frac{u_1^3}{u^2}+nu_1;$$

(2a) if $`q1`$ and $`f=nu+\frac{r}{q}u^{q+1}`$ then

$$h=u_3+\frac{\left(3q1\right)}{u}u_1u_2+q\left(q2\right)\frac{u_1^3}{u^2}+ru_1;$$
$`\left(13\right)`$
(2b) if $`q=2`$ or $`q=4/3`$ and $`f=nu+\frac{r}{q}u^{q+1}+mu^{q+3}`$ then $`h`$ is also given by (13);
(3) if $`q=\frac{1}{2}`$ and $`f=mu`$ then

$$h=u_3\frac{5u_1u_2}{2u}+\frac{5u_1^3}{4u^2}+r\mathrm{exp}\left(3mt/2\right)u^{5/2}+s\mathrm{exp}\left(mt/2\right)u^{1/2};$$

(4) if $`q=\frac{3}{2}`$ and $`f=nu+mu^{5/2}`$ then

$$h=u_3\frac{15u_1u_2}{2u}+\frac{35u_1^3}{4u^2}+r\mathrm{exp}\left(3nt/2\right)u^{5/2};$$

(5) if $`q=\frac{1}{2}`$ and $`f=mu2ku^{1/2}`$ then

$$h=u_3\frac{5u_1u_2}{2u}+\frac{5u_1^3}{4u^2}+ku_1+s\mathrm{exp}\left(mt/2\right)u^{1/2};$$

(6) if $`q=\frac{3}{2}`$ and $`f=nu`$ then

$$h=u_3\frac{15u_1u_2}{2u}+\frac{35u_1^3}{4u^2}+s\mathrm{exp}\left(7n/2t\right)u^{9/2}+r\mathrm{exp}\left(3nt/2\right)u^{5/2};$$

(7) if $`q=1`$ and $`f=mu`$ then

$$h=u_3\frac{4u_1u_2}{u}+\frac{3u_1^3}{u^2}+s\mathrm{exp}\left(2mt\right)u^2u_1;$$

with $`r,s,m,nR.`$ Here we also did not include functions $`h`$ leading to invariant solutions of (8).

3. Solutions of diffusion equations.

In this section we shall use the functions obtained above to construct solutions of diffusion equations (8). One can apply the method described in introduction.

We first take the function $`h=u_2+qu_1^2/u`$, where $`qR`$, corresponding to some cases mentioned above. Simply by equating this function to zero, we obtain the differential constraint

$$u_2+qu_1^2/u=0.$$
$`\left(14\right)`$
The equation (14) has two types of solutions:

$$u=\left(c_1x+c_2\right)^{\frac{1}{q+1}},q1,$$
$`\left(15\right)`$
$$u=c_1\mathrm{exp}\left(c_2x\right),q=1,$$
$`\left(16\right)`$
where $`c_1,c_2`$ are functions of $`t`$.

If we substitute the representation (16) into equation

$$u_t=\left(u_x/u\right)_x+ku\mathrm{ln}u$$
$`\left(17\right)`$
then this leads us to differential equations for $`c_1,c_2`$. From this equations it is easy to find $`c_1,c_2`$ and obtain the following solution of (17)

$$u=s_1\mathrm{exp}\left(s_2x\mathrm{exp}\left(kt\right)\right),s_1,s_2R.$$

Substituting (15) into equation

$$u_t=\left(u^qu_x\right)_x+su+ru^q,$$

we find the solution

$$u=\mathrm{exp}\left(st\right)\left(ax+b\frac{r}{s\left(q+1\right)}\mathrm{exp}\left(s\left(q+1\right)t\right)\right)^{\frac{1}{q+1}},$$

with $`a,bR`$.

It is easy to see that the differential constraint

$$u_2+\left(q1\right)\frac{u_1^2}{u}=0$$

for equation

$$u_t=\left(u^qu_x\right)_x+ru+su^{1q}$$

leads to solution

$$u=\left(qax\mathrm{exp}\left(qrt\right)+\frac{a^2}{r}\mathrm{exp}\left(2qrt\right)\frac{s}{r}\right)^{\frac{1}{q}},aR.$$

Now let us consider the differential constraints of the third order. We start with the equation

$$u_t=\left(u^qu_x\right)_x+su+ru^{1q}+n\frac{q+1}{q^2}u^{q+1},nR$$
$`\left(18\right)`$
As explained above this equation is compatible with the differential constraint

$$u_3+3\left(q1\right)\frac{u_1u_2}{u}+\left(q^23q+2\right)\frac{u_1^3}{u^2}+nu_1=0.$$
$`\left(19\right)`$
By a change of variable $`v=u^q`$ one may rewrite (18), (19) in the following way

$$v_t=vv_{xx}+\frac{1}{q}v_x^2+n\frac{q+1}{q}v^2+sqv+rq,$$
$`\left(20\right)`$
$$v_3+nv_1=0.$$
$`\left(21\right)`$
If $`n=1`$ then it follows from (21) that

$$v(t,x)=a+b\mathrm{exp}\left(x\right)+c\mathrm{exp}\left(x\right),$$
$`\left(22\right)`$
where $`a,b`$ and $`c`$ are some functions of $`t`$. Substituting this representation into equation (20) we obtain the system of ordinary differential equations for the function $`a,b`$ and $`c`$:

$$a_t=a^2\left(1+1/q\right)+asq4bc/q+rq,$$
$`\left(23\right)`$
$$b_t=ab\left(1+2/q\right)+bsq,$$
$`\left(24\right)`$
$$c_t=ac\left(1+2/q\right)+csq.$$
$`\left(25\right)`$
Using finite-dimensional invariant subspaces, Galaktionov  found representation (22).

From (24) and (25) we derive the first integral $`b=kc,kR`$. Therefore the system (23)-(25) can be reduced to nonlinear ordinary differential equation for the function $`a`$. In general we can not express solutions of (23)-(25) in terms of the elementary functions. We give one example of exact solution of equation (20), with $`q=1`$ and $`r=0`$. This solution has the representation (22) and the functions $`a,b,c`$ are

$$a=a_1\mathrm{sin}\left(pr_1\mathrm{exp}\left(r_1t\right)+m\right)\mathrm{exp}\left(r_1t\right)/\mathrm{cos}\left(pr_1\mathrm{exp}\left(r_1t\right)+m\right),$$

$$b=a_2\mathrm{exp}\left(r_1t\right)/\mathrm{cos}\left(pr_1\mathrm{exp}\left(r_1t\right)+m\right),$$

$$c=a_3\mathrm{exp}\left(r_1t\right)/\mathrm{cos}\left(pr_1\mathrm{exp}\left(r_1t\right)+m\right),$$

where $`a_1=pr_1^2`$, $`a_3=p^2r_1^4/4a_2`$, $`r_1=s`$ and $`p,a_2,m`$ are arbitrary constants.

Now let us consider the equation

$$u_t=\left(u^{1/2}u_x\right)_x+mu2k\sqrt{u},m,kR$$
$`\left(26\right)`$
and the differential constraint

$$u_3\frac{5u_1u_2}{2u}+\frac{5u_1^3}{4u^2}+ku_1+se^{mt/2}\sqrt{u}=0.$$
$`\left(27\right)`$
Using the equation (26), one can write (27) as

$$\left(\mathrm{ln}u\right)_{tx}+se^{mt/2}u^1=0.$$

Replacing $`\mathrm{ln}(e^{mt/2}u^1)`$ by $`w`$, the last equation is replaced by the Liouville equation

$$w_{tx}=s\mathrm{exp}\left(w\right).$$

Since the general solution of the Liouville equation is

$$w=\mathrm{ln}\frac{2T^{}X^{}}{s\left(T+X\right)^2},$$

it gives the representation

$$u=\frac{s\left(T+X\right)^2}{2T^{}X^{}}e^{mt/2},$$
$`\left(28\right)`$
where $`T`$ and $`X`$ are the arbitrary functions of $`t`$ and $`x`$ respectively. Substituting this representation into (26), we have

$$\sqrt{\frac{s}{2}}e^{mt/4}\left(2\left(T^{}\right)^{1/2}\left(T+X\right)\left(T^{}\right)^{3/2}T^{\prime \prime }\frac{m}{2}\left(T+X\right)\left(T^{}\right)^{1/2}\right)=$$

$$=\frac{3}{2}\left(X^{}\right)^{3/2}\left(X^{\prime \prime }\right)^2\left(X^{}\right)^{1/2}X^{\prime \prime \prime }2k\left(X^{}\right)^{1/2}.$$
$`\left(29\right)`$
Differentiating (29) with respect to $`t`$, it is easy to obtain the equation for $`T`$

$$2T^{\prime \prime \prime }T^{}3\left(T^{\prime \prime }\right)^2+\frac{m^2}{4}\left(T^{}\right)^2=0.$$

If $`m0`$ then the function

$$T=c_1\mathrm{tanh}\left(\frac{mt}{4}+c_3\right)+c_2,c_1,c_2,c_3R$$

is the general solution of this equation. Substituting the function $`T`$ into (29), we get the following equation for $`X`$

$$\sqrt{\frac{sm}{2c_1}}\left(c_1c_2X\right)=\frac{3}{2}\left(X^{}\right)^{3/2}\left(X^{\prime \prime }\right)^2\left(X^{}\right)^{1/2}X^{\prime \prime \prime }2k\left(X^{}\right)^{1/2}.$$

It should be noted that (28) is equivalent to the representation

$$u=\left(a_1+a_2e^{mt/2}\right)^2,$$

where $`a_1,a_2`$ are functions of $`x`$. This representation yields the following system for $`a_1`$ and $`a_2`$

$$a_{1_{xx}}+a_1^2m/2a_1k=0,$$
$`\left(30\right)`$
$$a_{2_{xx}}+a_1a_2m/2a_2k=0.$$
$`\left(31\right)`$
In general, it is possible to express $`a_1`$ in terms of the Weierstrass function $`\mathrm{}`$ and $`a_2`$ in terms of Lamé’s function  . However, if $`m=12`$ and $`k=4`$ , then the functions

$$a_1=\frac{1}{\mathrm{cosh}^2\left(x\right)},$$

$$a_2=\frac{a}{\mathrm{cosh}^2\left(x\right)}+\frac{b}{\mathrm{cosh}^2\left(x\right)}\left(\frac{\mathrm{sinh}4x}{32}+\frac{\mathrm{sinh}2x}{2}+\frac{3x}{8}\right),a,bR$$

satisfy the equations (30), (31).

According to our results in the previous section, as $`k=0`$ , the equation (26) is compatible with the differential constraint

$$u_3\frac{5u_1u_2}{2u}+\frac{5u_1^3}{4u^2}+re^{3mt/2}u^{5/2}+se^{mt/2}\sqrt{u}=0.$$
$`\left(32\right)`$
Using the equation (26), one can write (32) as

$$\left(\mathrm{ln}u\right)_{tx}+re^{3mt/2}u+se^{mt/2}u^1=0.$$
$`\left(33\right)`$
Replacing $`\mathrm{ln}(e^{3mt/2}u)`$ by $`w`$ in (33) yields

$$w_{tx}+e^w+sre^{wmt}=0.$$

If we set $`s=0`$ then from the last equation we find the following representation

$$u=\frac{2X^{}T^{}}{\left(X+T\right)^2}e^{3mt/2},$$

where $`T`$ and $`X`$ are the arbitrary functions of $`t`$ and $`x`$ respectively. Substituting this representation into (26) leads to equation

$$\sqrt{2/r}e^{3mt/4}\left(m\left(T^{}\right)^{1/2}\left(X+T\right)2\left(T^{}\right)^{1/2}T^{\prime \prime }\left(X+T\right)+4\left(T^{}\right)^{3/2}\right)=$$

$$=2X^{\prime \prime \prime }\left(X^{}\right)^{3/2}\left(X+T\right)^2+8X^{\prime \prime }\left(X^{}\right)^{1/2}\left(X+T\right)8\left(X^{}\right)^{3/2}+\left(X^{\prime \prime }\right)^2\left(X^{}\right)^{5/2}\left(X+T\right)^2$$

Introducing new functions

$$C\left(T\right)=\sqrt{2/r}e^{3mt/4}\left(m\left(T^{}\right)^{1/2}+2T^{\prime \prime }\left(T^{}\right)^{1/2}\right),$$

$$B\left(X\right)=\left(X^{\prime \prime }\right)^2\left(X^{}\right)^{5/2}2X^{\prime \prime \prime }\left(X^{}\right)^{3/2},$$

$$D\left(X\right)=2XB+8X^{\prime \prime }\left(X^{}\right)^{1/2},$$

$$Q\left(T\right)=\sqrt{2/r}e^{3mt/4}\left(m\left(T^{}\right)^{1/2}T+2\left(T^{}\right)^{1/2}T^{\prime \prime }T4\left(T^{}\right)^{3/2}\right),$$

$$R\left(X\right)=BX^2+8X^{\prime \prime }\left(X^{}\right)^{1/2}X8\left(X^{}\right)^{3/2},$$

one can write the last equation as

$$C\left(T\right)X+D\left(X\right)T+B\left(X\right)T^2+Q\left(T\right)+R\left(X\right)=0.$$
$`\left(34\right)`$
It is possible consider (34) as condition of orthogonality of two vector functions $`Z=(C,T,T^2,Q,1),W=(X,D,B,1,R).`$

Denote by $`\rho (Z)`$ and $`\rho (W)`$ the number of linearity independent functions among $`C,T,T^2,Q,1`$ and $`X,D,B,R,1`$ respectively. From orthogonality condition it follows that $`\rho (Z)+\rho (W)5.`$ It is possible to show that if $`T^{}0`$ then $`\rho (Z)=3`$ and $`\rho (W)=2`$. In this case we have

$$D\left(X\right)=a_1X+b_1,B\left(X\right)=a_2X+b_2,R\left(X\right)=a_3X+b_3,$$

with $`a_i,b_iR.`$ Because of (34) and definition of the functions $`D,B,R`$ we obtain equations

$$\left(X^{}\right)^3=\left(c_3X^3+c_2X^2+c_1X+c_0\right)^2,$$
$`\left(35\right)`$
$$\left(T^{}\right)^3=A\left(c_3T^3+c_2T^2c_1X+c_0\right)^2,$$
$`\left(36\right)`$
where $`c_3,c_2,c_1`$ and $`c_0`$ are arbitrary constants, $`A=(2r)^{1/3}`$.

The solutions of (35) and (36) can be expressed in the terms of the Weierstrass function $`\mathrm{}`$ . Indeed, one can write (35) and (36) as

$$\left(X^{}\right)^3=\left(c_3\left(X\alpha _1\right)\left(X\alpha _2\right)\left(X\alpha _3\right)\right)^2,$$
$`\left(38\right)`$
$$\left(T^{}\right)^3=A\left(c_3\left(T+\alpha _1\right)\left(T+\alpha _2\right)\left(T+\alpha _3\right)\right)^2.$$

Replacing $`X\alpha _1`$ by $`1/Y`$ in (38) yields

$$\left(Y^{}\right)^3+B^2\left(Yb_1\right)^2\left(Yb_2\right)^2=0,$$

where $`B=c_3(\alpha _2\alpha _1)(\alpha _3\alpha _1),b_1=\frac{1}{\alpha _2\alpha _1},b_2=\frac{1}{\alpha _3\alpha _1}.`$ Introducing new function $`Z`$ such that $`Z^3=B(Yb_1)(Yb_2)`$, we obtain equation

$$\left(Z^{}\right)^2+\frac{4B}{9}Z^3+\frac{B^2\left(b_1b_2\right)^2}{9}=0.$$
$`\left(39\right)`$
The solutions of the last equation are expressed in the terms of the Weierstrass function $`\mathrm{}`$. Applying the above process to (36), we obtain equation such as (39).

We shall omit here other cases and discuss briefly in the next section two-dimensional equation.

4. Two-dimensional equation.

We consider here the fast diffusion equation

$$u_t=\mathrm{\Delta }\mathrm{ln}\left(u\right).$$
$`\left(40\right)`$
Some applications of this equation can be found in . Galaktionov  used invariant subspaces to find some solutions of (40). If we set $`u=1/v`$, we obtain

$$v_t=v^2\mathrm{\Delta }\mathrm{ln}\left(v\right).$$
$`\left(41\right)`$
It is easy to check that one of elementary solutions of (41) is travelling wave given by

$$v=1+c\mathrm{exp}\left(mx+ny\left(m^2+n^2\right)t\right),$$
$`\left(42\right)`$
where $`c,m`$ and $`n`$ are arbitrary constants. Obviously, this is invariant solution. On the other hand, this solution satisfies differential constraints

$$v_xmv+m=0,$$

$$v_ynv+n=0.$$

It is possible to find other differential constraints that are linear with respect $`x,y`$ and $`v`$.

It can be shown that the differential constraints

$$v_x+\frac{c_1c_0tan\left(t\right)}{c_0^2+c_1^2}\left(vxc_0yc_1+t\left(c_0^2+c_1^2\right)\right)=c_0,$$

$$v_y\frac{c_0+c_1tan\left(t\right)}{c_0^2+c_1^2}\left(vxc_0yc_1+t\left(c_0^2+c_1^2\right)\right)=c_1$$

are compatible with the equation (41). Here $`c_0`$ and $`c_1`$ are arbitrary constants. The solution of (41) corresponding to these constraints is

$$v=c_0x+c_1y\left(c_0^2+c_1^2\right)t+c_2\mathrm{cos}\left(t\right)\mathrm{exp}\left(m_1x+m_2y+m_3t\right).$$

Here $`c_0,c_1,c_2`$ are arbitrary constants and

$$m_1=\frac{c_0tc_1}{c_0^2+c_1^2},m_2=\frac{c_1t+c_0}{c_0^2+c_1^2},m_3=tan\left(t\right).$$

We can derive other explicit solutions using invariant subspaces  or linear differential constraints. For example, from  one may extract the following representation

$$v=s_0+s_1\mathrm{cos}\left(x\right)+s_2\mathrm{sin}\left(x\right)+s_3\mathrm{exp}\left(y\right)+s_4exp\left(y\right),$$

where functions $`s_i(t)`$ satisfy ordinary differential equations

$$s_0^{}+s_1^2+s_2^24s_3s_4=0,$$
$`\left(43\right)`$
$$s_1^{}+s_1s_0=0,$$

$$s_2^{}+s_2s_0=0,$$
$`\left(44\right)`$
$$s_3^{}s_3s_0=0,$$
$`\left(45\right)`$
$$s_4^{}s_4s_0=0.$$

Because of (44) and (45) we find

$$s_3^{}/s_2^{}+s_3/s_2=0.$$

This yields

$$s_2=c_2/s_3,c_2R.$$

By arguments similar to that used above we have

$$s_1=c_1/s_3,s_4=c_4s_3,c_1,c_4R.$$

Substituting this into (43) leads to

$$s_0^{}+\left(c_1^2+c_2^2\right)s_3^24c_4s_3^2=0.$$

From (45) we express the function $`s_0`$ and obtain

$$\left(\mathrm{ln}s_3\right)^{\prime \prime }=as_3^2bs_3^2,$$
$`\left(46\right)`$
with $`a=4c_4,b=c_1^2+c_2^2.`$

Setting $`a=b=1`$ one can derive two elementary solutions

$$s_3=tanh\left(t\right),s_3=tan\left(t\right)$$

In general, the solutions of (46) can be expressed in terms of elliptic functions. It is easy to obtain the correspondent function $`u`$.

The more complicated representation is

$$v=s_0+s_1\mathrm{cos}\left(2x\right)+s_2\mathrm{sin}\left(2x\right)+s_3\mathrm{exp}\left(2y\right)+s_4\mathrm{exp}\left(2y\right)+$$

$$s_5\mathrm{sin}\left(x\right)\mathrm{exp}\left(y\right)+s_6\mathrm{sin}\left(x\right)\mathrm{exp}\left(y\right)+s_7\mathrm{cos}\left(x\right)\mathrm{exp}\left(y\right)+s_8\mathrm{cos}\left(x\right)\mathrm{exp}\left(y\right),$$

where $`s_i`$ are functions which satisfy some ordinary differential equations. The special case of this representation was found in .

It is important to note that the equation (40) is invariant under infinite-dimensional algebra of symmetry . Some solutions of (40) were obtained by means of these symmetries in . We shall describe other method of using symmetry. It is convenient to apply the complex conjugate variables $`z=x+iy`$, $`\overline{z}=xiy`$. Thus, we can write the equation (40) as

$$u_t=\frac{1}{4}\frac{^2u}{z\overline{z}}.$$
$`\left(47\right)`$
It is easy to check that (47) is invariant under the transformation

$$z^{}=A\left(z\right),\overline{z}^{}=B\left(\overline{z}\right),u^{}=u/\left(A_zB_{\overline{z}}\right),$$

where $`A(z)`$ and $`B(\overline{z})`$ are arbitrary functions. In other words, if the function $`f(t,z,\overline{z})`$ is a solution of (47) then $`f(t,A(z),B(\overline{z}))A_zB_{\overline{z}}`$ also satisfies (47).

For example, if we set $`m=n=1`$ then from (42) we can construct the solution of the equation (40)

$$u=\frac{A_z\overline{A_{\overline{z}}}}{1+c\mathrm{exp}\left(A+\overline{A}2t\right)},$$

where $`A`$ is an arbitrary function of $`z`$ and $`\overline{A}`$ is the complex conjugate function.

5. Conclusions.

In sections 2 and 3 we have shown how the method of the linear determining equations can be applied to find explicit solutions to nonlinear diffusion equations. We have found exact solutions of these equations, using only the simplest solutions of the linear determining equations. It is interesting to find solutions of the linear determining equations depending on derivatives of higher orders. A.Shmidt  applied this method to another parabolic equations and some systems; application to the elliptic equation is discussed in .

In section 4 we have considered the two-dimensional equation. Applications of systems of the linear determining equations to multi-dimensional equations briefly discussed in . Using results of section 3 one can find the following representation

$$u=\left(a+be^{mt/2}\right)^2.$$
$`\left(48\right)`$
of solution of the equation

$$u_t=\mathrm{\Delta }\left(u^{1/2}\right)+mu+nu^{1/2},m,nR,$$

where the functions $`a(x,y)`$ and $`b(x,y)`$ must satisfy the system

$$\mathrm{\Delta }a=ma^2+na,$$

$$\mathrm{\Delta }b=mab+nb.$$

It is easy to show that the differential constraint

$$u_{tt}=u_t^2/2u+mu_t/2$$

leads to the representation (48). The interesting reductions of some diffusion equations in several independent variables can be found in . It is important to explain these reductions by means of differential constraints.

Acknowledgements

This work was supported by RFBR grant 01 - 01 - 00850, ME of Russia grant E00 - 10 - 57 and SB of RAS grant 1.


format: text

text:
# Nature of mechanical instabilities and their effect on kinetic friction

## Abstract

It has long been recognized that the key to understand kinetic friction force $`F_k`$ is the analysis of microscopic instabilities that lead to sudden irreversible ”pops” of certain degrees of freedom. In this Letter, the nature of such instabilities is characterized with an emphasis on boundary lubricants. It is shown that there are certain critical values of the parameters defining our model Hamiltonian, where the behavior of the instabilities changes qualitatively. Simultaneously, the functional dependence of $`F_k`$ on the sliding velocity $`v_0`$ changes. The relevant parameters studied here are dimensionality of the interface, degree of commensurability, first higher harmonic in the lubricant wall potential, and temperature. Molecular dynamics simulations are carried out to test whether the predictions made on the basis of the simple model also hold in less idealized circumstances.

When a solid slider is moved laterally with respect to a substrate, the kinetic friction $`F_k`$ is usually almost independent of the sliding velocity $`v_0`$  with leading corrections in the order of $`\mathrm{ln}v_0`$ . This so-called Coulomb friction differs from drag or Stokes friction that states a linear relation between $`F_k`$ and $`v_0`$ and that can be understood from equilibrium statistical mechanics as a special case of the fluctuation dissipation theorem : For instance, the drag force experienced by a Brownian particle in solution arises from the many collisions between the Brownian particle and the solvent molecules. In general, the condition for a linear friction force to occur is the (anharmonic) coupling of a central degree of freedom like the Brownian particle or a phonon to (infinitely) many other discrete degrees of freedom. In this Letter, the question will be addressed what the condition is for Coulomb friction to occur.

It has long been recognized that Coulomb friction must be related to instabilities that occur on a microscopic scale . When a slider is moved laterally with respect to the substrate, fast motion (pops) of certain degrees of freedom become unavoidable even if the slider s center of mass velocity is extremely small. The microscopic peak velocities in such pops are rather independent of $`v_0`$ and consequently the energy dissipated via a Stokes-type mechanism also becomes almost independent of $`v_0`$. The prototypical instability leading to Coulomb friction was suggested by Prandtl and Tomlinson . In their one-dimensional model, a surface atom in the slider is coupled with a spring of stiffness $`k`$ to its ideal lattice site which moves at constant velocity $`v_0`$. Interactions with the rigid substrate are modeled with a potential energy surface $`V`$ that is periodic in the substrate’s lattice constant plus some drag force linear in the atom’s veloctity $`\dot{x}`$. If $`k`$ is sufficiently small, pops become unavoidable and neglecting thermal fluctuations, $`F_k`$ remains finite in the limit of zero $`v_0`$. Instabilities and small-velocity kinetic friction in more complex model systems have been investigated since the Prandtl Tomlinson model was introduced. In particular, the $`F_k(v_0)`$ relationship of elastic manifolds sliding in ordered and disordered media has been studied extensively . The functional dependence is commonly found to be

$$F_k(v)=F_k(0)+cv^\beta .$$
(1)
While the elastic instabilities considered in those models are important in various contexts, this does not seem to be the case for the atomistic explanation of solid friction: Many detailed calculations and atomistic computer simulations reveal that in almost all cases, interbulk interactions are too weak to lead to instabilities at the atomic scale and as a consequence $`F_k(0)=0`$ . If interbulk interactions are very strong, irreversible processes like plastic deformation, material mixing, cold-welding, etc. usually occur and prevent the instabilities from being elastic.

It has been suggested that the presence of adsorbed particles, i.e. a boundary lubricant, confined between two surfaces is a more likely explanation for the commonly observed presence of solid friction . The main argument is that molecules that are only weakly bound to either surface can accommodate the surface corrugation of both walls simultaneously, which locks the walls together. This argument leads to static friction, which is the minimum force to initiate sliding between two solids. However, in order to explain kinetic friction, it is necessary to analyze the adiabatic solution $`x_{\mathrm{ad}}(t)`$ of the lubricant atoms. $`x_{\mathrm{ad}}(t)`$ is the athermal trajectory of an atom that always relaxes to the closest mechanical equilibrium position at every instant of time. The analysis of $`x_{\mathrm{ad}}(t)`$ has become common practice in the context of the motion of elastic manifolds in disordered media . In the impurity limit, interactions between lubricant atoms can be neglected and consequently $`x_{\mathrm{ad}}(t)`$ merely depends on the initial condition and the relative motion of top and bottom wall. In the following, $`x_{\mathrm{ad}}(t)`$ and its connection to kinetic friction will be examined. To the best of the author’s knowledge such an analysis has not yet been done, although the present model has already been used extensively to predict successfully various tribological phenomena .

The equation of motion for a lubricant atom in the boundary regime reads:

$`\ddot{x}`$ $`=`$ $`\gamma _b\dot{x}\gamma _t(\dot{x}v_0)+{\displaystyle \frac{1}{m}}\mathrm{\Gamma }(t)`$ (3)
$`{\displaystyle \frac{1}{m}}{\displaystyle \frac{}{x}}\left\{V_b(x)+V_t(xv_0t)\right\},`$
where $`x`$ denotes the atom’s position, $`m`$ is the atom’s mass, $`\gamma _t`$ and $`\gamma _b`$ parametrize the damping forces from the top and the bottom wall, and $`V_t`$ and $`V_b`$ denote the interaction of the confined atom with slider and the substrate. $`\mathrm{\Gamma }(t)`$ is a Langevin type stochastic random force defining temperature. For the lubricant wall interactions, various choices will be considered. Centrosymmetric choices can all be written in the form:

$$V_{t,b}=V_{t,b}^{(0)}\mathrm{cos}(x/b_{t,b})+V_{t,b}^{(1)}\mathrm{cos}(2x/b_{t,b})+\mathrm{}$$
(4)
where $`2\pi b_t`$ and $`2\pi b_b`$ are the period of the top and the bottom wall, respectively. The relevant physical units are defined through the choice $`V_t^{(0)}=1`$, $`m=1`$, $`b_t=1`$, and Boltzmann’s constant $`k_B=1`$. Furthermore, we will only consider slightly underdamped dynamics ($`\gamma _b=1x(t)x_{\mathrm{ad}}(t)`$) and restrict ourselves to the symmetric choice of $`V_0:=V_{t,b}^{(0)}`$. The free parameters are thus the degree of lattice mismatch ($`b_tb_b`$), the value of the first higher harmonic ($`V_1=V_{t,b}^{(1)}`$), and temperature $`T`$.

We start the discussion of a commensurate (com.) system ($`b_t=b_b`$) without higher harmonics in the absence of thermal fluctuations. In that case, the net time-dependent potential is simply given by $`V(t)=2V_0\mathrm{cos}\left(\frac{1}{2b}vt\right)\mathrm{cos}\left[\frac{1}{b}(xvt/2)\right]`$. Thus, for times $`\mathrm{cos}\left(\frac{1}{2b}vt\right)0`$, the atoms move at velocity $`v_0/2`$ as shown in Fig. 1. An infinitely small moment after a situation where this inequality was not satisfied, an atom will not be able to find a stable position in the immediate vicinity of the previous stable position. As the slider moves on with respect to the substrate, the atoms should therefore slide rapidly towards a new mechanical equilibrium. However, in order for this to happen, one needs a symmetry breaking element such as thermal fluctuations, $`\gamma _b\gamma _t`$, $`b_tb_b`$, or round-off errors in numerical calculations. Otherwise all stable trajectories simply are $`x=n\pi +v_0t/2`$ where $`n`$ is an integer number and the average (conservative) force of the lubricant atom acting on either wall will average to zero. For the athermal com. system, we chose to break symmetry by violating Galilei invariance and set $`\gamma _t=0`$. This choice is arbitrary but convenient, and is a better-controlled procedure than relying on round-off erros. We note also that for inc. walls, there is no physical reason for assuming the equality $`\gamma _t=\gamma _b`$.

Owing to broken symmetry ($`\gamma _t=0`$), the atom can now slide to the next minimum indicated by the gray lines in Fig. 1b. However, the peak velocities $`\dot{x}_p`$ in this process remain much smaller than during a pop within the Prandtl Tomlinson model. The reason is that due to the symmetry of $`V_t`$ and $`V_b`$, the pops occur between equivalent positions in the limit of arbitrarily small $`v_0`$. Since no lower bound for the dissipated energy can be given, the zero-velocity $`F_k`$ will be zero. At the same time, there is no upper bound for the ratio $`\dot{x}_p/v_0`$, so that $`F_k`$ cannot simply vanish linearly with $`v_0`$ but only with some power $`\beta `$ smaller than unity.

The nature of the instabilities changes qualitatively when the first higher harmonic is different from zero. Hence, in the sense of Morse theory , which contains Landau’s theory of phase transition as a special case, the com. system without higher harmonics can be considered a (multi-/tri-) critical point. If $`V_1<0`$, the adiabatic solution $`x(t)`$ becomes continuous as shown in Fig. 1, however, the time derivative $`\dot{x}(t)`$ diverges. If $`V_1>0`$, the adiabatic solution is discontinuous and the pops occur between inequivalent positions as shown in Fig. 1c. In analogy to phase transitions, the pops occurring for $`V_1>0`$ shall be called first-order instabilities, those for $`V_10`$ second-order instabilities. Only for first-order instabilities can one expect finite energy dissipation and thus finite $`F_k`$ as $`v_0`$ approaches zero. A numerical analysis shows that all three cases can be described with Eq. (1). The results are shown in Fig. 2 and the expected trend is confirmed: The more discontinuous the adiabatic solution, the larger the kinetic friction. It is important to emphasize that the exponent $`\beta `$ (as determined for sufficiently small $`v_0`$) only depends on the sign of $`V_1`$ but not on its precise value as long as $`V_1`$ is not too large.

For 1-dimensional, inc. surfaces, the basic picture is similar. If $`V_1`$ is larger than a (positive) critical value $`V_1^{}`$, whose precise value depends on the lattice mismatch, then pops between inequivalent positions are present and $`F_k`$ remains finite as $`v_0`$ tends to zero. For $`V_1<V_1^{}`$, however, the time derivative of the (continuous) adiabatic solution remains finite at all times. Hence the microscopic (peak) velocities $`v_p`$ scale linearly with $`v_0`$, which implies Stokes-type friction in that regime (1-d, inc.). Again, the exponent $`\beta `$ only depends on the sign of $`V_1V_1^{}`$. More detailed results including an analysis of significantly more comples systems will be presented in a separate paper .

There is however an important difference between com. and inc. systems, when the analysis is extended to two-dimensional interfaces: In 2-d inc. systems, first-order instabilities can now occur without higher harmonics. The reason is that in 2-d, atoms can circumnavigate the points of maximum longitudinal force. In 2-d com. surfaces, the behavior remains qualitatively similar as in 1-d, because the interference of $`V_t`$ and $`V_b`$ remains similar. Large-scale molecular dynamics simulations by He and Robbins support this argument even if the boundary lubricant is not any longer in the impurity regime. They find small $`F_k`$ and large $`F_s`$ between com. surfaces, while no such gap is observed for inc. systems . Also experimentally, signs for the effects of increased static friction of com. surfaces in the presence of a lubricant were reported .

This has potentially measurable implications for the transition from stick-slip motion to smooth sliding. In the stick-slip regime, friction is dominated by $`F_s`$, while in the smooth sliding regime, only $`F_k`$ is relevant. We have extended previous simulations  to support this point further. For this purpose the same model has been employed as that in Ref. . A schematic of the simulation is shown in Fig. 3 together with the average friction, as defined by the energy dissipated per slid distance. The verification of this prediction requires smooth surfaces, because rough surfaces automatically lead to inhomogeneous energy landscapes . The data shown in Fig. 3 was produced at a normal pressure of 0.4 GPa and a velocity of about 1 m/s using the same conversion of units as in Ref. . Note that the exceedingly small inertia of the slider as compared to experiment moves the transition from stick-slip to smooth sliding to large velocities.

We will now turn to the discussion of the effects of thermal fluctuations. He and Robbins found that velocity dependent corrections in a 2-dimensional, lubricated, inc. interface satisfy :

$$F_k(v_0)=F_k(v_{\mathrm{ref}})+𝒪(\mathrm{ln}(v_0/v_{\mathrm{ref}}))$$
(5)
over several orders of magnitude in $`v_0`$. This is different from the behavior found in the Prandtl Tomlinson model. Although similar temperature corrections have been suggested, for instance by Prandtl himself, more rigorous treatments yield corrections of order $`(T\mathrm{ln}v)^{2/3}`$ . While these corrections describe atomic force microscope experiments of nanoscale single-asperity contacts fairly accurately, it seems that they cannot provide an explanation of the usually observed $`\mathrm{ln}v_0`$ corrections in a straightforward manner.

Here it will be shown that simple logarithmic corrections are obeyed even in the impurity limit. Moreover, the crossover to linear response at extremely small sliding velocities will be included in the discussion. Fig. 3 shows the normalized friction force obtained at thermal energy $`T=0.07`$ for the com. $`V_1>0`$ model. Three regimes can be identified. At very small velocities, friction is linear in $`v_0`$ and one may associate this regime with the creep regime. At intermediate $`v_0`$, Eq. (5) is rather well satisfied. At ”large” velocities, thermal fluctuations become less relevant and the motion is close to that of the athermal system. The data obtained at different temperatures can be collapsed on a master curve. The collapse requires two dimensionless scaling functions $`r(T)`$ and $`s(T)`$ that both depend on temperature $`T`$ only. The collapse is done via $`F_k(v,T)=r(T^{})F_k(v^{},T^{})/r(T`$ with $`s(T)\mathrm{ln}(v^{}(T^{})/v_1)=s(T^{})\mathrm{ln}(v(T)/v_1)`$ where $`v_1`$ is a constant, $`r(T)`$ is almost constant \[$`r(0.02)=0.98,r(0.2)=1.2`$\] and $`s(T)k_bT/\mathrm{\Delta }E`$ where $`\mathrm{\Delta }E`$ can be interpreted as an effective (free) energy barrier. Qualitatively similar crossover from the linear response regime to the activated regime are observed in many other systems such as single particles in a static periodic potential , driven thermal elastic manifolds , and shear-thinning fluids .

In conclusion, this Letter provides a classification scheme for instabilities that can occur when two solids are in relative sliding motion. First-order instabilities are defined as pops of atoms (or other degrees of freedom) between inequivalent positions. They lead to kinetic friction that remains finite when the sliding velocity $`v_0`$ goes to zero provided the system is athermal. The exponent $`\beta `$ of the velocity corrections $`v^\beta `$ depends on the details of the model, however, $`\beta `$ only changes its value at certain critical points in the parameter space defining interactions and geometry. At finite temperature, corrections in the order of $`\mathrm{ln}v_0`$ apply. Second-order instabilities are defined as pops between equivalent positions. They lead to a sub-linear power law $`F_kv^\beta `$. If, however, the adiabatic solution $`x_{\mathrm{ad}}`$ of the boundary lubricant moves with finite velocity at all times, then simple Stokes friction follows. Thus an important result of this analysis is the identification of critical points where the kinetic friction law changes qualitatively as a parameter describing the interactions or the geometry is varied.

While the present study is primarily concerned with dilute boundary lubricants, the concept itself seems to be more general. For instance $`x_{\mathrm{ad}}`$ can be a collective order parameter that fluctuates back and forth between two values. Such quasi-periodic phase transitions have been observed in computer simulations of Ni asperities moving over a Cu substrate . Of course, the situation in those simulations was more complex, because wear occurred as a side effect of the motion. In the other extreme, $`x_{\mathrm{ad}}`$ may merely denote the position of an electronic orbital.

I thank K. Binder for useful discussions. Support from the BMBF through Grant 03N6015 and from the Materialwissenschaftliche Forschungszentrum Rheinland-Pfalz is acknowledged.


text:
# Classical gluodynamics in curved space–time and the soft Pomeron

## 1 Introduction

At the classical level, QCD is a scale invariant theory. Quantum fluctuations, however, give rise to scale anomaly  which is associated with dimensionful scale $`\mathrm{\Lambda }0.2`$ GeV. The perturbative expansion is possible if the typical distance $`Q^1`$ at which the interaction occurs is much smaller than $`\mathrm{\Lambda }^1`$, so the coupling is weak : $`\alpha _s(Q)1`$. However, this does not mean that at short distances non-perturbative contributions can always be neglected. Indeed, the strength of the semi-classical vacuum fluctuations of the gluon field is inversely proportional to the coupling constant and thus can be large at short distances, when $`1/\alpha _s(Q)1`$. The goal of the present work is to derive an effective lagrangian for gluodynamics which makes it possible to systematically treat such contributions to the scattering amplitudes.

Gluodynamics is defined by the lagrangian

$$_{\mathrm{gluon}}=\frac{1}{4}F_{\mu \nu }^aF_{\mu \nu }^a.$$
(1)
If treated classically, this theory is invariant under the scale transformations $`x\lambda x`$. Transformation of scale is associated with the scale current $`s_\mu `$; its divergence is equal to the trace of the energy–momentum tensor:

$$^\mu s_\mu =\theta _\mu ^\mu (x).$$
(2)
Without quantum effects, $`\theta _\mu ^\mu (x)=0`$, and the theory is scale invariant. Quantum effects break scale invariance . However, the broken symmetry still manifests itself in the following set of low energy theorems for different Green functions involving operator $`\theta _\mu ^\mu (x)`$ which can be proven by using renormalization group arguments :

$$i^n𝑑x_1\mathrm{}𝑑x_n0|T\{\theta _{\mu _1}^{\mu _1}(x_1),\mathrm{},\theta _{\mu _n}^{\mu _n}(x_n)\}|0_{\mathrm{connected}}=\theta _\mu ^\mu (x)_{\mathrm{vac}}(4)^n$$
(3)
In Ref.  it was shown that these low–energy theorems entirely determine the form of the effective low energy lagrangian for gluodynamics, which describes semi-classical vacuum fluctuations of gluon field at large distances. It is formulated in terms of an effective scalar dilaton field:

$$_{\mathrm{dilaton}}=\frac{|ϵ_\mathrm{v}|}{m^2}\frac{1}{2}e^{\chi /2}(_\mu \chi )^2+|ϵ_\mathrm{v}|e^\chi (1\chi ),$$
(4)
where the real scalar field $`\chi `$ of mass $`m`$ represents the dilaton, and $`|ϵ_\mathrm{v}|`$ is the vacuum energy density.

The idea underlying the derivation can be formulated in the way which illustrates the analogy between gluodynamics and general relativity: QCD coupled to the conformally flat gravity is scale and conformally invariant in any number of dimensions, and thus such theory is anomaly–free. However, the vacuum energy density of flat–space gluodynamics $`|ϵ_\mathrm{v}|`$ (which plays a role of cosmological constant in general relativity) manifestly breaks scale and conformal symmetry of gluodynamics in conformally flat space-time. The low energy effective lagrangian (4) thus can be derived from the Einstein-Hilbert lagrangian for the described theory. Note that in the limit $`|ϵ_\mathrm{v}|0`$ lagrangian (4) is scale and conformally invariant, since another dimensionful parameter $`m`$ just sets the normalization of the dilaton field $`\chi `$. At small dilaton momenta $`pm`$ the lagrangian (4) satisfies the low energy theorems (3). It is remarkable that the effective lagrangian involves only one scalar field. The reason is that the gravity is excited by the trace of energy momentum-tensor, so there is only one independent Einstein equation.

In the present work we are going to extend the range of validity of the effective lagrangian (4) to shorter distances by including the interactions of dilaton field with gluons. Generally, it can be obtained from the lagrangian of gluodynamics by integrating out soft gluon modes up to the scale $`m`$ which thus becomes the dimensionful scale of the theory. This implies that the scale $`\mathrm{\Lambda }`$ and the dilaton mass $`m`$ are intimately connected. Actually, as we will show, it appears more convenient to construct the generalized effective lagrangian using Einstein equations once again. In the next section we employ this method to obtain the lagrangian (5). One can easily check that it is non-renormalizable. Therefore, we need to specify the additional dimensionful scale – the ultra-violet cutoff $`M_0`$. Obviously, $`M_0`$ is a function of $`|ϵ_\mathrm{v}|`$ and $`m`$. It can be calculated using the requirement of vacuum stability, i.e. all radiative corrections must cancel out to give correct value of vacuum energy density. Its numerical value turns out to be quite large, $`M_03`$ GeV. This is in agreement with Ref.  where it has been argued that the typical scale of the vacuum gluon density due to semi-classical fluctuations can indeed be that large. In sec. 2 we discuss the effective lagrangian (5) in detail and calculate the scale $`M_0`$.

In sec. 3 we apply the effective lagrangian (5) to calculate a contribution of the non-perturbative semi-classical fluctuations of the QCD vacuum to the pomeron intercept. Experimental data on different reactions at high energies and small momentum transfer support existence of the pole in the scattering amplitude responsible for high energy behavior of the total cross sections . This leads to the power-like behavior of the total cross section $`\sigma _{\mathrm{tot}}=\varsigma s^{\alpha _P1}`$, where $`\varsigma `$ depends upon particular reaction and $`\alpha _P`$ is the universal parameter called the intercept. After decades of hard work, understanding the nature of the object exchanged in $`t`$-channel (pomeron) is still a challenging problem. It turned out that the perturbative QCD is far from giving a clear answer even if the scales inherent to colliding particles are hard. Partly, this is due to diffusion of the transverse momenta in the partonic cascade towards the non-perturbative scales. This implies that the main contribution to the pomeron stems, perhaps, from the non-perturbative sector of QCD. Assuming that some non-perturbative fluctuations of vacuum dominate the scattering amplitude at high energies one has to explain the power-like behavior of the total cross section $`\sigma _{\mathrm{tot}}s^{\alpha _P1}`$ and to obtain the correct numerical value of $`\alpha _P`$. A few authors have reported on this problem. Contribution of instantons to the scattering amplitude at high energies was investigated in Ref. . In Ref.  it has been argued that fluctuations of strings in the framework of AdS/CFT theory can yield the desired result. In Refs.  the spectrum of glueballs in QCD was discussed as well as its relation to the structure of the pomeron. Possible contributions of extra compact dimensions was discussed in Ref. . A more phenomenological approach has been developed in framework of the BFKL gluon emission with additional assumption for such an emission in the non-perturbative QCD domain . In Ref.  two of us have used the fact that the non-perturbative QCD vacuum is dominated by the semi-classical fluctuations of the gluon field which is inversely proportional to the coupling $`F1/g`$. The importance of such fluctuations can be understood in the following way. Consider a gluon ladder in each rung of which two gluons in a singlet state are produced in $`s`$-channel. Each rung gives contribution of the order $`g^2F^2𝒪(1)`$ since the two-gluon singlet state gives contribution proportional to the trace of energy-momentum tensor. Such ladder diagrams can be summed up yielding a reasonable value of the pomeron intercept. In the present paper we suggest a method for systematic calculation of processes involving such non-perturbative modes of gluodynamics. In particular, we will calculate the total cross section at high energies, reproduce its power-like behavior and numerical value of the intercept. We also calculate the final state multiplicity in this approach.

We conclude the paper by presenting some speculative ideas on the possible link between QCD and gravity. Recall that the effective lagrangian (5) can be derived from the General Relativity. Thus the gravitational constant can be expressed in terms of the dilaton mass and the vacuum energy density, Eq. (11). Its numerical value $`G=0.2\mathrm{GeV}^2`$, of course, has nothing to do with the Planck gravitational constant. Nevertheless, it may be possible to interpret this result as a signature of the strong gravitational interactions at scale much lower than the Planck one. If we accept this interpretation, then it turns out that the gravitational effects are essential already at the GeV scale, in order to be responsible for the QCD vacuum structure. Several authors have suggested mechanisms by which gravitational interactions can become strong at scales much larger than the Planck one . Recently it has been suggested that quantum gravity can become strong at the TeV scale . In sec. 4 we briefly discuss those ideas.

## 2 Dilaton effective lagrangian

The dilaton effective lagrangian in gluodynamics describing the propagation and interactions of gluon and dilaton fields can be constructed as a generalization of Eqs. (4),(1). It reads

$$=\frac{|ϵ_\mathrm{v}|}{m^2}\frac{1}{2}e^{\chi /2}(_\mu \chi )^2+|ϵ_\mathrm{v}|e^\chi (1\chi )e^\chi (1\chi )\frac{1}{4}F_{\mu \nu }^aF^{a\mu \nu }.$$
(5)
Eq. (5) can be obtained utilizing Legendre transformation of the lagrangian describing coupling of gluodynamics to gravity

$$^{}=\left(\frac{1}{8\pi G}R\frac{1}{4}F_{\mu \nu }^aF^{a\mu \nu }+|ϵ_\mathrm{v}|\right)\sqrt{detg_{\mu \nu }}.$$
(6)
The first term in (6) corresponds to the gravitational lagrangian $`_G`$. It is expressed through the Ricci scalar $`R`$ which can be easily calculated for the conformally trivial metric $`g_{\mu \nu }=\eta _{\mu \nu }e^{h(x)}`$

$$_G=\frac{1}{8\pi G}R\sqrt{detg_{\mu \nu }}=\frac{1}{8\pi G}3e^h\left(_\mu ^2h+\frac{1}{2}(_\mu h)^2\right)=\frac{1}{8\pi G}\frac{3}{2}e^h(_\mu h)^2,$$
(7)
where we have used integration by parts to derive the last equation and $`\eta _{\mu \nu }=\mathrm{diag}\{1,1,1,1\}`$. In (7) $`G`$ is the gravitational constant. The gravitational field $`h`$ couples to the gluodynamics by means of the second term

$$\frac{1}{4}F_{\mu \nu }^aF^{a\mu \nu }e^{2h}.$$

These two terms are invariant under scale and conformal transformations in a space-time of arbitrary dimension. However, the scale symmetry in QCD is broken down by non-vanishing gluon condensate $`F_{\mu \nu }^aF_{\mu \nu }^a`$$`>0`$. This leads to the scale anomaly which manifests itself as an anomaly of the QCD energy-momentum tensor ($`\theta ^{\mu \nu }`$) trace

$$\theta _\mu ^\mu =\frac{\beta (g)}{2g}F_{\mu \nu }^aF^{a\mu \nu },$$
(8)
where the $`\beta `$-function in the leading order in the strong coupling $`g`$ reads

$$\beta (g)=\frac{11g^3}{(4\pi )^2}.$$
(9)
Its vacuum expectation value is

$$\theta _\mu ^\mu =4|ϵ_\mathrm{v}|.$$
(10)
In order for the energy-momentum tensor, calculated from the lagrangian (6) to satisfy (10), the third term in (6) which explicitly breaks the scale invariance has been introduced. This term can also be thought of as the cosmological constant. In the absense of the gravitational part (7) of the lagrangian (5) the condition (10) is satisfied trivially at the tree level. Performing the Legendre transformation of (6) one arrives at (5) . The mass of the dilaton can be expressed through the vacuum energy density and the gravitational constant as follows

$$m^2=\frac{64\pi }{3}|ϵ_\mathrm{v}|G.$$
(11)
Now we require the energy momentum tensor calculated from the lagrangian (5) to satisfy the vacuum normalization condition (10). This means cancelation of the dilaton and gluon contributions to the vacuum energy density. We can formulate (10) in terms of the vacuum expectation values. To this end let us calculate the trace of the energy-momentum tensor

$$\theta _\mu ^\mu =\eta ^{\mu \nu }\left(2\frac{}{\eta ^{\mu \nu }}\eta ^{\mu \nu }\right)+\frac{8|ϵ_\mathrm{v}|}{m^2}_\mu ^2e^{\chi /2},$$
(12)
where the last term in the right hand side is the total derivative . Equation of motion for the dilaton field is

$$\frac{|ϵ_\mathrm{v}|}{m^2}_\mu \left(e^{\chi /2}_\mu \chi \right)\frac{|ϵ_\mathrm{v}|}{4m^2}e^{\chi /2}\left(_\mu \chi \right)^2+\chi e^\chi |ϵ_\mathrm{v}|\chi e^\chi \frac{1}{4}F_{\mu \nu }^aF^{a\mu \nu }=0.$$
(13)
Using (13) in (12) we arrive at

$$\theta _\mu ^\mu =4|ϵ_\mathrm{v}|e^\chi +\chi e^\chi F_{\mu \nu }^aF^{a\mu \nu }.$$
(14)
Hence, (10) is equivalent to the following relation

$$\chi e^\chi F_{\mu \nu }^aF^{a\mu \nu }=4|ϵ_\mathrm{v}|e^\chi 1.$$
(15)
It is clear from equations (14) and (13) that in the limit $`|ϵ_\mathrm{v}|0`$ the classical equation $`\theta _\mu ^\mu =0`$ is recovered. However, quantum corrections can give non–vanishing contributions to $`\theta _\mu ^\mu `$ in the same limit by introducing dimensonal cutoff. Eq. (15) assures that this does not happen in our theory.

Now, we proceed to the derivation of Feynman rules. In addition to usual Fyenman diagrams of gluodynamics there are those corresponding to the dilaton propagator, dilaton self-interactions and dilaton-gluon interactions. In this work we are interested only in contributions of the order of $`g^0`$. They are shown in Fig. 1. Other vertices, which contain three and four gluons coupled to the dilaton are supressed by factors $`g^1`$ and $`g^2`$ respectively. In our approach the strong coupling constant is necessarily small since all gluon modes which enter the lagrangian (5) are hard by construction. All soft modes, which carry momenta $`p<m`$, have been integrated out. The dilaton Feynman propagator can be read from (5):

$$D_d(p)=\frac{m^2}{|ϵ_\mathrm{v}|}\frac{i}{p^2m^2i\epsilon }.$$
(16)
To completely specify our effective theory we still have to calculate the ultra-violet cutoff $`M_0`$. Contributions to the vacuum energy density involve loops and therefore they are divergent in the ultra-violet regime of the theory. Let us introduce the cut-off $`M_0`$ which controls divergences of vacuum excitations contributing to the vacuum expectation value of the trace of the energy-momentum tensor. The same parameter controls divergences in the dilaton phase space. Our requirement that the normalization (10), or equivalently (15), holds means that we impose a restriction on the possible values of $`M_0`$. Consider a contribution stemming from pure dilaton vacuum fluctuations, Fig. 2(a).

In the next section we will argue that for the processes we are interested in here the phase space is available for the emission of only two dilatons. In this two-dilaton approximation

$$\theta _\mu ^\mu _{\mathrm{vac}}=4|ϵ_\mathrm{v}|e^\chi _{\mathrm{vac}}4|ϵ_\mathrm{v}|\left(1+\frac{1}{2}\chi ^2_{\mathrm{vac}}\right).$$
(17)
Hence, the contribution to the vacuum energy-density is

$$2|ϵ_\mathrm{v}|\chi ^2_{\mathrm{vac}}=2|ϵ_\mathrm{v}|\frac{1}{\pi }_0^{M_0^2}\frac{dM^2}{M^2}\frac{m^2}{|ϵ_\mathrm{v}|}\frac{d^3k}{(2\pi )^32\omega }=\frac{1}{4\pi ^3}m^2M_0^2.$$
(18)
Another contribution of the same order in the strong coupling, $`𝒪(g^0)`$, comes from the first diagram shown in Fig. 2(b). We have

$$(N_c^21)\frac{1}{\pi }_0^{M_0^2}\frac{dM^2}{M^2}\frac{1}{2!}\frac{d^3q_1}{(2\pi )^32\omega _1}\frac{d^3q_2}{(2\pi )^32\omega _2}\mathrm{\Gamma }_2\left(q_1^\nu q_2^\mu (q_1q_2)\eta ^{\mu \nu }\right)^2\theta (M\omega _1\omega _2)$$

$$=\frac{M_0^8}{28048(2\pi )^6}(N_c^21)\left(\frac{m^2}{|ϵ_\mathrm{v}|}\right)^2,$$
(19)
where $`\mathrm{\Gamma }_2`$ is a two-dilaton phase-space factor, see (28). Contributions (18) and (19) cancel out when

$$M_0^2=8\pi \mathrm{\hspace{0.17em}2}^{2/3}\mathrm{\hspace{0.17em}105}^{1/3}(N_c^21)^{1/3}\left(\frac{|ϵ_\mathrm{v}|}{m}\right)^{2/3}.$$
(20)
The pure gluon vacuum diagrams are proportional to the strong coupling constant $`\alpha _s`$ as can be easily seen from (8) and (9). Inasmuch as the strong coupling constant is small, those contributions to the vacuum energy density can be neglected.

## 3 Derivation of the soft pomeron

The leading contribution to the total cross section at high energies stems from the ladder diagrams . Such a diagram is shown in Fig. 3. Note that the energy dependence of the Born amplitude is determined by the spin of particles exchanged in $`t`$-channel. Thus the leading contribution at high energies comes from exchange of vector particles – gluons . These gluons prefer to form ladders, since each rung gets enhanced by large logarithm of energy $`\mathrm{ln}s1`$. By virtue of the optical theorem, the total cross section $`\sigma _{tot}`$ can be calculated as a sum of all such diagrams

$$\sigma _{tot}=\underset{r0}{}\underset{n2}{}\sigma _r^n=\frac{2}{s}\underset{r0}{}\underset{n2}{}\mathrm{Im}_{\mathrm{el}}(t=0)_r^n,$$
(21)
i.e. we sum over all possible numbers $`n2`$ of dilatons produced in $`s`$-channel in each rung for a given number $`r`$ of rungs and then sum over all $`r`$.

The most convenient way to calculate the amplitudes associated with ladder diagrams is to employ the Weizsäcker-Williams approximation. This means replacing almost longitudinally polarized Coulomb gluons exchanged in $`t`$-channel by transversely polarized (real) ones. In this approximation the gluon propagator reads :

$$D_g(q^2)\frac{2}{s}\frac{q_{}^\mu q_{}^\nu }{\alpha _q\beta _qq_{}^2},$$
(22)
where $`\alpha _q`$ and $`\beta _q`$ are Sudakov variables defined by

$$q=\alpha _qp_1+\beta _qp_2+q_{}.$$
(23)
Weizsäcker-Williams approximation is valid as long as the following conditions hold

$`1\alpha _{q_1}\mathrm{}\alpha _{q_i}\mathrm{}\alpha _{q_{r+1}},`$ (24)
$`1\beta _{q_{r+1}}\mathrm{}\beta _{q_i}\mathrm{}\beta _{q_1},`$ (25)
$`|q_i^2|s.`$ (26)
In the Born approximation (two gluon exchange) there are no produced dilatons. Contribution of a one-rung diagram to the total cross section reads

$$\sigma _{r=1}^n=\frac{2}{s}\frac{d^4q_1}{(2\pi )^4}\frac{d^4q_2}{(2\pi )^4}\left|\frac{2g^2(q_1q_2)_{}\mathrm{Tr}(t^at^b)}{\beta _{q_1}\alpha _{q_2}s}\right|^2\frac{1}{4N_c^2}(2\pi )^2\delta ((p_1q_1)^2)\delta ((p_2+q_2)^2)\mathrm{\Gamma }_n(M),$$
(27)
where $`\mathrm{\Gamma }_n(M)`$ represents contribution of $`n`$ dilatons to the phase space and $`M^2=(q_1q_2)^2`$ is the four-momentum squared available for production in the dilaton vertex. The phase space $`\mathrm{\Gamma }(M)`$ has been evaluated in Ref. :

$$\mathrm{\Gamma }_n(M)=\frac{1}{n!}\left(\frac{m^2}{|ϵ_\mathrm{v}|}\right)^n(n1)^2(2\pi )^4\delta ^4\left(\underset{i=1}{\overset{n}{}}k_ik\right)\underset{i=1}{\overset{n}{}}\frac{d^3k_i}{(2\pi )^32\omega _i}$$

$$=\left(\frac{m^2}{|ϵ_\mathrm{v}|}\right)^n(2\pi )^{43n}(n1)^2\frac{(\pi /2)^{n1}(M^2)^{n2}}{n!(n1)!(n2)!},$$
(28)
where $`(n1)^2`$ comes from the vertex shown in Fig. 1(a). Using (27) and (28) we derive the contribution of one-ladder diagram

$$\underset{n2}{}\sigma _{r=1}^n=\sigma _B\frac{\mathrm{ln}s}{(2\pi )^6}\pi ^2M_0^2_0^{M_0^2}𝑑M^2\underset{n=2}{\overset{\mathrm{}}{}}\mathrm{\Gamma }_n(M)=\sigma _B\mathrm{\Delta }_P\mathrm{ln}s,$$
(29)
where $`\sigma _B\alpha _s`$ is a Born amplitude,

$$\mathrm{\Delta }_P=\frac{a^2}{4\pi }_0F_2(;1,3;a),$$
(30)
$`{}_{0}{}^{}F_{2}^{}`$ is a generalized hypergeometric function, and the relevant parameter of our theory for high energy scattering is

$$a=\frac{m^2M_0^2}{16\pi ^2|ϵ_\mathrm{v}|}.$$
(31)
Expansion of function $`{}_{0}{}^{}F_{2}^{}`$ in powers of parameter $`a`$

$$_0F_2(;1,3;a)=1+\frac{a}{3}+\frac{a^2}{48}+\frac{a^3}{2160}+𝒪(a^4)$$
(32)
is equivalent to expansion in the number of emitted dilatons. We will argue below that the series (32) are rapidly convergent and actually, there is not enough phase space for emission of large number $`n`$ of dilatons since $`\mathrm{\Gamma }_n(M)`$ is a rapidly decreasing function of $`n`$. Thus, we can safely neglect emission of more than two dilatons.

Now, as we know the contribution of one-rung diagram, calculation of the total cross section is straightforward. Indeed, the contributions of different rungs factorize out, leading to the Regge-like behavior of the total cross section (21)

$$\sigma _{tot}=\sigma _B\underset{r=0}{\overset{\mathrm{}}{}}\frac{1}{r!}\mathrm{ln}^rs\mathrm{\Delta }_P^r=\sigma _Bs^{\mathrm{\Delta }_P}.$$
(33)
Corrections to the pomeron intercept coming from perturbative QCD are proportional to the strong coupling constant $`g^2`$ and neglected in this approach. Loop corrections to the dilaton propagator give small contributions of the order of

$$\mathrm{\Delta }D_d\frac{am^2}{M_0^2}1.$$

Finally, upon substitution of (31) and (20) into (30) we obtain the following result for the soft pomeron intercept

$$\alpha _P=1+\mathrm{\Delta }_P=1+\frac{(105/2)^{2/3}}{4\pi ^3}(N_c^21)^{2/3}\left(\frac{m^4}{|ϵ_\mathrm{v}|}\right)^{2/3}.$$
(34)
To estimate $`|ϵ_\mathrm{v}|`$ in pure gluodynamics we use the fact that it is related to the vacuum energy density in QCD $`|ϵ_\mathrm{v}|_{\mathrm{QCD}}`$ by

$$|ϵ_\mathrm{v}|\frac{11}{9}\mathrm{\hspace{0.17em}3}|ϵ_\mathrm{v}|_{\mathrm{QCD}}.$$
(35)
Indeed, as was argued in Ref. ,

$$\alpha _sF_{\mu \nu }^aF^{a\mu \nu }_{\mathrm{glue}}\mathrm{\hspace{0.17em}3}\alpha _sF_{\mu \nu }^aF^{a\mu \nu }_{\mathrm{QCD}}.$$
(36)
Also the beta function (9) in QCD scales with number of col-ours and flavors as $`\beta _{\mathrm{QCD}}(11N_c2N_f)/3`$ whereas in pure gluodynamics $`\beta _{\mathrm{glue}}11N_c/3`$. Since by (8) $`|ϵ_\mathrm{v}|\beta `$ the vacuum energy density in pure gluodynamics gets enhanced by the factor $`11/9`$. Together with (36) this yields (35). Sum rules analysis of  makes it possible to estimate the QCD vacuum energy density $`|ϵ_\mathrm{v}|_{\mathrm{QCD}}=(0.24\mathrm{GeV})^4`$. Owing to (35) we obtain the following estimate:

$$|ϵ_\mathrm{v}|=0.012\mathrm{GeV}^4.$$
(37)
Another parameter of our effective theory is the dilaton mass. Note, that due to the isospin conservation each dilaton can produce only even number of pions in the final state. Let us assume that in full QCD with light quarks the dilaton mixes with the broad physical scalar $`\sigma `$–resonance  at

$$m=0.4÷0.5\mathrm{GeV}.$$
(38)
Note that the spectral density of the dilaton excitations can be spread over several scalar resonances of different mass (see, e.g., ). However, for our purposes it suffices to keep only the lightest resonance, which gives the dominant contribution to the intercept because of the phase space constraints.

Let us now proceed to the numerical estimates. In view of (37) and (38) Eq. (20) implies that $`M_0=3.0÷2.8\mathrm{GeV}`$. Substituting these numbers to (34) we obtain for the pomeron intercept $`\alpha _P=1.047÷1.085`$ which is in reasonable agreement with the phenomenology of high energy scattering . Using these numbers in (32), we see that the emission of three and more dilatons contributes $`1{}_{0}{}^{}F_{2}^{}(;1,3;0)/_0F_2(;1,3;a)`$=$`21÷27\%`$ to the final result; this can be used as an estimate of the accuracy of our two-dilaton approximation.

We can calculate the inclusive spectrum of dilaton production by observing that Eq. (28) is nothing but the number distribution of the dilatons produced per unit of rapidity. Indeed, each rung of the ladder diagram corresponds to one unit in rapidity (multiplied by $`\mathrm{\Delta }_P`$) in the leading logarithmic approximation. Define the probability $`P_n(a)`$ to find $`n`$ dilatons in the final state

$$P_n(a)=\frac{_0^{M_0^2}𝑑M^2\mathrm{\Gamma }_n(M)}{_{n=2}^{\mathrm{}}_0^{M_0^2}𝑑M^2\mathrm{\Gamma }_n(M)}.$$
(39)
Then the average number of dilatons produced per unit of rapidity is given by

$$\frac{dN_\chi }{d(\mathrm{\Delta }_Py)}=\underset{n=2}{\overset{\mathrm{}}{}}nP_n(a)=\frac{2_0F_2(;1,2;a)}{{}_{0}{}^{}F_{2}^{}(;1,3;a)}.$$
(40)
From (31) we get an estimate of the dilaton phase space parameter $`a=0.768÷1.033`$. Substituting this value into (40) one arrives at $`dN_\chi /dy=0.1÷0.2`$. Since we have identified the dilaton with the physical $`\sigma `$-resonance, we have to assume that it decays mostly in two pions. Charged pions give approximately 2/3 of the total number of pions. Therefore, the average multiplicity of charged pions per unit of rapidity is estimated as

$$\frac{dN_\pi ^{\mathrm{ch}}}{dy}=\frac{2}{3}\mathrm{\hspace{0.17em}2}\frac{dN_\chi }{dy}=\mathrm{\hspace{0.17em}0.14}÷0.26.$$
(41)
This agrees with the result obtained in: the soft pomeron is responsible for less then 10% of the observed multiplicity. In Fig. 4 we show the number distribution (39) of produced charged pions.

## 4 Summary and discussion

Summarizing, we have started with the derivation of the QCD effective lagrangian (5) which incorporates scale anomaly and represents propagation and interaction of gluons and dilatons. The dilaton field saturating the sum rules (3) describes the QCD vacuum structure. The pure dilaton part of lagrangian (5) (first and second terms) is a low energy effective gluodynamics obtained by integrating out low energy modes with momenta $`q`$ such that $`qm0.4÷0.5`$ GeV. The third term of (5) represents pure gluodynamics and the interaction of gluon field with the dilaton one. It is important as far as momentum transfer does not exceed $`M_0`$. It has been already pointed out in Ref.  that the scale $`M_0`$, determined by the semi-classical fluctuations of the gluon fields in vacuum is quite large, $`M_0\mathrm{\Lambda }`$.

Normalizing the vacuum energy density by (10) we obtained an expression (20) for the ultra-violet cutoff $`M_0`$ of the theory in terms of vacuum energy density $`|ϵ_\mathrm{v}|`$ and the dilaton mass $`m`$. Divergent integrals that we have encountered in Eqs. (18),(19),(29) should not surprise the reader since the dilaton acts as the gravitational field. Indeed, quantum gravity is a non-renormalizable theory since the coupling of gravity to matter fields $`G`$ has canonical dimension (mass)<sup>-2</sup>. In technical terms, each emitted dilaton contributes a factor of $`m^2/|ϵ_\mathrm{v}|`$ coming from the normalization of the dilaton propagator which must be compensated for by $`M_0^2`$ to maintain the scattering amplitude dimensionless.

In the framework of this effective theory we addressed the problem of soft pomeron. We calculated the high energy scattering amplitude in the leading logarithmic approximation by summing the ladder diagrams, see Fig. 3. Generally, any number of dilatons can be emitted in the $`t`$-channel in each ladder rung, but actually most of the phase space is occupied by only two of them. Owing to the isospin conservation the dilaton in the final state decays into an even number of pions. By identifying the dilaton with the physical $`\sigma `$ resonance, we conclude that it should decay mostly in two pions. We estimate the pomeron intercept $`\alpha _P=1.047÷1.085`$. This value agrees well with the analysis of the data on high energy scattering .

We also estimate the final state multiplicity per unit of rapidity $`\frac{dN_\chi }{dy}0.2`$. The low particle multiplicity resulting from the soft pomeron exchange has been noted already in . It implies that the bulk of particles is produced by another mechanism. Such a mechanism must involve multiple production of particles per unit of rapidity and perhaps closely related to high parton density regime of QCD .

The reader familiar with the previous work of two of us, Ref. , might be confused by the fact that there the opposite to (34) dependence of $`\mathrm{\Delta }`$ on the vacuum energy density $`|ϵ_\mathrm{v}|`$ and the scalar particle (glueball) mass $`m_R`$ has been found. The reason is simple: a glueball mediates a short-range strong interaction, so the cross section for its production is proportional to $`m_R^4`$. On the other hand, gravitation is a long-range interaction, the strength of which is proportional to masses of interacting particles (dilatons). In other words, the effective lagrangian (5) is constructed to obey scale invariance on the classical level. This requirement translates into the absence of the coupling of the dilaton to two gluons, because a term linear in the dilaton field $`\chi `$ in (5) would induce a non-zero vacuum expectation value $`\chi `$ which would violate the scale invariance already on the classical level. (Even though the field $`\chi `$ is dimensionless, its couplings are not, and so its v.e.v. would violate scale invariance of the vacuum state). This is why the coupling of the scalar glueball to two gluons, used in , does not appear in the approach of this paper, which explains the difference in the dependence of the final results on the scalar particle mass.

The reader may also wonder about the meaning of coupling QCD to gravity. One way of interpreting this is to say that coupling to gravity is a purely mathematical trick which makes it possible to take all symmetries into account in the most elegant way .

However, there might be another interpretation which goes beyond the standard model. We treated the gravitational constant $`G`$ as a free parameter. We can find its numerical value using (11): $`G=0.2`$ GeV<sup>-2</sup>. One might thus assume that gravitational interactions become strong already at the scale $`G^{1/2}`$2.2 GeV. In this case cutting divergent integrals associated with dilaton fields at the scale of $`M_0`$ means that the discussed theory (5) is an effective one. It implies that the full renormalizable theory must be recovered at some higher energies. In this sense our scale $`M_0`$ can be derived from some fundamental scale $`M^{}`$. In recent years we have learnt that the existence of strong gravity regime at experimentally accessible energies can be realized in theoretical models if one assumes that our world has dimension larger then four. The Standard Model is confined to the four-dimensional manifold, while gravity propagates in all dimensions. One further assumes that there is no specific gravity scale, but rather a unique fundamental scale for the Standard Model and gravity around the electroweak scale $`M^{}1`$ TeV . Thus the observed Plank scale is small since the size of extra dimensions is large (up to mm). This does not seem to contradict experimental data since the gravitational potential has been tested at distances down to a cm. Existence of strong gravitational interactions at TeV scale has a number of remarkable experimental signatures . On the contrary, one may assume that the scale of gravity is different from the Standard Model one, but the gravitational lagrangian involves higher derivatives of the metric tensor . In this approach the quantum gravity scale can be as low as $`10^3`$ eV. However, we explicitly neglected all higher order derivatives of metric in (5). Even if the quantum gravity indeed becomes strong at $`M^{}M_{\mathrm{Pl}}`$ the probability that a dilaton escapes into higher dimensions is suppressed by factor of order of $`𝒪(s/M^2)`$, which means that our four-dimensional effective lagrangian (5) is a good approximation at energies $`\sqrt{s}M^{}`$. So, at these energies gravity can be considered as a background classical field, since quantum corrections are small. It would be interesting to address the question of what is the mechanism by which the full theory at the fundamental scale yields the conformally flat GeV-scale gravitational effects. However, we leave this problem aside.

Whatever the underlying mechanism is, the broken scale invariance can be encoded in the effective low–energy theory of gluodynamics in the form of the lagrangian (5). We hope that it provides a useful tool for a systematic approach to non-perturbative processes at high energies.

Acknowledgements We are grateful to J. Ellis, E. Gotsman, A.B. Kaidalov, D.B. Kaplan, Yu. Kovchegov, U. Maor, L. McLerran, M. Pospelov, E. Shuryak, C-I Tan and I. Zahed for illuminating discussions on the subject and helpful comments. The work of D.K. was supported by the U.S. Department of Energy under Contract No. DE-AC02-98CH10886. This research was supported in part by the BSF grant #98000276, by the GIF grant # I-620-22.14/1999 and by Israeli Science Foundation, founded by the Israeli Academy of Science and Humanity. The work of K. T. was sponsored in part by the U.S. Department of Energy under Grant No. DE-FG03-00ER41132. The research of D.K. and K.T. was supported in part by the National Science Foundation under Grant No. PHY94-07194.


text:
# Creation of a molecular condensate by dynamically melting a Mott-insulator

## Abstract

We propose creation of a molecular Bose-Einstein condensate (BEC) by loading an atomic BEC into an optical lattice and driving it into a Mott insulator (MI) with exactly two atoms per site. Molecules in a MI state are then created under well defined conditions by photoassociation with essentially unit efficiency. Finally, the MI is melted and a superfluid state of the molecules is created. We study the dynamics of this process and photoassociation of tightly trapped atoms.

The generation of Bose Einstein condensates (BEC) of dilute atomic gases has resulted in a remarkable series of experiments demonstrating various properties of quantum degenerate gases BEC . One of the next major goals in this effort is the realization of a molecular BEC. A promising route towards a molecular condensate is the conversion of an atomic BEC to molecules via photoassociation, a process discussed so far for conditions of quasihomogeneous trapping of atomic gases Heinzen ; Javanainen ; Holland ; TwoColor . In this Letter we describe a novel path to create condensates of composite atomic objects, in particular a molecular BEC, based on photoassociation via a Mott insulator state of bosonic atoms trapped in an optical lattice Mott ; MottExp . This provides an efficient way of generating a molecular BEC, avoiding some of the problems encountered in the quasihomogeneous case Heinzen . It also touches upon fundamental questions related to the formation of a BEC by “melting” of a Mott-insulator (MI) state in a quantum phase transition, as opposed to the familiar growth from a thermal cloud of atoms growth .

Experimental advances in manipulating BECs BEC , and in particular the loading of a BEC into an optical lattice generated by interfering laser beams have recently led to a seminal experiment by I. Bloch and collaborators MottExp . They demonstrated a quantum phase transition from a BEC or superfluid (SF) state into a MI by varying the lattice laser intensity, as proposed theoretically in Mott . While a SF phase has long range order, the MI phase corresponds to the loading of a precise number of atoms into each lattice site, i.e. Fock state occupation of lattice sites. Among the proposed applications of this new atomic quantum phase are the study of ultracold controlled collisions and quantum computing with neutral atoms CCC . In the present context, the MI phase opens the possibility to efficiently create a molecular BEC in the following four steps: (i) an atomic BEC is loaded into an optical lattice, (ii) the depth $`V_0`$ of the optical lattice is increased to create a MI with exactly two particles per lattice site, (iii) a molecular MI state is produced by two-color photoassociation of the atoms under tight trapping conditions, and (iv) by decreasing the depth of the optical lattice the MI state is “melted”, and thus a molecular BEC is created in a quantum phase transition.

At the end of step (i) above, we have an ensemble of bosonic atoms illuminated by orthogonal, standing wave laser fields tuned far from atomic resonance. These laser fields generate a potential for atomic motion of the form $`V(\stackrel{}{x})=_{i=1}^3V_{0i}\mathrm{sin}^2(kx_i)`$ with $`k=2\pi /\lambda `$ the wave-vector of the light and lattice period $`a=\lambda /2`$. The dynamics of bosonic atoms occupying the lowest Bloch band of an optical lattice is well described by the Bose-Hubbard model (BHM) Mott  which includes the interaction $`U_a`$ between particles occupying the same lattice site and the tunneling $`J_a`$ of particles from one site to the next. The BHM Hamiltonian is given by

$$H_a=J_a\underset{i,j}{}a_i^{}a_j+\underset{i}{}ϵ_i\widehat{n}_i+\frac{1}{2}U_a\underset{i}{}\widehat{n}_i(\widehat{n}_i1)$$
(1)
where $`a_i`$ is a bosonic destruction operator of a particle at site $`i`$. The number of particles in site $`i`$ is given by the operator $`\widehat{n}_i=a_i^{}a_i`$, and $`ϵ_i`$ is an energy offset due to an external trapping potential. The first sum in Eq. (1) runs over all nearest neighbors denoted by $`i,j`$. Increasing the laser intensity of the trapping laser tends to compress atoms near the nodes of the lattice field, and thus leads to an increased on-site interaction $`U_a`$, while the atomic tunneling rate $`J_a`$ decreases Mott . The BHM predicts a quantum phase transition from the SF phase to the MI state: according to mean field theory this occurs at the critical value $`U_a^{(c)}5.8zJ_a`$ MottOld  with $`z`$ the number of nearest neighbors of each site. This corresponds to a (moderate) potential depth of $`V_0=10E_R`$ where $`E_R=\mathrm{}^2k^2/2m`$ is the recoil energy for atoms with mass $`m`$.

We will first illustrate the dynamics of the BHM Hamiltonian with a time dependent depth $`V_0(t)`$ of the optical lattice controlled by the laser intensity, leading to a variation $`U_a(t)`$ and $`J_a(t)`$ in Eq. 1. We assume the system initially to be in the SF ground state and calculate its time evolution for the time dependence shown in Fig. 1a Mott . For a model problem of $`N`$ particles, where $`N`$ is small ($`10`$), in a few lattice sites, the time dependent Schrödinger equation for the wave function $`\mathrm{\Psi }(t)`$ can be solved exactly. Fig. 1b plots the eigenvalues $`e_l`$ of the one particle density matrix $`\rho _{i,j}=a_i^{}a_j`$ ($`\mathrm{\Psi }(t)|a_i^{}a_j|\mathrm{\Psi }(t)`$). As expected, there is one large eigenvalue for $`U_a<U_a^{(c)}`$ of the order of the number of particles. The corresponding wave function is approximately given by $`|\psi _{\mathrm{SF}}(_ia_i^{})^N|\mathrm{vac}`$ for $`N`$ particles in $`M`$ sites with $`|\mathrm{vac}`$ the vacuum state. All the other eigenvalues are small and are associated with the quantum depletion of the SF state. As $`U_a`$ increases and crosses the critical point $`U_a^{(c)}`$ all eigenvalues tend towards one, corresponding to a diagonal single particle density operator (MI). Upon ramping $`U_a`$ down again the SF is restored.

An extension of the BHM (1) in an optical lattice describes the situation with atoms and molecules present. Denoting the annihilation operator of a molecule by $`b_i`$, we add to the Hamiltonian (1) a tunneling $`J_b`$ and on-site interaction term $`U_b`$ for molecules, and an atom-molecule interaction term of the form $`H_{ab}=U_{ab}_ib_i^{}b_ia_i^{}a_i`$. The underlying assumption is that the laser beams generate an optical lattice for atoms and molecules with the same structure of nodes and antinodes, although both lattices can have different depth. The process of atom-molecule conversion by photoassociation is described by the Hamiltonian $`H_{\mathrm{conv}}=\mathrm{\Omega }(t)_i(b_i^{}a_ia_i+\mathrm{h}.\mathrm{c}.)/\sqrt{2}`$, where $`\mathrm{\Omega }(t)`$ is an effective Rabi frequency which is turned on at time $`t=0`$ for a short time interval to convert two atoms at one lattice site into one molecule. In practice, this process will consist of several Raman steps to go to the molecular ground state (see Fig. (2)).

We emphasize several distinguishing key features of the atom-molecule conversion process in an optical lattice. First conversion is most efficient under tight trapping conditions (in a regime where tunneling between lattice sites is negligible). The high atomic densities associated with the strong compression of the atoms opens inelastic collision channels that typically quench all lattice sites with three or more atoms. Thus we assume that only lattice sites with occupation of two (or one) atom survive the lattice compression. These loss processes are added to our Hubbard dynamics by writing down a master equation for a density matrix that includes these decay channels. However, since in the MI phase atomic number fluctuations are small, there are only a very few lattice sites that are actually depleted by this loss process. Second, for (exactly) two atoms trapped at one lattice site we have a complete microscopic understanding of the two-atom dynamics, and the conversion to molecules by photoassociation Mies  beyond the effective description contained in the Hubbard model.

The two-atom Schrödinger equation Mies  at a given lattice site separates for harmonic confinement into center-of-mass and relative coordinates. The potential curves for the relative motion of the two atoms in the trap are schematically shown in Fig. 2: for small distances $`R`$ we have the familiar molecular Born-Oppenheimer potentials while the large $`R`$-behavior is dominated by the trap confinement. The trapping potential discretizes the molecular continuum (scattering) states to form a series of harmonic oscillator trap states with frequency $`\nu =(4V_0E_R)^{1/2}/\mathrm{}`$ with $`V_0`$ the potential depth. The goal is to perform a Rabi flop $`\mathrm{\Omega }T=\pi `$ from the lowest trap state of two atoms (i.e. the state associated with the lowest atomic Bloch band) to a bound molecular state, with $`\mathrm{\Omega }`$ the two-photon Rabi frequency. The discreteness of the trap states makes this a bound-bound Raman transition. The condition for not exciting any other trap states, i.e. to avoid heating, is $`\mathrm{\Omega }\nu `$. On the other hand, for incoherent processes (as spontaneous decay from the intermediate state) with effective decay rate $`\gamma `$ to be small, we must have $`\gamma T1`$. Thus we require a large two-photon Rabi frequency and tight trapping, $`\gamma \mathrm{\Omega }\nu `$. Note that $`\mathrm{\Omega }`$ involves a matrix element from the bound trap state to a molecular state, which results in a scaling $`\mathrm{\Omega }\nu ^{3/4}`$ Mies , i.e. under conditions of tight trapping the two-photon Rabi frequency will be significantly enhanced. We will give specific numbers for these parameters for the case of 87Rb below. Thus it is the preparation of the two-atom MI phase together with strong confinement which guarantees the coherent conversion of atoms to molecules with essentially unit efficiency.

Fig. 1c shows results from the exact integration of the Schrödinger equation for three sites with six atoms. Starting from the SF phase, the atoms are driven to the MI phase and converted into molecules at $`t=0`$. Melting of the molecular MI phase then produces a molecular SF.

To describe the dynamics with a large number of particles in 2D and 3D we employ a time dependent mean field approximation based on a Gutzwiller ansatz MottOld . For simplicity of writing we consider for the moment the case of atoms alone, where we write the wave function as the product of superposition states at the various lattice sites, $`|G(t)=_{i=1}^M\left(_{n=0}^{\mathrm{}}f_n^{(i)}(t)|n_i\right)`$. This ansatz is motivated by the success and simplicity of time-independent Gutzwiller mean field theory to model the ground state and phase diagram of the BHM MottOld . The ground state is obtained from the variational principle $`G|H|G\mu G|\widehat{N}|G\mathrm{min}`$, where $`\mu `$ is a chemical potential introduced to enforce a given mean particle number Mott . From the time-dependent variational principle, $`G(t)|i\mathrm{}\frac{}{t}H(t)|G(t)\mathrm{min}`$, the following time dependent equation is readily derived:

$`i\dot{f}_n^{(i)}`$ $`=`$ $`{\displaystyle \frac{U_a}{2}}n(n1)f_n^{(i)}`$
$`J_a{\displaystyle \underset{i,j}{}}\left(\mathrm{\Phi }_j^{}f_{n+1}^{(i)}\sqrt{n+1}+\mathrm{\Phi }_jf_{n1}^{(i)}\sqrt{n}\right),`$
where $`\mathrm{\Phi }_iG|a_i|G=_nf_{n1}^{(i)}\sqrt{n}f_n^{(i)}`$ is the atomic SF density. Eq. (Creation of a molecular condensate by dynamically melting a Mott-insulator) is a nonlinear equation for the amplitudes $`f_l^{(i)}`$, which preserves both normalization of the wave function and the mean particle number. In the SF limit and for a coherent state distribution of $`f_n^{(i)}=\psi _i^n\mathrm{exp}(\left|\psi _i\right|^2/2)/\sqrt{n!}`$, Eq. (Creation of a molecular condensate by dynamically melting a Mott-insulator) reduces to a time-dependent Gross-Pitaevskii equation for $`\psi _i`$ on a lattice.

By projection to a state with definite particle number $`N`$, $`|G_N=𝒫_N|G/𝒫_N|G_0^{2\pi }𝑑\phi \mathrm{exp}(iN\phi )_{i=1}^M(_{n=0}^{\mathrm{}}\mathrm{exp}(in\phi )f_n^{(i)}|n_i)`$ a more consistent description is obtained. The resulting time dependent Schrödinger equation for the amplitudes $`f_n^{(i)}`$ is significantly more complex. It can be shown, however, that $`f_n^{(i)}`$ of the number projected Gutzwiller wave function again obey Eq. (Creation of a molecular condensate by dynamically melting a Mott-insulator), provided the variance of the particle numbers at each lattice site satisfies $`\mathrm{\Delta }n_i1/\sqrt{N}`$ (where $`\mathrm{\Delta }n_i^2n_i^2n_i^2`$ with $`n_i^2=_nn^2|f_n^{(i)}|^2`$ and $`n_i=_nn|f_n^{(i)}|^2`$). Note that this excludes the regime where we have a *precise* locking of the particle number, i.e., $`f_n^{(i)}=\delta _{n,n_0}`$, as in the MI obtained from the non-number conserving Gutzwiller for a homogeneous situation, when the number of particles $`N`$ is commensurate with the lattice sites $`M`$ and $`n_0=N/M`$. However, a precise Fock state is never realized in the time evolution we consider since the initial superfluid density is not completely destroyed while ramping the optical lattice up. Also, an additional trap potential confining the system to a certain region in space ensures the existence of a remnant superfluid component. Below we model the evolution of an initial SF to an (approximate) MI while changing $`U_a`$ and $`J_a`$ by integrating mean field equations of the type (Creation of a molecular condensate by dynamically melting a Mott-insulator). Fig. 1d gives the results obtained from Gutzwiller theory with initial state given by the time-independent Gutzwiller wave function, to be compared with the exact integration of the time dependent Schrödinger equation in $`1`$D for a few particles in Fig. 1b. As expected, mean field theory shows a more pronounced phase transition than the few atom $`1`$D calculation.

Using a generalization of the Gutzwiller ansatz to include superposition states of atoms and molecules, $`|G(t)=_{i=1}^M\left(_{n_a,n_b=0}^{\mathrm{}}f_{n_a,n_b}^{(i)}(t)|n_a,n_b_i\right)`$, where $`n_a`$ and $`n_b`$ refer to the atomic and molecular occupation, respectively, we have numerically investigated the creation of a molecular BEC in a 2D lattice with a superimposed harmonic trapping potential. The results are shown in Fig. 3. As expected we find that molecules are only created in sites with an atomic occupation of two before the Raman process (cf. Fig. 3b). They are surrounded by a ring of atoms which originates from those site with an atomic occupation of one before the Raman transition (see Fig. 3c). Finally, Fig. 3d shows the superfluid molecular density after ramping the optical lattice down.

We now turn to the description of the coherent Raman transitions involved in creating molecules in step (iii) of our scheme. There are several constraints on the choice of detunings and laser intensities (see Fig. 2) to maximize $`\mathrm{\Omega }`$. First, the detuning $`\mathrm{\Delta }`$ from an allowed excited photoassociation resonance $`|\mathrm{\Psi }_{\mathrm{PR}}`$ (cf. Fig. 2) must be large compared to the natural linewidth $`\gamma _{PR}`$ of the intermediate state to suppress spontaneous Raman scattering. Also, one should not detune half-way between two resonances since interference between these two intermediate states will lead to a minimum in the effective two-photon Rabi frequency $`\mathrm{\Omega }`$. For appropriate choices of $`\mathrm{\Delta }`$ and the intermediate state $`|\mathrm{\Psi }_{\mathrm{PR}}`$, the coupling to all intermediate states other than $`|\mathrm{\Psi }_{\mathrm{PR}}`$ can be neglected. Second, as outlined above, we have to ensure that the process where a $`\omega _2`$ photon is absorbed before a $`\omega _1`$ photon (schematically shown in the inset of Fig. 2) has negligible probability. This process causes trap excitations of single atoms and thus leads to heating. If both of these conditions are fulfilled the effective Rabi frequency $`\mathrm{\Omega }`$ on Raman resonance for the first Raman step is given by $`\mathrm{\Omega }=\mathrm{\Omega }_1\mathrm{\Omega }_2/2\mathrm{\Delta }`$, where $`\mathrm{\Omega }_{1,2}`$ are the Rabi frequencies for the first (second) step, and the effective spontaneous emission rate is $`\gamma =\gamma _{\mathrm{PR}}\mathrm{\Omega }_1^2/4\mathrm{\Delta }^2`$.

For the case of 87Rb, the following Raman pathway is viable for producing $`X^1\mathrm{\Sigma }_g^+(v=0,J=0)`$ molecules:

$`X^1\mathrm{\Sigma }_g^+(v_{\mathrm{trap}}=0)A^1\mathrm{\Sigma }_u(v=213)X^1\mathrm{\Sigma }_g^+(v=120)`$
$`X^1\mathrm{\Sigma }_g^+(v=120)A^1\mathrm{\Sigma }_u(v=185)X^1\mathrm{\Sigma }_g^+(v=52)`$
$`X^1\mathrm{\Sigma }_g^1(v=52)A^1\mathrm{\Sigma }_u(v=24)X^1\mathrm{\Sigma }_g^+(v=0).`$ (3)
A two-step Raman pathway also exists for producing $`a^3\mathrm{\Sigma }_u^+(v=0,J=0)`$ molecules. For the first step of the singlet pathway, the vibrational spacing near the $`A^1\mathrm{\Sigma }_u(v=213)`$ level is 110 GHz, $`\gamma _{PR}12`$ MHz, and the binding energy of the $`X^1\mathrm{\Sigma }_g^+(v=120)`$ is 31.9 GHz. Given a trap frequency $`\nu `$ of 1 MHz and intensities $`I_1=1\text{W}/\text{cm}^2`$ and $`I_2=10^3\text{W}/\text{cm}^2`$ for the first Raman step of Eq. (3), one obtains $`\mathrm{\Omega }_1=0.71`$ MHz and $`\mathrm{\Omega }_2=3.7`$ MHz. For a red detuning of 200 linewidths we get $`\mathrm{\Omega }=1.1`$ kHz and and an effective spontaneous Raman scattering rate of $`\gamma =5.5`$ Hz. For the process described in the preceding paragraph one finds that the undesirable $`\omega _2\omega _1`$ process has a blue detuning of $`29.5`$ GHz off the same intermediate level and a Rabi frequency that is more than a factor of 100 below $`\mathrm{\Omega }`$. Although the pathway for producing singlet molecules is partially optimized, no attempt to optimize the intensities and detunings has been made. However we have chosen values to show that the undesirable $`\omega _2\omega _1`$ process can be sufficiently suppressed. The rate limiting step in the overall pathway is determined by the matrix element of the first Raman step. The matrix elements for all subsequent steps are at least 3-orders of magnitude larger than this one.

For experimental purposes, molecules in higher lying vibrational levels of the $`X^1\mathrm{\Sigma }_g^+`$ or $`a^3\mathrm{\Sigma }_u^+`$ state should also allow the formation and detection of a molecular BEC. Vibrational relaxation caused by inelastic molecule-molecule collisions, that ultimately limits the lifetime of the resulting molecular condensate, will in general be strongly suppressed as the kinetic energy in the exit channels increases. This occurs as the vibrational spacing increases and should allow for observations of molecular BEC for moderately bound vibrational levels.

In conclusion, generation of a MI phase of atoms allows efficient conversion of atoms to molecules, and to obtain a molecular condensate via “melting” in a quantum phase transition. This idea can be immediately generalized to e.g. heteronuclear molecules, or one could use laser chemistry to build more complex composite objects (trimers etc.) and, possibly, corresponding condensates by quantum melting. The key to designing these processes is the fact that the MI phase provides us with a given small number of particles per lattice site (reaction partners) whose few body dynamics can be understood in all detail and controlled via laser interactions.

Discussions with E. Tiesinga, P.S. Julienne and R. Grimm are gratefully acknowledged. Work supported in part by the Austrian Science Foundation, EU Networks, and by the U.S. Office of Naval Research.