text:
# 1 A derivation of the filtration-absorption equation

## 1 A derivation of the filtration-absorption equation

A permeable rock layer on a horizontal impermeable bed is considered. It is well known that the gently sloping groundwater flows in a purely porous medium with infiltration are described by the Boussinesq equation (see refs. (1–3)):

$$m_th=\kappa \left(h_{xx}^2h+(_xh)^2\right)+q.$$
$`[1]`$
Here $`h`$ is the groundwater level, $`t`$ is the time, $`x`$ is the horizontal space coordinate along the impermeable bed, $`q`$ is the intensity of the groundwater inflow or outflow. Furthermore, $`\kappa =\rho gk/\mu `$ is a joint property of the pair-rock/fluid, assumed here to be constant, $`k`$ is the rock permeability, $`m`$ the rock porosity, $`g`$ the gravity acceleration, $`\rho `$ and $`\mu `$ are fluid density and dynamic viscosity.

Assume now that when the fluid level is decreasing, some part of the fluid remains being absorbed by rock, e.g. due to capillary imbibition to the micropores. Then for a fixed fluid particle the rate of fluid absorption $`q`$ will be proportional to the individual time derivative of the fluid level $`dh/dt`$, so that

$$q=\alpha m\frac{dh}{dt}=\alpha m(_th+u_p_xh).$$
$`[2]`$
Here $`\alpha `$ is also a joint rock/fluid property which we also assume here to be constant, and $`u_p`$ is the actual fluid velocity.

According to the Darcy law the filtration velocity (fluid flux per unit area) is equal to

$$u_f=\frac{\rho gk}{\mu }_xh.$$
$`[3]`$
The actual fluid velocity is different (see (1)) from the filtration velocity $`u_f`$. For a purely porous medium $`u_p=u_f/m`$. It is very important that if the rock is fissurized, i.e. contains a connected network of cracks, then

$$u_p=\frac{u_f}{m_1}$$
$`[4]`$
where $`m_1`$ (“fissure porosity”) is much less than block porosity $`m`$: fluid is contained in pores but moves through cracks which are much wider than the pores but occupy much less space. Therefore

$$q=\alpha m\left(_th\frac{\rho gk}{\mu m_1}(_xh)^2\right).$$
$`[5]`$
Substitution of  to the water balance equation  gives the equation

$$m(1\alpha )_th=\frac{\rho gk}{\mu }(h_{xx}^2h+\left(1\alpha \frac{m}{m_1}\right)(_xh)^2),$$

which can be reduced to the form

$$_th=\frac{\kappa }{m(1\alpha )}(h_{xx}^2h(c1)(_xh)^2),$$
$`[6]`$
where

$$c=\alpha \frac{m}{m_1}.$$
$`[7]`$
Here there is a specially important point to be mentioned: If $`m_1=m`$ (non-fissurized purely porous medium), then $`c=\alpha `$ and obviously $`c`$ is always less than one because the absorption cannot exceed the available amount of fluid. However, if the fluid is fissurized, $`m_1`$ can be substantially less than $`m`$, and $`c=\alpha m/m_1`$ can be substantially larger than one. Replacing $`x`$ by $`x/\sqrt{\kappa /m(1\alpha )}`$ and leaving previous notation $`x`$ for transformed space coordinate, we reduce the basic equation to the canonic form

$$_th=h_{xx}^2h(c1)(_xh)^2,$$
$`[8]`$
which will be investigated further.

Note that the proposed model has some common features with the Mirzadzhan-Zadeh model (4) of filtration of the gas-condensate mixture, but is not identical to this model which leads to a different basic equation.

## 2 Self-similar solutions

We look for self-similar solutions with shrinking support and with finite time to collapse (total annihilation): For sufficiently large $`c`$ (in fact, see later, for $`c>3/2`$)

$$h=A(t_0t)^\lambda F\left(\frac{xx_0}{B(t_0t)^\mu }\right).$$
$`[9]`$
Here $`x_0`$ is the point where the solution is collapsing at $`t=t_0`$. We assume the structure of the groundwater “dome” is symmetric, so that $`F(\xi )`$, $`\xi =(xx_0)/B(t_0t)^\mu `$ is an even function. Here $`A,B`$ and $`\lambda `$ are constants, but as we will see later, of a different nature. We will determine the function $`F(\xi )`$ in the interval $`0\xi 1`$, so that $`F(\xi )0`$ at $`\xi 1`$. The quantity

$$x_f=B(t_0t)^\mu $$
$`[10]`$
is the contracting half-width of the groundwater dome. As is well known (5) for the degenerate parabolic differential equations of the type under consideration, the support of the solution remains compact if it is compact initially. Furthermore, the function $`F`$ can be normalized arbitrarily, so we can assume $`A=B^2\mu `$. Substitute  into the basic equation . Bearing in mind that the coefficients of the resulting equation for $`F(\xi )`$ cannot contain the time $`t`$ explicitly, we obtain $`\lambda =2\mu 1`$, and the equation for $`F(\xi )`$ assumes the form:

$$F\frac{d^2F}{d\xi ^2}(c1)\left(\frac{dF}{d\xi }\right)^2\xi \frac{dF}{d\xi }+\frac{2\mu 1}{\mu }F=0.$$
$`[11]`$
We turn now to the boundary conditions. The first one is the condition of symmetry

$$F^{}(0)=0.$$
$`[12]`$
The second boundary condition follows from the continuity of the groundwater level at the free boundary $`x=x_f`$:

$$F(1)=0.$$
$`[13]`$
The last boundary condition follows from the continuity of the groundwater flux at the free boundary $`x=x_f`$. The solution close to free boundary can be considered as quasi-steady: $`h=h(\zeta )`$, $`\zeta =xx_f`$. We obtain from :

$$V\frac{dh}{d\zeta }=h\frac{d^2h}{d\zeta ^2}(c1)\left(\frac{dh}{d\zeta }\right)^2$$
$`[14]`$
where $`V=dx_f/dt`$.

For the solutions of equation  having the quantity $`d(F^2)/d\xi `$ equal to zero at $`\xi =10`$, which is needed to have continuous flux at the free boundary, the quantity $`hd^2h/d\zeta ^2`$ tends to zero at $`x=x_f0`$, so that at $`x=x_f`$

$$V=(c1)\frac{dh}{d\zeta }.$$
$`[15]`$
Bearing in mind that for the self-similar solution , $`V=dx_f/dt=B\mu (t_0t)^{\mu 1}`$, and $`(dh/d\zeta )_{\zeta =0}=(_xh)_{x=x_f}`$, we obtain the relation

$$B\mu (t_0t)^{\mu 1}=B\mu (t_0t)^{\mu 1}(c1)F^{}(1)$$

from which the third boundary condition follows, to be satisfied by the function $`F`$

$$F^{}(1)=\frac{1}{c1}.$$
$`[16]`$
So, a nonlinear eigenvalue problem is obtained: We have to find for the second order equation  in the interval $`[0,1]`$ the solution satisfying three boundary conditions ,,, and the eigenvalue $`\mu `$.

We consider the special case when the function $`F(\xi )`$ has a maximum at $`\xi =0`$. This allows one to search the solution in the form of an expansion

$$F(\xi )=a(1\xi ^2)+\underset{n=2}{\overset{\mathrm{}}{}}a_n(1\xi ^2)^n,$$

so that the terms of the sum (except the first one) do not contribute to all three boundary conditions of the eigenvalue problem. The result is unexpectedly simple:

$$a=\frac{1}{2(c1)},\mu =\frac{c1}{2c3},a_n=0(n2)$$
$`[17]`$
so that, if $`c>3/2`$, the solution to the nonlinear eigenvalue problem is obtained in the form

$$F=\frac{1}{2(c1)}(1\xi ^2),\mu =\frac{c1}{2c3}$$
$`[18]`$
and the self-similar solution under consideration is represented by the relation

$$h=\frac{1}{2(2c3)}B^2(t_0t)^{\frac{1}{2c3}}\left[1\frac{(xx_0)^2}{B^2(t_0t)^{\frac{2(c1)}{2c3}}}\right].$$
$`[19]`$
## 3 Investigation of the self-similar solutions

In spite of its very simple form, the solution  is a typical self-similar solution of the second kind (see (6)). The exponent $`\mu =(c1)/(2c3)`$ cannot be obtained using some conservation laws, but only by solving a nonlinear eigenvalue problem, and the constants $`B`$ and $`t_0`$, as well as $`x_0`$, are obtained by matching the self-similar solution with the solution to the Cauchy problem at the non-self-similar stage. To a certain extent this problem is similar to the problem of the evolution of a turbulent burst (see (6), section 10.2.4). The solution  has essentially different behavior in various intervals of the values of the absorption coefficient $`c`$:

$$0<c<1;1<c<\frac{3}{2};\frac{3}{2}<c.$$
$`[20]`$
The form  is appropriate for the last interval $`3/2<c`$ where the collapse time $`t_0`$ is finite. It is instructive to investigate the limiting behavior of the solution  at $`c3/2`$ from above. Putting $`c=3/2+\epsilon `$ ($`\epsilon >0`$ is a small parameter), we obtain $`\mu =1/4\epsilon +\frac{1}{2}`$, so that

$$(t_0t)^\mu =t_0^{\frac{1}{4\epsilon }+\frac{1}{2}}\left(1\frac{t}{t_0}\right)^{\frac{1}{4\epsilon }+\frac{1}{2}},(t_0t)^{2\mu 1}=t_0^{\frac{1}{2\epsilon }}\left(1\frac{t}{t_0}\right)^{\frac{1}{2\epsilon }},$$

and

$$x_f(t)=Bt_0^{\frac{1}{4\epsilon }+\frac{1}{2}}\left(1\frac{t}{t_0}\right)^{\frac{1}{4\epsilon }+\frac{1}{2}},h(x_0,t)=\frac{B^2}{4\epsilon }t_0^{\frac{1}{2\epsilon }+1}\frac{1}{t_0}.$$
$`[21]`$
Therefore if at $`\epsilon 0`$ the quantity $`4\epsilon t_0`$ tends to a certain constant $`\mathrm{\Theta }`$, and the quantity $`B^2t_0^{\frac{1}{2\epsilon }+1}`$ to another constant which we denote by $`C^2\mathrm{\Theta }`$, the solution  tends to a finite limit:

$$h=C^2e^{2t/\mathrm{\Theta }}\left[1\frac{(xx_0)^2}{C^2\mathrm{\Theta }e^{2t/\mathrm{\Theta }}}\right],x_f=C\sqrt{\mathrm{\Theta }}e^{t/\mathrm{\Theta }}.$$
$`[22]`$
In the interval $`1<c<3/2`$ the exponent $`\mu `$ becomes negative, and it is convenient to replace $`\mu `$ by $`\mu `$, and $`t_0`$ by $`t_0`$. Solution  may be represented in a different form

$$h=\frac{1}{2(32c)}B^2(t_0+t)^{\frac{1}{32c}}\left[1\frac{(xx_0)^2}{B^2(t_0+t)^{\frac{2(c1)}{32c}}}\right],$$
$`[23]`$
so that $`h(x_0,t)=h_{\mathrm{max}}(t)`$ and $`x_f`$ decay with time according to the power laws

$$h_{\mathrm{max}}(t)=\frac{1}{2(32c)}B^2(t_0+t)^{\frac{1}{32c}},x_f(t)=B(t_0+t)^{\frac{c1}{32c}}.$$
$`[24]`$
The time of collapse is infinite and $`t_0`$ becomes simply an additive constant. In the limit $`c3/2`$ from below, $`c=3/2\epsilon `$, $`\epsilon >0`$, we obtain $`\mu =1/4\epsilon 1/2`$, and

$$h_{\mathrm{max}}(t)=\frac{1}{4\epsilon }B^2t_0^{\frac{1}{2\epsilon }+1}\left(1+\frac{t}{t_0}\right)^{\frac{1}{2\epsilon }+1}\frac{1}{t_0},x_f(t)=Bt_0^{\frac{1}{4\epsilon }+\frac{1}{2}}\left(1+\frac{t}{t_0}\right)^{\frac{1}{4\epsilon }+\frac{1}{2}}.$$

Assuming again that at $`\epsilon 0`$ the quantity $`4\epsilon t_0`$ tends to a constant $`\mathrm{\Theta }`$, and $`Bt_0^{\frac{1}{4\epsilon }+\frac{1}{2}}`$ tends to another constant $`C\sqrt{\mathrm{\Theta }}`$, we obtain the same limiting formula .

In the interval $`0<c<1`$ (weak absorption) the compact support extends, not contracts, although slower than in the case of “porous medium equation” $`c=0`$. In this special case $`c=0`$ solution  is reduced to a known self-similar solution fo the first kind ((7),(8); see also (9),(6)). The degenerate special case $`c=1`$ was considered previously; the papers by J. R. King (10) and P. Rosenau (11) should be mentioned specially. It is instructive to compare the results obtained above with those obtained in the paper by B. Meerson et al (12).

## 4 Nonuniqueness of solutions of the Cauchy problem

We consider solutions of equation  with initial condition

$$h(x,0)=h_0(x)\text{for}x.$$
$`[25]`$
where $`h_0(x)`$ is a continuous function which is positive in an interval $`(x_L(0),x_R(0))`$ and which vanishes elsewhere. Let $`x_L(t)`$ and $`x_R(t)`$ be two continuous functions for $`t0`$ such that $`x_L(t)`$ is nondecreasing, $`x_R(t)`$ is nonincreasing, and $`x_L(t)x_R(t)`$ for $`t0`$. It is known that if $`c1`$ for any such pair $`x_L(t)`$ and $`x_R(t)`$ there exists a solution $`h(x,t)`$ of the Cauchy problem , such that $`h(x,t)`$ is positive if $`x_L(t)<x<x_R(t)`$, $`t0`$ and $`h(x,t)`$ vanishes elsewhere. For the proof, the definition of solution, and further references we refer to (13).

Of special interest is the choice of steady interfaces: $`x_L(t)=x_L(0)`$ and $`x_R(t)=x_R(0)`$ for all $`t0`$. The corresponding solution is larger than any other solution, and in (14) a numerical scheme has been introduced which leads to this unique solution. From the modeling point of view it is interesting to observe that this solution can be obtained from the following limiting procedure: replace $`h_0(x)`$ by $`h_0(x)+\epsilon `$ $`(\epsilon >0)`$, solve problem ,, and let $`\epsilon 0`$.

Angenent (15) has constructed a solution of , if $`h_0(x)`$ has nonzero slope at $`x_L(0)`$ and $`x_R(0)`$ (for technical reasons the construction is local in time), which is unique in the class of solutions which can be expanded in a Taylor series of suffcieintly high degree near the interfaces:

$$h(x,t)=\underset{k=0}{\overset{N}{}}c_k(t)(xx_L(t))^k+o((xx_L(t))^N)\text{as}xx_L(t);$$

a similar expression holds at the right interface $`x=x_R(t)`$. Here the uniqueness not only refers to $`h`$, but also to the interfaces. We observe that the self-similar solution  belongs to this class of solutions.

In section 5 we shall construct a numerical scheme which yields solutions converging to the self-similar solution for large times, and it is natural to ask if these solutions belong to the class introduced by Angenent. In particular this scheme yields solutions which are different from the ones obtained by the scheme in (14). We conjecture that the solutions which we construct in the following section are physically relevant, but undoubtedly future research is needed to provide definite answers to the nonuniquness question.

## 5 Numerical experiment

The goal of the numerical experiment was to indicate that the self-similar solution obtained above attracts the solutions to non-self-similar Cauchy problems having the initial condition of compact support, generally speaking a non-symmetric one:

$$h(x,0)=h_0(x),x_L(0)xx_R(0),$$

and $`h(x,0)0`$ outside the interval $`x_L(0)xx_R(0)`$. The basic equation  can be transformed to a form convenient for numerical calculations

$$_th=\frac{\dot{x}_0(t)\dot{x}_f(t)\xi }{x_f(t)}_\xi h+\frac{1}{[x_f(t)]^2}[h_{\xi \xi }^2h(c1)(_\xi h)^2],$$

where $`\dot{x}_0=dx_0/dt`$, $`\dot{x}_f=dx_f/dt`$, and

$$\xi =\frac{xx_0(t)}{x_f(t)},x_f(t)=\frac{x_R(t)x_L(t)}{2},x_0(t)=\frac{x_R(t)+x_L(t)}{2},$$

so that the interval of new space variable $`\xi `$ becomes fixed: $`1\xi 1`$, whereas the solution $`h(x,t)`$ is different from zero in the time dependent interval $`x_L(t)xx_R(t)`$. Again, assuming naturally the quasi-steadiness of the level distribution in the vicinities of free boundaries, we obtain the conditions

$$\dot{x}_R(t)=(c1)\frac{_\xi h(1,t)}{x_f(t)},\dot{x}_L(t)=(c1)\frac{_\xi h(1,t)}{x_f(t)}$$
$`[26]`$
at $`t>0`$, and the basic equation takes the form:

$$_th=\frac{1}{x_f^2(t)}\left[(c1)_\xi h\frac{(\xi +1)_\xi h(1,t)(\xi 1)_\xi h(1,t)}{2}+h_{\xi \xi }h(c1)(_\xi h)^2\right]$$
$`[27]`$
with the initial condition

$$h(\xi ,0)=h_0(\xi ),|\xi |1;h(\xi ,0)0,|\xi |1.$$
$`[28]`$
Two numerical schemes were used in our computations performed by finite-difference approximations: (i) a forward-in-time, centered-in-space explicit approximation, and (ii) a forward-in-time, centered-in-space implicit approximation. For the most part numerical calculations have been run with the time step $`\mathrm{\Delta }t=10^5`$ for the explicit scheme and $`\mathrm{\Delta }t=10^4`$ for the implicit one. The number $`N`$ of subintervals of length $`\mathrm{\Delta }\xi `$, $`N=2/\mathrm{\Delta }\xi `$ was equal to 202 for both schemes. The results obtained by using these numerical approximations coincided with good accuracy. The absorption coefficient $`c`$ was always equal to 1.75.

The first initial condition was taken as a “smoothed block”: a homogeneous water level distribution smoothly going to zero at the edges. The results of the computation are presented in Figure 1 in the form of the distribution of the scaled level: level $`h(x,t)`$ divided by maximum level at each time $`h_{\mathrm{max}}(t)`$. It is seen that the curves corresponding to different times collapse to the parabola, corresponding to the self-similar solution . The time of collapse $`t_0`$ and the constant $`B`$ were determined in the following way (Figure 2): According to the intermediate asymptotics  at small $`t_0t`$,

$$\frac{x_f^2(t)}{h_{\mathrm{max}}(t)}=\frac{t_0t}{\mu F(0)}.$$
$`[29]`$
so that the quantity $`x_f^2(t)/h_{\mathrm{max}}(t)`$ should be a linear function of time, and the intersection of its graph with the time axis (Figure 2,a) gives the value of $`t_0`$. We obtain from  a linear relation between $`\mathrm{ln}x_f`$ and $`\mathrm{ln}(t_0t)`$, i.e. in the coordinates $`\mathrm{ln}x_f,\mathrm{ln}(t_0t)`$ a straight line with the slope $`\mu `$. It gives us the value $`B`$ (Figure 2,b) and an additional possibility of checking the asymptotics. Naturally $`t_0`$ and $`B`$ depend on the initial condition. For our case we found $`t_0=0.345`$, $`B=3.73`$ and $`\mu =1.499`$, which agrees well with the analytic value $`\mu =1.5`$.

The next computation was performed for a nonsymmetric initial condition:

$$h(x,0)=\{\begin{array}{cc}4x^2+4x,\hfill & 0<x<\frac{1}{2}\hfill \\ \frac{4}{9}x^2+\frac{4}{9}x+\frac{8}{9},\hfill & \frac{1}{2}<x<2.\hfill \end{array}$$

Figures 3,a and 3,b demonstrate the behavior of the numerical solution for different times. It is clearly seen that the solution becomes symmetric and tends to the self-similar asymptotics . The values of $`t_0=1.138`$, $`B=0.63`$ and $`\mu =1.5027`$ have been calculated as before; the calculated and analytic values of $`\mu `$ agree with high precision.

For comparison we have taken the solution  for a certain $`t`$ as an initial condition, and computed the solution to the partial differential equation further using the same algorithm. The results are presented in Figure 4,a for different times. Being plotted in scaled coordinates (Figure 4,b) they collapse to a single curve, giving us an additional check of the numerical procedure.

## 6 Conclusions

We presented a new derivation of the filtration-absorption equation based on a model of groundwater flow with partial absorption. It is shown that for a sufficiently large absorption constant the time of collapse is finite. A family of self-similar solutions to this equation is obtained. Numerical experiments indicate that these self-similar solutions obtained are self-similar intermediate asymptotics for the solutions to the Cauchy problems having the initial conditions data with compact support, but due to the nonuniqueness of the solution of the Cauchy problem future research is needed to provide more definite conclusions.

We express our gratitude to Professor A. J. Chorin and Professor R. Dal Passo for their valuable comments. We thank Professor S. Abarbanel for his interest in our work. This work was supported in part by the National Science Foundation under Grant DMS 97–32710, and in part by the Applied Mathematics subprogram of the U.S. Department of Energy under contract DE–AC03–76–SF00098.

1. Polubarinova-Kochina, P. Ya. (1962). Theory of Groundwater Movement, Princeton University Press, Princeton.
2. Bear, J. (1972). Dynamics of Fluids in Porous Media, Dover, New York.
3. Barenblatt, G. I., Entov, V. M., and Ryzhik, V. M. (1990). Flow of Fluids Through Natural Rocks, Kluwer Academic, Dordrecht, The Netherlands.
4. Magerramov, N. Kh., and Mirzadzhan-Zadeh (1960). Journal of Applied Mathematics & Mechanics (Prikl. Mat. Mekh.), vol. XXIV, no. 6.
5. Kalashnikov, A. S. (1987). Russian Math. Surveys 42, 169–222.
6. Barenblatt, G. I. (1996). Scaling, Self-similarity, and Intermediate Asymptotics, Cambridge University Press, Cambridge.
7. Zeldovich, Ya. B., and Kompaneets, A. S. (1950). On the theory of propagation of heat with thermal conductivity depending on temperature. In Collection of Papers Dedicated to the 70th Birthday of A. F. Ioffe, pp.61–71 Izd. Akad. Nauk USSR, Moscow.
8. Barenblatt, G. I. (1952). Prikl. Mat. Mekh. 16, no. 1, 67–78.
9. Zeldovich, Ya. B., and Raizer, Yu. P. (1967). Physics of Shock Waves and High Temperature Hydrodynamic Phenomena, vol. 1, Academic Press, New York.
10. King, J. R. (1993). Journ. Engng. Math. 27, no. 2, 31–72.
11. Rosenau, P. (1995). Phys. Rev. Letters 74, no. 7, 1056–1059.
12. Meerson, B., Sasorov, P.V., and Sekimoto, K. (2000) Phys. Rev. Letters E 61, no. 2, 1403–1406.
13. Bertsch, M., Dal Passo, R., and Ughi, M. (1992). Annali di Matematica pura ed applicata (IV), vol. CLXI, 57–81.
14. Bertsch, M. and Dal Passo, R. (1990). Quarterly Appl. Math. 48, 133–152.
15. Angenent, S. (1980). Math. Annalen 280, 465–482.

Figure Captions

Figure 1. The numerical solution to the Cauchy problem for $`c=1.75`$ with the initial condition of a “smoothed block” type for different times in the scaled coordinates. The solution is collapsing to the self-similar asymptotics.

Figure 2.
(a) The determination of the parameter $`t_0`$.

(b) The determination of the parameter $`B`$.

Figure 3.
(a) Evolution of the nonsymmetric initial distribution to a symmetric self-similar
asymptotics.

(b) The same evolution presented in scaled coordinates.

Figure 4. The numerical solution preserves the self-similarity.


format: text

text:
# On finitization of the Gordon identities

## Abstract

In this paper we prove the identity that generalizes the Andrews–Gordon identity. Also we discuss the relation of our formula to the geometry of affine flag varieties and to the geometry of polyhedra.

This work is partially supported by grants RFBR 99-01-01169, INTAS-OPEN-97-1312

1. Introduction.

Let $`\widehat{SL_2}`$ be the group corresponding to the affine Lie algebra $`\widehat{sl_2}`$. This group is a central extension of the group of matrices $`2\times 2`$ with determinant one over the ring of functions on a circle. By $`\widehat{𝒩}`$ denote the subgroup of $`\widehat{SL_2}`$ which consists of unipotent upper–triangular matrices, by $`\widehat{n}`$ denote the Lie algebra of $`\widehat{𝒩}`$. In other words, the group $`\widehat{𝒩}`$ consists of currents on the circle with values in the maximal nilpotent subgroup of $`SL_2`$. The algebra $`\widehat{n}`$ is abelian and it can be identified with the space of functions on the circle. Let us choose a basis $`\{e_i\}`$, $`i`$ in $`\widehat{n}`$, such that the element $`e_m`$ corresponds to the function $`\mathrm{e}^{2\pi im\varphi }`$, where the coordinate $`\varphi `$ is the angle parameter on the circle. Let $`Gr`$ be the affine grassmanian for the group $`\widehat{SL_2}`$, i. e. the quotient of $`\widehat{SL_2}`$ by a maximal parabolic subgroup $`P`$. Suppose that $`e_i`$ belongs to the Lie algebra of $`P`$ if $`i0`$. By $`pGr`$ denote the image of $`P`$ in the variety $`Gr`$.

In the paper  the closure of the $`\widehat{𝒩}`$–orbit of the point $`p`$ on $`Gr`$ is studied. Denote this variety by $``$. The geometric definition of the variety $``$ is rather complicated, but the description of $``$ in terms of its coordinate ring seems more simple and clear to us.

Recall that there exists a line bundle $`\xi `$ on the variety $`Gr`$ such that the space of sections of $`\xi ^k`$ is a realization of the irreducible vacuum representation $`L_{0,k}`$ at level $`k`$ of the algebra $`\widehat{sl_2}`$. More precisely, $`H^0(Gr,\xi ^k)(L_{0,k})^{}`$. So the coordinate ring of the variety $`Gr`$ is the direct sum of the spaces $`(L_{0,k})^{}`$ with the natural product $`(L_{0,k_1})^{}(L_{0,k_2})^{}(L_{0,k_1+k_2})^{}`$.

The variety $``$ is a subvariety of $`Gr`$, so its coordinate ring is a quotient of the coordinate ring of $`Gr`$. Namely, let $`v_k`$ be the highest vector of $`L_{0,k}`$ with respect to the subgroup $`P`$. The space $`L_{0,k}`$ contains the subspace $`D(k)=U(\widehat{n})v_k`$. Then the coordinate ring of $``$ is the direct sum of the spaces $`D(k)^{}`$ with the induced product.

In integrable representations of $`\widehat{sl_2}`$ at level $`k`$ we have the relation $`e(t)^{k+1}=0`$, when $`e(t)=e_it^i`$. In  it is shown that this relation leads to the following functional realization of the spaces $`D(k)^{}`$.

$$D(k)^{}=\mathrm{\Sigma }_0(k)\mathrm{\Sigma }_1(k)\mathrm{\Sigma }_2(k)\mathrm{},$$

where $`\mathrm{\Sigma }_m(k)`$ is the space of symmetric polynomials in $`m`$ variables $`t_1,\mathrm{}t_m`$ which are divisible by the product $`t_1\mathrm{}t_m`$ and have a zero if $`t_1=\mathrm{}=t_{k+1}`$.

In this notation the product on the direct sum of $`D(k)^{}`$ is defined as the set of maps $`\mathrm{\Sigma }_{m_1}(k_1)\mathrm{\Sigma }_{m_2}(k_2)\mathrm{\Sigma }_{m_1+m_2}(k_1+k_2)`$. The image of a pair of polynomials $`P_1(t_1,\mathrm{}t_{m_1})`$ and $`P_2(t_1,\mathrm{},t_{m_2})`$ under this map is the symmetrization with respect to the variables $`t_1,\mathrm{},t_{m_1+m_2}`$ of their product:

$$(P_1P_2)(t_1,\mathrm{},t_{m_1+m_2})=\mathrm{Symm}(P_1(t_1,\mathrm{},t_{m_1})P_2(t_{m_1+1},\mathrm{},t_{m_1+m_2})).$$
(1)
Let $`N`$ be an integer. Let $`\widehat{n}_N`$ be the subalgebra of $`\widehat{n}`$ spanned by $`e_i`$ for $`i>N`$, let $`\widehat{𝒩}_N`$ be the corresponding subgroup. We want to define the variety $`_N`$ as the quotient of the variety $``$ by the group $`\widehat{𝒩}_N`$. From the geometric point of view, to do it one should remove all the non–stable orbits from the variety $``$ and factorize the resulting space by the action of the group. From the algebraic point of view, the coordinate ring of $`_N`$ is the subring of $`\widehat{𝒩}_N`$–invariants (or, just in other words, $`\widehat{n}_N`$–invariants) in the coordinate ring of $``$.

Thus the coordinate ring of $`_N`$ is the direct sum of spaces $`D_N(k)^{}`$ dual to the spaces $`D_N(k)=\widehat{n}_ND(k)`$ of coinvariants of $`D(k)`$ with respect to the algebra $`\widehat{n}_N`$. The product on these spaces is defined in the similar way. Now define the variety $`_N`$ as the projective spectrum of the graded algebra $`D_N(k)^{}`$. In this way we obtain a line bundle $`\xi _N`$ on $`_N`$, such that $`H_0(_N,\xi _N^k)D_N(k)^{}`$.

On the variety $`Gr`$, there is an action of the two–dimensional torus $`T`$, whose tangent space in the unit is spanned by the Cartan element $`h_0`$ and the energy element $`L_0`$. It is easy to show that the action of $`T`$ on the varieties $``$ and $`_N`$ is also defined and the action of $`T`$ on $`D(k)`$ and $`D_N(k)`$ is defined by $`[h_0/2,e_i]=e_i`$, $`[L_0,e_i]=ie_i`$.

Let

$$d_N(k;q,z)=\mathrm{Tr}_{D_N(k)}q^{L_0}z^{h_0/2}$$

be the Hilbert polynomial of the graded space $`D_N(k)`$. There are two ways to compute it, and it is shown in  that they both succeed for the spaces $`D(k)`$. The first way is purely algebraic. It is based on the functional realization of $`D_N(k)^{}`$ and the Gordon filtration on the space of symmetric polynomials. This computation was done in , and the result is written in the left hand side of (2). The second way is an application of the holomorphic Lefschetz formula to the bundle $`\xi _N^k`$. Let us discuss it in detail.

It can be shown that for even $`N`$ the variety $`_N`$ is non–singular (but for odd $`N>3`$ it contains exactly one singular point) and $`H^i(_N,\xi _N^k)=0`$ for $`i>0`$. In this case the fixed points of the action of the torus $`T`$ on $`_N`$ can be enumerated by two non–negative integers $`m`$ and $`n`$, such that $`m+nN/2`$. The fixed point $`x_{0,0}`$ is the image of the point $`pGr`$ in $`_N`$. For $`t`$ consider the point

$$x_{m,n}(t)=\underset{i=0}{\overset{m1}{}}\mathrm{exp}(te_{2i+1})\underset{i=0}{\overset{n1}{}}\mathrm{exp}(te_{N2i})x_{0,0}.$$

Then the fixed point $`x_{m,n}`$ is the limit $`\underset{t\mathrm{}}{lim}x_{m,n}(t)`$. Adding the contributions of these points up we obtain the right hand side of (2). This expression is also written in .

So for even $`N`$ we have the following formula.

$$\underset{v_1,\mathrm{},v_k_+}{}q^{\underset{i,j}{}v_iv_j\mathrm{min}(i,j)}z^{\underset{i}{}iv_i}\underset{i=1}{\overset{k}{}}\left[\genfrac{}{}{0pt}{}{\left(N+1\right)i2\underset{j}{}v_j\mathrm{min}(i,j)+v_i}{v_i}\right]_q=$$
(2)
$$=\underset{\genfrac{}{}{0pt}{}{m,n_+}{m+nN/2}}{}\frac{z^{k\left(m+n\right)}q^{k\left(m^2+\left(N+1\right)nn^2\right)}}{\underset{i=1}{\overset{m}{}}\left(1q^i\right)\left(1q^{im+1}z^1\right)\underset{i=2m+1}{\overset{N2n}{}}\left(1q^iz\right)\underset{i=1}{\overset{n}{}}\left(1q^i\right)\left(1q^{N+i+n2}z^1\right)},$$

where $`\left[\genfrac{}{}{0pt}{}{n}{m}\right]_q`$ is the $`q`$–binomial coefficient and the notation $`\frac{1}{1q^\alpha z^\beta }`$ should be considered as $`\underset{i0}{}q^{i\alpha }z^{i\beta }`$ when $`\alpha ,\beta 0`$ and as $`\left(\underset{i<0}{}q^{i\alpha }z^{i\beta }\right)`$ when $`\alpha ,\beta 0`$.

The right hand side of this formula resembles the alternating sum in the Weyl formula. Recall that the Weyl formula for characters of integrable representations of a Kac–Moody algebra can be interpreted as the Lefschetz formula applied to bundles on the corresponding flag variety. The fixed points of the action of the torus on the flag variety can be enumerated by the elements of the Weyl group, so they correspond to the extremal vectors of an integrable representation. And the contribution of a fixed point to the Lefschetz formula can be interpreted as the character of the representation ”in a small vicinity” around the corresponding extremal vector. Let us describe it for representations of $`\widehat{SL_2}`$.

Consider the integrable representation $`L_{l,kl}`$ of $`\widehat{SL_2}`$, let $`v_{k,l}`$ be its highest vector. Consider the extremal vector $`w(v_{k,l})`$, when $`w`$ is an element of the Weyl group of $`\widehat{SL_2}`$. Denote the degrees of this vector by $`c_q(w;k,l)`$ and $`c_z(w;k,l)`$. We call the limit

$$q^{c_q(w;k,l)}z^{c_z(w;k,l)}\underset{r\mathrm{}}{lim}\left(q^{c_q(w;rk,rl)}z^{c_z(w;rk,rl)}\mathrm{ch}L_{rl,r(kl)}\right)$$
(3)
the character of the representation $`L_{l,kl}`$ ”in a small vicinity” around the vector $`w(v_{k,l})`$. It is easy to show that it coincides with the image of the character of the Verma module $`M_{l,kl}`$ under the action of $`w`$. To obtain the contribution of the vector $`w(v_{k,l})`$ we should expand the formula (3) in positive powers of $`q`$. In this way we obtain the character of a certain Verma module with the sign equal to $`(1)^{l(w)}`$. The Weyl formula states that the character of the representation $`L_{l,kl}`$ is the sum of such contributions.

In this paper we use this idea to calculate the characters of the spaces $`D_N(k)`$. Namely, we decompose the character of $`D_N(k)`$ into the sum of characters ”in small vicinities” of some vectors. In this way we generalize the formula (2) and prove the generalization by combinatorial methods. Indeed, the bosonic part of our formula is also the Lefschetz formula applied to bundles on a certain variety.

2.1 Consider the space

$$D_N(k,l,r)=[e_1,\mathrm{},e_N]/I_N(k,l,r),N2,$$

where the ideal $`I_N(k,l,r)`$ is generated by the elements

$$e_1^{l+1},e_N^{r+1},\underset{\genfrac{}{}{0pt}{}{1\alpha _1,\mathrm{},\alpha _{k+1}N}{\alpha _1+\mathrm{}+\alpha _{k+1}=i}}{}e_{\alpha _1}\mathrm{}e_{\alpha _{k+1}},\text{when}k+1iN(k+1).$$

In other words, $`I_N(k,l,r)`$ is generated by the elements

$$e_1^{l+1},e_N^{r+1},e(t)^{k+1},t,\text{when}e(t)=\underset{i=1}{\overset{N}{}}e_it^i.$$

In particular, the space $`D_N(k)`$ discussed in the introduction is isomorphic to $`D_N(k,k,k)`$.

On the space $`D_N(k,l,r)`$ we introduce a bi–grading $`(\mathrm{deg}_q,\mathrm{deg}_z)`$. To do it we are bi–grading the algebra $`[e_1,\mathrm{},e_N]`$ by setting $`\mathrm{deg}_qe_i=i`$, $`\mathrm{deg}_ze_i=1`$. As the generators of $`I_N(k,l,r)`$ are homogeneous, this bi–grading on $`[e_1,\mathrm{},e_N]`$ defines a bi–grading on $`D_N(k,l,r)`$. By $`d_N(k,l,r;q,z)`$ denote the Hilbert polynomial of this space:

$$d_N(k,l,r;q,z)=\underset{i,j0}{}q^iz^jdim(D_N(k,l,r))^{i,j}.$$

2.2 The space dual to $`D_N(k,l,r)`$ has the following functional realization.

Let $`\mathrm{\Sigma }_m^N`$ be the space of symmetric polynomials of the form

$$t_1t_2\mathrm{}t_mf(t_1,t_2,\mathrm{},t_m)$$
(4)
with degree not exceeding $`N`$ in each of $`m`$ variables. By $`\mathrm{\Sigma }_m^N(k,l,r)`$ denote the subspace of $`\mathrm{\Sigma }_m^N`$ which consists of polynomials of the form (4), such that

1. $`f(t_1,\mathrm{},t_m)=0`$ if $`t_1=\mathrm{}=t_{k+1}`$.
2. $`f(t_1,\mathrm{},t_m)=0`$ if $`t_1=\mathrm{}=t_{l+1}=0`$.
3. A polynomial $`f(t,\mathrm{},t,t_{r+2},t_{r+3},\mathrm{},t_m)`$ has the degree less then $`(r+1)(N1)`$ with respect to the variable $`t`$.

###### Proposition 1

The space $`D_N(k,l,r)`$ is dual to the direct sum

$$\mathrm{\Sigma }_0^N(k,l,r)\mathrm{\Sigma }_1^N(k,l,r)\mathrm{\Sigma }_2^N(k,l,r)\mathrm{}.$$

In these settings the subspace $`(D_N(k,l,r))^{i,j}`$ is dual to the subspace of $`\mathrm{\Sigma }_j^N(k,l,r)`$ which consists of homogeneous polynomials of degree $`i`$.

The construction of the duality is the following. Let $`S_N^m`$ be the space of homogeneous polynomial of homogeneous degree $`m`$ in $`N`$ variables $`e_1,\mathrm{},e_N`$. Define the map $`𝒫:(S_N^m)^{}\mathrm{\Sigma }_m^N`$ by setting for $`\varphi (S_N^m)^{}`$

$$𝒫(\varphi )(t_1,\mathrm{},t_m)=\varphi (e(t_1)\mathrm{}e(t_m)).$$

To complete the prove of the proposition it is enought to observe that this map is an isomorphism and that the restriction 1) on the subspace $`\mathrm{\Sigma }_m^N(k,l,r)`$ corresponds to factorizing the dual space by $`e(t)^{k+1}`$, the restriction 2) corresponds to factorizing by $`e_1^{l+1}`$ and the restriction 3) corresponds to factorizing by $`e_N^{r+1}`$.

2.3 Note that we can define an associative commutative product on the space $`D_N^{}=\mathrm{\Sigma }_m^N(k,l,r)`$. Namely, the product of polynomials $`P_1\mathrm{\Sigma }_{m_1}^N(k_1,l_1,r_1)`$ and $`P_2\mathrm{\Sigma }_{m_2}^N(k_2,l_2,r_2)`$ can be defined by the formula (1). It is clear that $`P_1P_2\mathrm{\Sigma }_{m_1+m_2}^N(k_1+k_2,l_1+l_2,r_1+r_2)`$.

The structure of commutative algebra on the space $`D_N^{}`$ implies the following statement.

###### Proposition 2

There exists a variety $`_N^{}`$ together with line bundles $`\xi _K`$, $`\xi _L`$, $`\xi _R`$ and vector fields $`L_q`$, $`L_z`$ on $`_N^{}`$, such that the space $`D_N(k,l,r)`$ is dual to the space $`H^0(_N^{},\xi _K^k\xi _L^l\xi _R^r)`$ and $`d_N(k,l,r;q,z)`$ is the trace of the action of $`q^{L_q}z^{L_z}`$ on this space.

The construction of the variety $`_N^{}`$ is the following. Let $`S`$ be the spectrum of the algebra $`D_N^{}`$. The algebra $`D_N^{}=D_N(k,l,r)^{}`$ is polygraded by integers $`k`$, $`l`$ and $`r`$, so there is an action of three–dimensional complex torus $`T_3`$ on $`S`$. Let $`\stackrel{~}{}_N^{}`$ be the maximal subvariety of $`S`$ where the torus $`T_3`$ acts freely. Set $`_N^{}=\stackrel{~}{}_N^{}/T_3`$. Then we obtain the bundles $`\xi _K`$, $`\xi _L`$, $`\xi _R`$ as the bundles associated to the corresponding one–dimensional representations of the torus $`T_3`$.

2.4 Using the Gordon filtration on the spaces $`\mathrm{\Sigma }_m^N`$ (see , ) we obtain the following formula that generalize the left hand side of (2).

###### Theorem 1

Let $`Q`$ be the $`k\times k`$ matrix with $`Q_{i,j}=\mathrm{min}(i,j)`$; let $``$ and $``$ be the vectors in $`^k`$ such that

$$_i=\{\begin{array}{cc}il\hfill & il\hfill \\ 0\hfill & il\hfill \end{array}_i=\{\begin{array}{cc}ir\hfill & ir\hfill \\ 0\hfill & ir\hfill \end{array}$$

Then $`d_N(k,l,r;q,z)=`$

$$=\underset{\genfrac{}{}{0pt}{}{v=(v_1,\mathrm{},v_k)}{v_i_+}}{}q^{v^tQv+^tv}z^{{\scriptscriptstyle iv_i}}\underset{i=1}{\overset{k}{}}\left[\genfrac{}{}{0pt}{}{(N+1)i(2Qv++v)_i}{v_i}\right]_q,$$

where the symbol t is just the transposition and

$$\left[\genfrac{}{}{0pt}{}{n}{m}\right]_q=\frac{\underset{i=1}{\overset{n}{}}(1q^i)}{\underset{i=1}{\overset{m}{}}(1q^i)\underset{i=1}{\overset{nm}{}}(1q^i)}.$$

The proof of this formula repeats the proof of similar statements in  and .

2.5 Next we construct a basis in the space $`D_N(k,l,r)`$.

###### Theorem 2

The monomials of the form $`e_1^{a_1}\mathrm{}e_N^{a_N}`$ such that

$$a_1l,a_Nr,a_i+a_{i+1}k\text{for any}1i<N,$$
(5)
form a basis in $`D_N(k,l,r)`$.

Proof.

In  this statement is proved for $`r=k`$ and, as $`D_N(k,l,0)D_{N1}(k,l,k)`$, for $`r=0`$.

First of all, from the fact that the monomials of the form $`e_1^{a_1}\mathrm{}e_N^{a_N}`$ with the restrictions $`a_i+a_{i+1}k`$ span the space $`D_N(k,k,k)`$ it follows that the monomials with the restrictions (5) span $`D_N(k,l,r)`$. So to complete the proof it is sufficient to compare the dimension of the space $`D_N(k,l,r)`$ and the number of elements of the presumed monomial basis.

By $`\stackrel{~}{D}_N(k,l,r)`$ denote the set of monomials with the restrictions (5). Also let $`d_N(k,l,r)=dimD_N(k,l,r)`$ and let $`\stackrel{~}{d}_N(k,l,r)`$ be the number of elements in $`\stackrel{~}{D}_N(k,l,r)`$. We have

###### Lemma 1

1. $`\stackrel{~}{d}_N(k,l,r)=\stackrel{~}{d}_N(k,l,r1)+\stackrel{~}{d}_{N1}(k,l,kr)`$.
2. $`d_N(k,l,r)d_N(k,l,r1)+d_{N1}(k,l,kr)`$.

Proof of the lemma.

(i) Let $`m=e_1^{a_1}\mathrm{}e_N^{a_N}\stackrel{~}{D}_N(k,l,r)`$. Then either $`a_N<r`$ and $`m\stackrel{~}{D}_N(k,l,r1)`$, or $`a_N=r`$ and $`m^{}=e_1^{a_1}\mathrm{}e_{N1}^{a_{N1}}\stackrel{~}{D}_{N1}(k,l,kr)`$. And conversely, for any element of $`\stackrel{~}{D}_N(k,l,r1)`$ or $`\stackrel{~}{D}_{N1}(k,l,kr)`$ we can find the corresponding element in $`\stackrel{~}{D}_N(k,l,r)`$.

(ii) Introduce the map $`\stackrel{~}{i}_N(r):[e_1,\mathrm{},e_{N1}]D_N(k,l,r)`$ defined by $`\stackrel{~}{i}_N(r)(P)=Pe_N^r`$. It is easy to see that the image of the ideal $`I_{N1}(k,l,kr)`$ under this map is zero, so $`\stackrel{~}{i}_N(r)`$ induces the map $`i_N(r):D_{N1}(k,l,kr)D_N(k,l,r)`$. In particular we obtain $`dim\mathrm{Im}(i_N(r))d_{N1}(k,l,kr)`$.

It is also clear that the image of the ideal $`I_N(k,l,r1)`$ under the natural projection $`D_N(k,l,r)\mathrm{Coker}(i_N(r))`$ is zero, so $`dim\mathrm{Coker}(i_N(r))d_N(k,l,r1)`$. As $`dimD_N(k,l,r)=dim\mathrm{Im}(i_N(r))+dim\mathrm{Coker}(i_N(r))`$, we have the statement of the lemma.

From the equality $`d_N(k,l,r)=\stackrel{~}{d}_N(k,l,r)`$ for $`r=0,k`$ and also for $`N=2`$ it follows that the inequality in statement (ii) of lemma 1 is in fact an equality. So we have the equality $`d_N(k,l,r)=\stackrel{~}{d}_N(k,l,r)`$ for all $`r`$ and the statement of the theorem.

###### Corollary 1

$`d_N(k,l,r;q,z)=d_N(k,l,r1;q,z)+(q^Nz)^rd_{N1}(k,l,kr;q,z)`$.

2.6 Let us present another corollary from the theorem. Let $`P_N(k,l,r)`$ be the polyhedron in $`^N=\{(x_1,\mathrm{},x_N)\}`$ defined by the inequalities

$$x_1l;x_Nr;x_i+x_{i+1}k\text{for any}1i<N.$$

Consider the functionals $`\phi _z=x_i`$, $`\phi _q=ix_i`$ on $`^N`$.

###### Corollary 2

The dimension of the space $`D_N(k,l,r)`$ coincides with the number of integer points of the polyhedron $`P_N(k,l,r)`$, and

$$d_N(k,l,r;q,z)=\underset{𝐱^NP_N(k,l,r)}{}q^{\phi _q(𝐱)}z^{\phi _z(𝐱)}.$$

Following the ideas described in  we decompose the polynomial $`d_N(k,l,r;q,z)`$ into the sum of contributions of the vertices of the polyhedron $`P_N(k,l,r)`$. Namely, we propose (see theorem 3) an identity of the form

$$d_N(k,l,r;q,z)=\underset{M\text{ is a vertex of }P_N(k,l,r)}{}d_M.$$

Let $`M`$ be a vertex. In a small vicinity of $`M`$ our polyhedron looks like a cone over a certain $`N1`$–dimensional polyhedron. Denote this cone by $`C(M)`$. In our case for the vertex $`M`$ we set

$$d_M=\underset{𝐱^NC(M)}{}q^{\phi _q(𝐱)}z^{\phi _z(𝐱)}.$$

According to  we say that a vertex $`M`$ is simple if $`C(M)`$ is a cone over a simplex. In this case there exists a simple formula for the contribution. It is easy to show that for generic $`k,l,r`$ (namely, for $`k>l>0`$, $`k>r>0`$, $`lr`$, $`l+rk`$) the simple vertices of $`P_N(k,l,r)`$ are

$$M_{m,n}=(\underset{m}{\underset{}{l,kl,l,kl,\mathrm{}}},\underset{Nmn}{\underset{}{0,\mathrm{},0}},\underset{n}{\underset{}{\mathrm{},r,kr,r}}).$$

In this formula $`n+mN`$, and if $`n+m=N`$ then the sum of $`m`$th and $`(m+1)`$th coordinates should not exceed $`k`$.

###### Proposition 3

Suppose that $`k`$, $`l`$ and $`r`$ are generic. Then any integer point of $`C(M_{m,n})`$ can be obtained from $`M_{m,n}`$ by adding a linear combination with non–negative integer coefficients of the vectors $`v^1,\mathrm{},v^N`$, when

$$v^i=\{\begin{array}{cc}(\underset{mi}{\underset{}{0,\mathrm{},0}},\underset{i}{\underset{}{1,1,\mathrm{},\left(1\right)^{i1},\left(1\right)^i}},\underset{Nm}{\underset{}{0,\mathrm{},0}})\hfill & im,\hfill \\ (\underset{i1}{\underset{}{0,\mathrm{},0}},1,\underset{Ni}{\underset{}{0,\mathrm{},0}})\hfill & m<iNn,\hfill \\ (\underset{Nn}{\underset{}{0,\mathrm{},0}},\underset{Ni+1}{\underset{}{\left(1\right)^{Ni+1},\left(1\right)^{Ni},\mathrm{},1,1}},\underset{n+iN1}{\underset{}{0,\mathrm{},0}})\hfill & Nn<i\hfill \end{array}$$

So

$$d_{M_{m,n}}=q^{\phi _q(M_{m,n})}z^{\phi _z(M_{m,n})}\underset{i=1}{\overset{N}{}}\frac{1}{1q^{\phi _q(v^i)}z^{\phi _z(v^i)}}.$$

The quantities $`\phi _q(M_{m,n})`$ and $`\phi _z(M_{m,n})`$ are linear functions in $`k`$, $`l`$, $`r`$. The quantities $`\phi _q(v^i)`$ and $`\phi _z(v^i)`$ are independent on $`k`$, $`l`$, $`r`$ and equal to the following numbers.

|  | $`mi`$ | | $`Nni`$ | $`i>Nn`$ | |
| --- | --- | --- | --- | --- | --- |
|  | $`i`$ is even | $`i`$ is odd | $`i>m`$ | $`Ni`$ is even | $`Ni`$ is odd |
| $`\phi _z\left(v^i\right)`$ | 0 | -1 | 1 | -1 | 0 |
| $`\phi _q\left(v^i\right)`$ | $`i/2`$ | $`\left(i1\right)/2m`$ | $`i`$ | $`n1\left(3Ni\right)/2`$ | $`\left(Ni+1\right)/2`$ |

Thus $`d_{M_{m,n}}=`$

$$=\frac{q^{\phi _q\left(M_{m,n}\right)}z^{\phi _z\left(M_{m,n}\right)}}{\underset{i=1}{\overset{\left[m/2\right]}{}}\left(1q^i\right)\underset{i=\left[m/2\right]+1}{\overset{m}{}}\left(1q^iz^1\right)\underset{i=m+1}{\overset{Nn}{}}\left(1q^iz\right)\underset{i=1}{\overset{\left[n/2\right]}{}}\left(1q^i\right)\underset{i=\left[n/2\right]+1}{\overset{n}{}}\left(1q^{iN1}z^1\right)}$$
(6)
To consider this formula as a series in positive powers of $`q`$ and $`z`$ let us set the term $`\frac{1}{1q^\alpha z^\beta }`$ equal to $`\underset{i0}{}q^{i\alpha }z^{i\beta }`$ if $`\alpha ,\beta 0`$ and equal to $`\left(\underset{i<0}{}q^{i\alpha }z^{i\beta }\right)`$ if $`\alpha ,\beta 0`$ (these are all the possible cases). Denote the resulting series by $`d_N^{m,n}(k,l,r;q,z)`$. Formula (6) defines it for any integer $`k`$, $`l`$, $`r`$.

2.7 It turns out that in our case it is enought to sum up the contributions of the simple vertices ($`M_{m,n}`$), that is, the sum of the contributions of the other vertices is zero. Let us set

$$d_N^e(k,l,r;q,z)=\underset{\genfrac{}{}{0pt}{}{m+n<N\text{or}}{m+n=N,m\text{is even}}}{}d_N^{m,n}(k,l,r;q,z),$$

$$d_N^o(k,l,r;q,z)=\underset{\genfrac{}{}{0pt}{}{m+n<N\text{or}}{m+n=N,m\text{is odd}}}{}d_N^{m,n}(k,l,r;q,z).$$

###### Theorem 3

1. If $`N`$ is odd and $`lr`$ then $`d_N(k,l,r;q,z)=d_N^o(k,l,r;q,z)`$,
2. If $`N`$ is odd and $`lr`$ then $`d_N(k,l,r;q,z)=d_N^e(k,l,r;q,z)`$,
3. If $`N`$ is even and $`l+rk`$ then $`d_N(k,l,r;q,z)=d_N^o(k,l,r;q,z)`$,
4. If $`N`$ is even and $`l+rk`$ then $`d_N(k,l,r;q,z)=d_N^e(k,l,r;q,z)`$.

Proof. First, prove some lemmas.

###### Lemma 2

1. $`d_N^e(k,l,r;q,z)=d_N^e(k,l,r1;q,z)+(q^Nz)^rd_{N1}^e(k,l,kr;q,z)`$
2. $`d_N^o(k,l,r;q,z)=d_N^o(k,l,r1;q,z)+(q^Nz)^rd_{N1}^o(k,l,kr;q,z)`$

Proof of the lemma. It follows from the identity

$$d_N^{m,n}(k,l,r;q,z)=d_N^{m,n}(k,l,r1;q,z)+(q^Nz)^rd_{N1}^{m,n1}(k,l,kr;q,z)$$

for $`n>0`$ and the identity $`d_N^{m,0}(k,l,r;q,z)=d_N^{m,0}(k,l,r1;q,z)`$.

###### Lemma 3

$`d_N^e(k,l,1;q,z)=0`$ for odd $`N`$, $`d_N^o(k,l,1;q,z)=0`$ for even $`N`$.

Proof of the lemma. It follows from the identity

$$d_N^{m,2n}(k,l,1;q,z)=d_N^{m,2n+1}(k,l,1;q,z).$$

###### Lemma 4

If $`l+r=k`$ then $`d_{2N}^e(k,l,r;q,z)=d_{2N}^o(k,l,r;q,z)`$.

Proof of the lemma.

We want to prove that $`\underset{m=0}{\overset{2N}{}}(1)^md_{2N}^{m,2Nm}(k,a,ka;q,z)=0`$ for any $`k`$ and $`a`$. Let

$$\mathrm{\Phi }(n)=\underset{i=1}{\overset{n}{}}\frac{1}{1q^i},\mathrm{\Psi }(m,n)=\underset{i=m}{\overset{n}{}}\frac{1}{1q^iz^1},\mathrm{\Phi }^{}(n)=\underset{i=1}{\overset{n}{}}\frac{1}{1q^i}.$$

Then

$$d_{2N}^{2s,2N2s}(k,a,ka;q,z)=P_{q,z}\mathrm{\Phi }(s)\mathrm{\Psi }(s+1,N+s)\mathrm{\Phi }^{}(Ns),$$

$$d_{2N}^{2s+1,2N2s1}(k,a,ka;q,z)=P_{q,z}\mathrm{\Phi }(s)\mathrm{\Psi }(s+1,N+s+1)\mathrm{\Phi }^{}(Ns1),$$

where $`P_{q,z}`$ is the common numerator of these fractions that doesn’t depend on $`s`$. It easy to prove by induction that for $`s<N`$ we have

$$\underset{m=0}{\overset{2s}{}}\left(1\right)^md_{2N}^{m,2Nm}(k,a,ka;q,z)=\frac{P_{q,z}}{1q^N}\mathrm{\Phi }\left(s\right)\mathrm{\Psi }(s+1,N+s)\mathrm{\Phi }^{}\left(Ns1\right),$$

$$\underset{m=0}{\overset{2s+1}{}}\left(1\right)^md_{2N}^{m,2Nm}(k,a,ka;q,z)=\frac{P_{q,z}}{1q^N}\mathrm{\Phi }\left(s\right)\mathrm{\Psi }(s+2,N+s+1)\mathrm{\Phi }^{}\left(Ns1\right).$$

Using the second formula for $`s=N1`$ we obtain the statement of the lemma.

Now prove the theorem. For $`N=2`$ the proof is straightforward. Here we already have two cases depending on $`l`$ and $`r`$. Suppose that all the statements of the theorem are proved for $`N1`$ and prove them for $`N`$.

Let $`N`$ be odd. By lemmas 2 and 3 we have $`d_N^e(k,l,0;q,z)=d_{N1}^e(k,l,k;q,z)`$. As $`D_N(k,l,0)D_{N1}(k,l,k)`$, the induction hypothesis implies the equality $`d_N(k,l,0;q,z)=d_N^e(k,l,0;q,z)`$. Induction on $`r`$ based on lemma 2 and corollary 1 shows that $`d_N(k,l,r;q,z)=d_N^e(k,l,r;q,z)`$ if $`rl`$, which is statement (ii).

Statement (i) follows from (ii) by applying the transformation of the space $`D_N(k,l,r)`$ mapping $`e_i`$ to $`e_{Ni+1}`$.

Now let $`N`$ be even. The proof of statement (iii) is similar to the proof of (ii). For $`r=kl`$ statement (iv) follows from statement (iii) and lemma 4. Then statement (iv) for an arbitrary $`rkl`$ can be proved by similar induction on $`r`$.

2.8 To obtain the right hand side of (2) consider the particular case $`l=r=k`$. In this case the vertices $`M_{2m,2n}`$, $`M_{2m1,2n}`$, $`M_{2m,2n1}`$ and $`M_{2m1,2n1}`$ glue into one vertex (let us denote it by $`M_{m,n}^{}`$) and, moreover, for odd $`N`$ all the vertices of the form $`M_{n,Nn}`$ glue into one vertex (let us denote it by $`M^{}`$).

We set

$$d_{M_{m,n}^{}}=\underset{i,j=0,1}{}d_N^{2mi,2nj}(k,k,k;q,z),$$

assuming that $`d_N^{i,j}(k,k,k;q,z)=0`$ for $`i<0`$ and for $`j<0`$. Then the following statement implies the right hand side of (2).

###### Proposition 4

$$d_{M_{m,n}^{}}=\frac{z^{k\left(m+n\right)}q^{k\left(m^2+\left(N+1\right)nn^2\right)}}{\underset{i=1}{\overset{m}{}}\left(1q^i\right)\left(1q^{im+1}z^1\right)\underset{i=2m+1}{\overset{N2n}{}}\left(1q^iz\right)\underset{i=1}{\overset{n}{}}\left(1q^i\right)\left(1q^{N+i+n2}z^1\right)}$$

Also let $`d_M^{}=\underset{i+j=N,i\text{is even}}{}d_N^{i,j}(k,k,k;q,z)`$ for odd $`N`$. We propose the following statement that describes the contribution of the singular point on $`_N`$.

###### Conjecture 1

Let $`N`$ be odd. Then

$$d_M^{}=\frac{P_N(q,z)}{\underset{i=1}{\overset{N}{}}(1q^iz^1)},$$

when $`P_N(q,z)`$ is a polynomial.


text:
# YOUNG CIRCUMSTELLAR DISKS NEAR EVOLVED MASSIVE STARS AND SUPERNOVAE

## 1 INTRODUCTION

There is growing evidence that low mass stars form in association with massive stars. Infrared studies of Orion show evidence for a cluster of about 2,000 low mass stars within 2 pc of the Trapezium; the core radius of the cluster is $`0.2`$ pc (Hillenbrand & Hartmann 1998). Stars in Orion within $`0.3`$ pc of the central star $`\theta ^1`$ C Ori show observable disks because of the effects of photoionization by the massive star (O’Dell, Wen, & Hu 1993). The disk ages may be limited by mass loss driven by a photodissociation front (Störzer & Hollenbach 1999). However, the age of the association of $`10^6`$ years and the fact that disks are apparently present within $`0.01`$ pc of $`\theta ^1`$ C Ori shows that the disks can survive, perhaps because their host stars are on radial orbits (Störzer & Hollenbach 1999). Absorption studies of $`\pi `$ Sco in the upper Sco association give evidence for a population of circumstellar disks that have properties like those in Orion (Bertoldi & Jenkins 1992). The estimated age of the upper Sco association is $`5\times 10^6`$ years for both low and high mass stars, with a small age dispersion (de Geus 1992; Preibisch & Zinnecker 1999). de Geus (1992) presents evidence that a supernova occurred in the association $`(11.5)\times 10^6`$ years ago; it is possible that these circumstellar disks have been exposed to a nearby supernova.

There have been estimates of the lifetime of the solar nebula $`10^6`$ years, but Podosek & Cassen (1994) argue that a lifetime $`10^7`$ years or more is plausible. These longer lifetimes, together with evidence for clusters of low mass pre-main sequence stars near massive stars, suggests that many gaseous disks may be exposed to the late evolutionary phases and supernovae of massive stars. The lifetime of a $`60M_{}`$ (initial mass), or O5 main sequence, star is $`4\times 10^6`$ years and that of a $`20M_{}`$ (O9) star is $`1\times 10^7`$ years. In their late phases, these stars go through a red supergiant and/or a Wolf-Rayet phase, with strong stellar winds. They end their lives as supernovae.

In view of the fact that many circumstellar disks may be exposed to these unusual environments, the aim of this paper is to investigate their effects on the disks. This is of interest both for possible observational consequences of the disks and for possible implications for the early solar system, if it went through such a phase. The effects of radiation on disks and of flows past a disk are discussed in § 2 and § 3, respectively. The possible implications for the early solar system are in § 4.

## 2 EFFECTS OF RADIATION

The effects of a main sequence stellar radiation field on a disk have been treated in detail for the case of Orion. Störzer & Hollenbach (1999) find that the FUV (far ultraviolet) stellar flux determines the mass loss through a photodissociation front. In this regime the mass loss rate does not have a strong dependence on the stellar radiation flux. However, once the radiation flux drops below the flux at which the FUV emission dominates, the EUV (extreme ultraviolet) emission dominates and the mass loss rate drops with the stellar flux. For the case of $`\theta ^1`$ C Ori, Störzer & Hollenbach (1999) estimate that the transition occurs at a distance of $`0.3`$ pc from the O star, which is the outer limit at which disks are observed. The rate of ionizing photon production, $`S_{}`$, from $`\theta ^1`$ C Ori is $`2.6\times 10^{49}`$ s-1, within the range of $`2\times 10^{48}`$ to $`8\times 10^{49}`$ s-1 for O9 to O4 main sequence stars (Panagia 1973). The Orion case may thus be typical of expectations around an O star.

Similar conditions are relevant to the case of Wolf-Rayet stars, if that is the next phase of evolution of the massive star. The case of $`\gamma ^2`$ Velorum is of interest because it is the brightest Wolf-Rayet star and X-ray observations show evidence for a surrounding cluster of pre-main sequence stars (Pozzo et al. 2000). There is no evidence at present for disks around these stars, the presence of which would be important for establishing disk lifetimes. For a distance of 410 pc, these stars extend to a radial distance from $`\gamma ^2`$ Velorum of 2 pc. $`\gamma ^2`$ Velorum is a WC8 star in a binary system with an O9 I star (Conti & Smith 1972). The ultraviolet radiation field is dominated by the O9 I star, which has $`S_{}1.3\times 10^{49}`$ s-1 (Panagia 1973). Although this is comparable to the case of $`\theta ^1`$ C Ori, the stronger wind expected for the Wolf-Rayet star could affect the appearance of nebulae around disks (see § 3).

When the massive star becomes a supernova, the highest luminosity occurs when the shock front breaks out of the stellar surface and there is a burst of hard radiation. In the case of a red supergiant progenitor, $`2\times 10^{48}`$ ergs of hard ultraviolet/soft X-ray radiation is emitted over $`10^3`$ s (Klein & Chevalier 1978; Matzner & McKee 1999). For a Type Ib or Ic supernova from a Wolf-Rayet star, the smaller surface area leads to $`3\times 10^{44}2\times 10^{46}`$ ergs of soft X-ray radiation emitted over $`220`$ s (Matzner & McKee 1999). The case of the unusual supernova SN 1987A is intermediate between these. The initial burst can be followed by prolonged X-ray emission if the shock wave runs into a dense circumstellar wind.

The supernova radiation can heat and ionize gas in the disk, leading to a delayed burst of line radiation. I consider a disk with mass $`M_d=0.01M_{d2}M_{}`$ and radius $`d=10^{15}d_{15}`$ cm, so that the typical surface density is $`\sigma _d=M_d/(\pi d^2)=6.4M_{d2}d_{15}^2`$ g cm-2. The disk scale height at a radius $`0.5d`$ is $`H0.7\times 10^{14}d_{15}^{3/2}T_2^{1/2}(M_{}/M_{})^{1/2}`$ cm, where $`T`$ is the temperature of the disk in units of $`10^2`$ K and $`M_{}`$ is the mass of the host star, so that the typical density in the disk is $`\rho _d1\times 10^{14}d_{15}^{7/2}T_2^{1/2}M_{d2}(M_{}/M_{})^{1/2}`$ g cm-3. For the explosion of a red supergiant with a $`2\times 10^{48}`$ erg initial burst of radiation, the amount of radiative energy intercepted by a disk is $`1\times 10^{42}(r/0.2\mathrm{pc})^2`$ ergs, where $`r`$ is the distance of the disk from the massive star. At $`r=0.2`$ pc, the supernova radiation is absorbed in the outer layers of the disk and reradiated; if $`n_H=5\times 10^9`$ cm-3, an ionized layer $`2\times 10^{12}`$ cm thick can come into equilibrium with the ionizing luminosity. If present, dust grains would compete with gas for the ionizing photons, but they are likely to be evaporated by the high luminosity (see below). Based on calculations by Chevalier & Fransson (1994), I estimate that several % of the reradiated luminosity is in the H$`\alpha `$ line. The timescale for the emission at high density is determined by the light travel time, $`7\times 10^3`$ s, leading to an H$`\alpha `$ luminosity of $`6\times 10^{36}(r/0.2\mathrm{pc})^2`$ ergs s-1 which could be observed up to $`1.3(r/0.2\mathrm{pc})`$ yr after the initial explosion.

Although this emission is faint compared to the supernova luminosity, there is some chance of detection because it is in a narrow line. The case we have been able to observe in most detail is SN 1987A, which apparently had a $`20M_{}`$ progenitor star with a lifetime of $`1\times 10^7`$ years (Arnett et al. 1989). This is close to the age limit at which substantial disks are expected to be present. Efremov (1991) finds that the supernova is at the edge of a star cluster with the expected age, but it is sparse (see also Panagia et al. 2000). A number of short-lived, discrete sources of H$`\alpha `$ emission were observed near the supernova in the first two years after the explosion (Cumming & Meikle 1993 and references therein). Although its parameters ($`d2\times 10^{15}`$ cm and $`r0.5`$ pc) are reasonable for a circumstellar disk, the best observed emission knot (Cumming & Meikle 1993) appears unlikely to be a disk because of the $`11\mathrm{km}\mathrm{s}^1`$ redshift of the narrow H$`\alpha `$ line compared to the centroid of the bright ring emission. Cumming & Meikle (1993) estimate a density of $`(12)\times 10^5`$ cm-3 from the 60 day timescale of the emission; this is low for a circumstellar disk, but it may be that only the surrounding, low density parts of the disk were heated and ionized. SN 1987A was a compact and relatively low luminosity supernova so that dust should survive at this distance from the supernova, but the dust opacity through the ionized region is not expected to be large. Cumming & Meikle (1993) find that the energy in the H$`\alpha `$ line was consistent with the expected ionizing radiation from the supernova.

The supernova radiation can process dust in the disk (by vaporization and melting), and the surviving dust can emit scattered light as well as reradiate infrared light from absorbed supernova radiation. Vaporization is caused by the highest luminosity radiation that is absorbed by dust. From the results of Dwek (1983) for $`n=1`$, where the dust opacity is $`\kappa \lambda ^n`$, and an evaporation temperature of $`T_v=1500`$ K, the dust is evaporated out to $`r_v0.3L_{45}^{1/2}`$ pc from the supernova, where $`L_{45}`$ is the luminosity in units of $`10^{45}\mathrm{erg}\mathrm{s}^1`$. The initial ionizing burst has $`L_{45}2`$ so that a disk at $`r=0.2`$ pc has its dust withing the ionized region evaporated. The luminosity rapidly drops, so that a layer just interior to the evaporated region is expected to have its dust melted. The properties of this layer depend on the disk position and the grain properties.

Dust in circumstellar disks can reprocess the non-ionizing supernova light, which typically has a total energy $`2\times 10^{49}`$ ergs and $`L_{45}<0.01`$ for a Type II supernova, to scattered and reradiated infrared light. The problem with detecting this emission is that it is not transient and is likely to be difficult to distinguish from dust emission from the presupernova wind (Dwek 1983). The disks are optically thick to the supernova light, so the fraction of the supernova radiative energy that is reprocessed depends on the area covering factor of circumstellar disks. Due to the initial dust evaporation, the massive star wind is likely to be optically thin to the supernova light, but it has a 100% covering factor.

## 3 EFFECTS OF FLOWS

In addition to strong radiation fields, the environment of an evolved massive star can have flows which affect a nearby circumstellar disk. In the case of main sequence stars, the properties of the cometary nebulae observed around disks in Orion can be explained by the effects of radiation from the massive star (e.g., McCollough et al. 1995; Störzer & Hollenbach 1999). Although there is evidence for a stand-off bow shock in the wind from $`\theta ^1`$ C Ori in some cases (McCollough et al. 1995), the wind does not play a role in driving mass loss from the disks. The mass loss rate from Wolf-Rayet stars is typically 100 times larger than that from massive main sequence stars. Their wind velocities are comparable so that the wind ram pressure is 100 times larger for the Wolf-Rayet case at a given distance from the stars. If the radiatively driven flows are comparable in the two cases, the wind bow shock is then 10 times closer to the disk in the Wolf-Rayet case, at a radius of $`10^{15}`$ cm, which places it at about the same radius as the photoionization front for nebulae like those in Orion (McCollough et al. 1995). The wind can thus affect the optical appearance of the nebula, by making it somewhat more compact.

For a strong, fast flow, like that of a nearby supernova, matter may be stripped from a disk. If the timescale for the flow interaction is longer than the dynamical timescale for the disk, the ram pressure of the flow can come into equilibrium with the gravitational forces maintaining the disk. If the interaction is rapid, the question is whether the momentum in the flow can cause disk material to reach escape velocity. The escape velocity from the outer parts of the disk is $`v_{esc}=(2GM_{}/d)^{1/2}=5.2(M_{}/M_{})^{1/2}d_{15}^{1/2}\mathrm{km}\mathrm{s}^1`$, where where $`G`$ is the gravitational constant. The disk dynamical timescale is $`t_d=d/v_{esc}=61(M_{}/M_{})^{1/2}d_{15}^{3/2}`$ yr.

In ram pressure stripping, the disk is disrupted if the ram pressure in the flow, $`\rho _fv_f^2`$ where $`\rho _f`$ is the local density in the flow and $`v_f`$ is the velocity, exceeds the gravitational force per unit area, $`p_{grav}`$, that keeps the disk bound to the central pre-main sequence star. An estimate of the gravitational force per unit area is $`p_{grav}GM_{}\sigma _d/d^22\times 10^3\left(M_{}/M_{}\right)M_{d2}d_{15}^4\mathrm{dynes}\mathrm{cm}^2.`$

The supernova case is complicated by the uncertainties in the density distributions of the supernova and the surrounding medium. In order to obtain an estimate of the effects, I assume a constant density supernova with ejecta mass $`M_{ej}`$ and energy $`E`$ in a constant density medium with H density $`n_0`$ (assuming a 10:1 H to He ratio by number). The ram pressure effect is largest if the disk is directly exposed to the freely expanding supernova ejecta, i.e. it becomes placed within the reverse shock wave of the supernova remnant. For a constant density supernova, the maximum radius of the reverse shock is $`r_{rs}=4.4\left(M_{ej}/10M_{}\right)^{1/3}n_0^{1/3}\mathrm{pc}`$ (Truelove & McKee 1999). For the values of $`r`$ considered here, the circumstellar disk is exposed to the expanding ejecta. The maximum ram pressure that can be exerted by the ejecta occurs if the disk is hit by the outer, undecelerated edge of the ejecta, which moves with a velocity $`v_e=(10E/3M_{ej})^{1/2}=4.1\times 10^3E_{51}^{1/2}(M_{ej}/10M_{})^{1/2}\mathrm{km}\mathrm{s}^1`$, where $`E_{51}`$ is in units of $`10^{51}`$ ergs. The supernova interaction time is $`t_{SN}r/v_e=240(r/\mathrm{pc})E_{51}^{1/2}(M_{ej}/10M_{})^{1/2}`$ yr, which is equal to the disk dynamical timescale, $`t_d`$, for $`r0.25`$ pc. For $`r0.25`$ pc, the ram pressure stripping arguments are relevant, but for a close supernova, the momentum in the ejecta is the important factor. The peak ram pressure is $`p_{ram}=5E/(2\pi r^3)=3\times 10^5E_{51}(r/\mathrm{pc})^3\mathrm{dynes}\mathrm{cm}^2,`$ which is to be compared to $`p_{grav}`$ above. A disk at 1 pc from the supernova can survive even this extreme case. At distances smaller than 0.25 pc, momentum transfer causes stripping to occur, beginning at the outer parts of the disk. The criterion for stripping is now that $`M_{ej}v_e/(4\pi r^2)>\sigma _dv_{esc}`$. It can be seen that the ram pressure and momentum stripping criteria are roughly the same when the age equals the disk dynamical time, or $`r/v_ed/v_{esc}`$.

Stripping is the most likely consequence of a strong flow past a disk, but there may be circumstances under which some of the flow is accreted to the disk. The problem with the fast winds from main sequence or Wolf-Rayet stars or with supernova ejecta is that the wind gas is shock heated to a high temperature ($`10^7[v_f/1000\mathrm{km}\mathrm{s}^1]^2`$ K for a H rich gas) and the radiative cooling time is longer than the flow time. Unless the gas can mix with the cool disk gas and share heat with it, it is implausible that fast flow gas can be added to a disk. The slow wind from a red supergiant star can marginally cool in a flow time and gravitational effects can be significant for its accretion. The Bondi-Hoyle accretion rate for a supersonic flow is $`\dot{M}4\pi (GM_{}/v_f^2)^2\rho _wv_f`$, where $`\rho _w`$ is the local density in the wind and $`v_f`$ is the wind velocity. Integrating over the time of accretion, the total accreted mass is $`2\times 10^9(M_{}/M_{\mathrm{}})^2(v_f/10\mathrm{km}\mathrm{s}^1)^4(r/\mathrm{pc})^2M_l`$, where $`M_l`$ is the total mass lost during the red supergiant phase.

## 4 THE EARLY SOLAR SYSTEM

In view of the evidence for the birth of low mass stars close to massive stars, it is interesting to speculate on the consequences for the proto-solar system if it were born is such an environment. Störzer & Hollenbach (1999) note that the disks close to $`\theta ^1`$ C Ori ($`r0.3`$ pc) can have their outer parts photoevaporated by the strong radiation field from the star. They further note that if this applied to the solar system, it could help explain why Uranus and Neptune have considerably less hydrogen than Jupiter and Saturn.

There are further consequences of the solar system being born at $`0.2`$ pc from a massive star. When the massive star became a red supergiant, at an age of $`(34)\times 10^6`$ years, the solar system would have been enveloped by a slow dense wind. Gravitational accretion of the wind would lead to $`1.4\times 10^7M_{}`$ of dusty gas being accreted on the solar system, if the massive star lost about $`3M_{}`$ of material during this phase. The late accretion of grains would contribute matter that would not have to pass through the possibly destructive environment of the early solar nebula formation, although the accretion process could lead to grain destruction. The wind material would contain radioactive 26Al, which is known to have been present at the birth of the solar system (Lee et al. 1976). The abundance of 26Al in the wind from the H envelope is $`2\times 10^6`$ by mass for a $`25M_{}`$ star (Meyer, Woosley, & Weaver 1995), so that $`3\times 10^{13}M_{}`$ of 26Al might be accreted. This is less than the $`3\times 10^9M_{}`$ of 26Al that is required if the 26Al was mixed throughout the proto-solar system including the sun (Cameron et al. 1995), but it is sufficient to contaminate $`1\times 10^4M_{}`$ of disk gas to the observed level. The initial mass of the protoplanetary disk must be $`10^2M_{}`$, but only part of it may contain 26Al. Other extinct radioactivities (e.g., Cameron et al. 1995 and references therein) are produced deeper in the massive star and are not ejected in the red supergiant wind unless there is a mechanism to mix them into the outer layers. The inner layers would be ejected past the solar system during the Wolf-Rayet and/or supernova phases, but I have argued above that the high velocity of the gas in these phases makes accretion difficult. If the gas could stick, the picture would resemble the early “fly-paper model” of T. Gold (see p. 267 in Clayton 1977).

As discussed in § 3, a layer of melted grains may be formed soon after the time of supernova shock wave breakout, depending on the disk position and grain properties. This is of interest for solar system chondrules, which are grains that have been melted (Hewins 1997). In the supernova case, the grains would be rapidly heated by the arrival of radiation at the time of shock breakout, but would cool more slowly because of the continued supernova radiation. There is some evidence that the chondrules were rapidly heated and subsequently cooled on a longer timescale (Hewins 1997). The determination of whether supernova radiation is a suitable heat source will require more detailed calculations.

The supernova gas would have reached the solar nebula tens of years after the explosion. If the disk had already been limited to a size of 10 a.u. by the earlier action of photoevaporation, the passage of the ejecta would not further disrupt the disk, but it would drive a shock front into the disk. The ratio of $`p_{grav}`$ discussed above to the thermal pressure in the disk $`p_{grav}/p_{th}v_{esc}/c_d`$ is typically $`>10`$, where $`c_d`$ is the sound speed in the disk. The pressure due to ejecta approaches $`p_{grav}`$ and can drive a shock wave. The shock wave propagation depends on the detailed density distribution in the disk. Shock waves in the solar nebula have been a leading explanation for the formation of chondrules (Boss 1996) and provide an alternative to the radiative heating. If cooling of the shocked material is slow, some ablation of the disk gas is possible. The supernova shock wave in the red supergiant wind generates a high pressure region that envelops the circumstellar disk for tens of years. For a wind mass loss rate of $`3\times 10^6M_{}\mathrm{yr}^1`$, $`v_f=10\mathrm{km}\mathrm{s}^1`$, and a shock velocity of $`v_{sh}=4,000\mathrm{km}\mathrm{s}^1`$, the pressure is $`\rho _wv_{sh}^2=6\times 10^6`$ dynes cm-2, $`10^7`$ times the present interstellar pressure. The cosmic ray energy density is likely to be similarly enhanced because of shock acceleration of particles; relativistic electrons in such shocked layers are observed in the radio supernova phenomenon. Cosmic rays have been suggested as a source of isotopic anomalies in meteorites (e.g., Clayton & Jin 1995), but the energy available in the current situation does not appear to be sufficient for interesting effects.

This speculative scenario of coeval formation of the solar system and a nearby massive star is an alternative to the hypothesis that the formation of the solar system was triggered by a supernova and that radioactive isotopes were injected at that time (e.g., Cameron et al. 1995; Foster & Boss 1996, 1997). A major difference is that in the present scenario, the early solar system is considerably closer to the supernova and is exposed to a more extreme environment. It is able to survive because of the gravitational binding of a disk to its central star. The unusual environment has a bearing on several perplexing properties of solar system material. The probability that the solar system went through such a phase is small but perhaps non-negligible. If the Orion Nebula Cluster has a limiting radius of 2 pc, $`0.1`$ of the solar mass stars are within 0.2 pc of the center (Hillenbrand & Hartmann 1998). The cluster is likely to become sparser as is ages (Hillenbrand & Hartmann 1998), so the probability of a nearby massive star may decline by the time of a supernova. This probability will be better determined by further observations of circumstellar disks in massive star clusters and by studies aimed at finding disks near massive stars in their final evolutionary phases. The proposal made here of a search for transient narrow H$`\alpha `$ line emission in the spectra of Type II supernovae is one possibility, although it is only near the more massive supernovae ($`20M_{}`$ initial mass) that young disks are likely to be present.

I am grateful to Donald Clayton and Zhi-Yun Li for comments on the manuscript and to the referee for an especially useful report. Support for this work was provided in part by NASA grant NAG5-8088.


text:
# Charmonium + Light Hadron Cross Sections

## 1 Introduction

Many unusual subjects have been studied in the name of QCD. One of the more unusual, which has arisen in the field of heavy ion collisions, is the size of cross sections of charmonia on light hadrons. This has attracted attention because of the proposal by Matsui and Satz  that suppression of $`J/\psi `$ production could be used as a signature for the formation of a quark-gluon plasma.

This suggestion, like many signatures proposed for the quark-gluon plasma, is perhaps excessively intuitive. The idea is that a QGP will screen the linear confining interaction between quarks, so that a $`c\overline{c}`$ pair produced within a QGP will be less likely to form a bound $`c\overline{c}`$ charmonium resonance, as in Fig.1, but instead will more likely separate to form open-charm mesons.

Even if this simple picture of $`c\overline{c}`$ production in a QGP is qualitatively correct, it can only be confirmed easily if the competing direct charm production and scattering by the initial relativistic nucleons is understood  and if there is little subsequent dissociation of the charmonia by the many other “comoving” light hadrons produced in such a collision. To summarize the last point, if charmonium + light hadron “comover” dissociation cross sections are small (case 1, top of Fig.1) and the background of direct charm production from the initial nucleons is understood, one may have a useful signature of QGP formation, but if the comover dissociation cross sections are large (case 2, bottom of Fig.1) one must distinguish a QGP-reduced charmonium production amplitude from subsequent dissociative scattering, and the interpretation of this signal will therefore be ambiguous.

Thus it is of great relevance to the interpretation of RHIC physics to establish the approximate size of these low energy $`c\overline{c}`$ \+ light hadron cross sections.

## 2 Experiment, or what passes for it

Unfortunately we have no charmonium beams or targets, so the experimental cross sections must be inferred indirectly and are poorly known. The earliest estimates of lower energy charmonium hadronic cross sections came from $`J/\psi `$ photoproduction experiments in the mid 1970s, which were interpreted in terms of a $`J/\psi +N`$ cross section given additional theoretical assumptions. Early Fermilab and SLAC photoproduction experiments gave rough estimates of $`1`$ mb for $`\sigma _{J/\psi +N}`$, assuming vector dominance, for photon energies from $`E_\gamma 13`$ to $`200`$ GeV . A subsequent SLAC photoproduction experiment in 1977 used the $`A`$ dependence of $`J/\psi `$ absorption to infer a rather larger cross section of $`\sigma _{J/\psi +N}=3.5(0.8)`$ mb at $`E_\gamma 17`$ GeV ($`\sqrt{s}6`$ GeV) . The vector dominance hypothesis may have lead to an underestimate of the cross section in the earlier references .

In heavy ion collisions these cross sections may be estimated from the ratio of lepton pairs produced in the $`J/\psi `$ peak to “background” Drell-Yan pairs nearby in energy. Since the $`J/\psi `$ must reach the exterior of the nuclear target to decay into a sharp mass peak, this ratio gives us an estimate of the absorption cross section through a classical survival probability formula,

$$\sigma (J/\psi \mu ^+\mu ^{})/\sigma (DrellYan\mu ^+\mu ^{})$$

$$=\mathrm{exp}(\rho \sigma _{J/\psi +N}^{abs}L)$$
(1)
where $`\rho `$ is the mean nucleon density and $`L`$ is the estimated mean path length in the experimental nuclear system. A “naive” interpretation of the $`J/\psi `$ production data from collisions of various nuclear species using this formula gives $`\sigma _{J/\psi +N}^{abs}6`$ mb at $`\sqrt{s}10`$ GeV , with a numerically similar result for the $`\psi ^{}`$.

Of course one may raise many questions about the validity of this simple estimate, such as the importance of shadowing in Drell-Yan, the neglect of $`J/\psi `$ scattering by other light hadrons formed in the collision (such as $`\pi `$ and $`\rho `$), and the assumption of a single, constant $`J/\psi +N`$ cross section in all circumstances.

Recently, concerns have been expressed that the $`J/\psi `$ and $`\psi ^{}`$ wavefunctions have not had sufficient time to form within the nucleus in these collisions, so experiment may instead be measuring the cross section for a small initial $`c\overline{c}`$ “premeson” on a nucleon. One can increase the time spent in the interior of the nuclear system by selecting small and even negative $`x_F`$ events, as has been done by E866 at Fermilab. As discussed by He, Hüfner and Kopeliovich , this leads one to infer $`\sigma _{J/\psi +N}^{abs}=2.8(0.3)`$ mb and $`\sigma _{\psi ^{}+N}^{abs}=10.5(3.6)`$ mb respectively, also at $`\sqrt{s}10`$ GeV. This is rather more satisfying to people who have an intuitive notion that the larger $`\psi ^{}`$ should have a larger reaction cross section. Actually the connection between cross section and physical extent is less direct (compare $`KN`$ and $`\overline{K}N`$), and in any case the relative proximity of inelastic thresholds alone would suggest a larger $`\psi ^{}`$ cross section. These experiments also indicate a preference for dissociation over elastic cross sections in this energy region by roughly a factor of 30 .

## 3  Theory: Introduction and pQCD

To quote B.Müller, “…the state of the theory of interactions between $`J/\psi `$ and light hadrons is embarrassing \[\]. Only three serious calculations exist (after more than 10 years of intense discussions about this issue!) and their results differ by at least two orders of magnitude in the relevant energy range \[\]. There is a lot to do for those who would like to make a serious contribution to an important topic.” .

The theoretical situation has improved in the subsequent year, at least in terms of the quantity of calculations if not in the understanding of the scattering mechanism. A list of $`c\overline{c}`$ \+ light hadron cross section calculations is given in Table I.

The most cited work, albeit furthest in its predictions from a low-energy “theoretical mean”, is the pQCD calculation of Kharzeev and Satz . This reference is basically a restatement of the color-dipole scattering model developed in the late 1970s by Peskin and Bhanot  for scattering of light hadrons by Coulombic bound states of very massive quarks. According to Peskin, the criterion for validity of this approach is “…not met even for the $`b\overline{b}`$ system.” , so there may be large systematic errors at the $`c\overline{c}`$ mass scale. This approach certainly makes marginal approximations for charmonium, such as the use of Coulombic wavefunctions (which are far from accurate for $`c\overline{c}`$) and the introduction of a $`Q\overline{Q}`$ binding energy (which is hard to interpret for charmonium, and is taken to be $`2M_DM_{J/\psi }`$ by Kharzeev and Satz). These color-dipole scattering formulas also implicitly assume that charmonia are small relative to the natural QCD length scale. Since potential models actually find rms $`c\overline{c}`$ separations of about 0.4 fm for the $`J/\psi `$, 0.6 fm for the $`\chi _c`$ states and 0.8 fm for the $`\psi ^{}`$ , this is also a dubious assumption.

Although this approach has problems with justification for $`c\overline{c}`$, the predictions are nonetheless interesting as estimates of the scale of these cross sections assuming a color-dipole scattering mechanism, and the approach could presumably be extended to lower energies by generalizing the wavefunctions and interaction. The formula for the $`J/\psi +N`$ cross section quoted by Kharzeev and Satz  is

$$\sigma _{J/\psi +N}=2.5\text{ mb}\left(1\frac{\lambda _0}{\lambda }\right)^{6.5},$$
(2)
where $`\lambda =(sM_{J/\psi }^2M_N^2)/2M_{J/\psi }`$, the constant $`\lambda _0`$ is defined to be $`\mathrm{`}\mathrm{`}M_N+ϵ_0\mathrm{"}`$ according to the text following Eq.(24) of  (we assume the equality), and the “binding energy” $`ϵ_0`$ is set to $`2M_DM_{J/\psi }`$. The result is shown in Fig.3, together with the single lower-energy SLAC experimental point .

Evidently the Kharzeev-Satz cross section is smaller than this SLAC point (which was an inferred cross section and certainly needs confirmation) by about two orders of magnitude, and falls precipitously as $`\sqrt{s}`$ is decreased. This approach actually has no direct information about physical thresholds, so Kharzeev et al. typically leave their curves “dangling” just below $`\sqrt{s}=5`$ GeV. (See Fig.2 of  for example.) If we actually plot their formula (3.1) for $`J/\psi +N`$ at low energy (Fig.3), we find the unphysical prediction of a zero cross section near 4.5 GeV, whereas the physical $`\mathrm{\Lambda }_c+\overline{D}`$ threshold is at 4.15 GeV. Obviously the method is inapplicable at low energies, which is unfortunate because this regime is of greatest interest for RHIC studies of QGP production. Cross section calculations using this approach with more realistic $`c\overline{c}`$ wavefunctions are currently being carried out by Kopeliovich et al. .

## 4  Theory: meson exchange

Recently, several calculations of charmonium + light hadron cross sections have been reported assuming $`t`$-channel charmed meson exchange. Of course this picture is also problematic, since the range of the exchanged charmed meson would be only about $`1/M_D0.1`$ fm, and the assumption of nonoverlapping hadrons at this separation is clearly invalid. (This is the Isgur-Maltman  argument as to why vector meson exchange is unjustifed as the source of the short ranged NN core interaction.) Nonetheless it is again interesting to see what scale of cross section is predicted by this type of model, since it might at least incorporate the correct scales and degrees of freedom, and it assumes a different scattering mechanism from the pQCD color-dipole model advocated by Kharzeev and Satz.

The first such meson exchange calculation, due to Matinyan and Müller , considered $`t`$-channel $`D`$ exchange as the mechanism for the reactions $`J/\psi +\pi D^{}\overline{D}+h.c.`$ and $`J/\psi +\rho D\overline{D}`$. The cross sections found by this reference are shown in Fig.4 below. Note that 500 MeV above threshold these cross sections lie in the 0.5 to 1 mb range. Similarly, a recent $`D`$ and $`D^{}`$ meson exchange calculation of $`J/\psi +N\mathrm{\Lambda }_c+\overline{D}`$ by Sibirtsev, Tsushima and Thomas  finds a peak cross section of about 2 mb near $`\sqrt{s}=4.6`$ GeV. (An earlier calculation of $`J/\psi +N\mathrm{\Lambda }_c+\overline{D}`$ by Haglin , assuming $`D`$ exchange but no hadronic form factors, found a somewhat larger peak of about 7 mb near $`\sqrt{s}=4.3`$ GeV.) In comparison, Kharzeev and Satz predict nanobarn-scale $`J/\psi +N`$ cross sections 500 MeV above threshold, six orders of magnitude smaller! The scatter of theoretical predictions for these processes is remarkably wide.

This discrepancy between theoretical cross sections may appear discouraging. One can instead consider it an opportunity to learn something important about QCD: the predictions differ because they come from very different assumptions regarding the scattering mechanism, and as they are so far apart, for once we have a clear possibility of distinguishing between different hadron-hadron scattering models.

Meson exchange calculations of $`c\overline{c}`$ \+ $`q\overline{q}`$ dissociation cross sections have since been reported by Lin and Ko  and Haglin and Gale . The models differ considerably in detail, due to different choices for the diagrams included, the effective meson Lagrangian used, and the coupling constants and form factors assumed. With pointlike “hard” vertices, $`J/\psi +\pi `$ cross sections in the 10s of mb not far above threshold are typical. After introducing plausible hadronic form factors, these are usually reduced to peak values of $`1`$-$`10`$ mb near $`\sqrt{s}4.0`$-$`4.2`$ GeV. Unfortunately there is considerable “guessing” regarding hadronic couplings constants, which may be unnecessary because these can be calculated using well established quark model techniques, for example the 3P0 meson decay model. Similarly, hadronic form factors can be derived from the 3P0 model, presumably with sufficient accuracy for these purposes.

## 5  Theory: constituent interchange

We advocate an approach which uses nonrelativistic quark model wavefunctions and calculates these cross sections assuming a constituent interchange scattering mechanism, driven by the Born-order matrix element of the standard quark model Hamiltonian. This technique, which has no free parameters once quark model wavefunctions and the interquark Hamiltonian are specified, has been shown to compare reasonably well with experimental low energy hadron-hadron scattering data near threshold in a wide range of annihilation-free reactions . In meson-meson scattering there are four distinct diagrams (see Fig.6), each of which has an associated overlap integral of the nonrelativistic quark model external meson wavefunctions convolved with the interquark Hamiltonian. Constituent interchange is forced at Born-order because $`H_I\lambda ^a\lambda ^a`$ changes each initial color-singlet $`q\overline{q}`$ meson into a color octet, which has overlap with color-singlet final state mesons only after quark line interchange. The Feynman rules for these diagrams were given by Barnes and Swanson .

The first quark model calculation using this approach was reported by Martins, Blaschke and Quack , who considered the reactions $`J/\psi +\pi D^{}\overline{D}+h.c.`$ and $`D^{}\overline{D}^{}`$ (The amplitude for $`J/\psi +\pi D\overline{D}`$ is zero in the nonrelativistic quark model without spin-orbit forces, and has been found to be quite weak in a relativized calculation .) Martins et al. found that these exclusive final states have numerically rather similar cross sections (except for their different thresholds), and give a maximum total cross section of about 7 mb at $`\sqrt{s}4.1`$ GeV. A constituent interchange calculation of $`J/\psi +N`$ and $`\psi ^{}+N`$ cross sections using a simplified quark+diquark model of the nucleon  also found several-mb peak cross sections not far above threshold.

Our collaboration recently carried out quite similar constituent interchange quark model calculations (Wong et al. ). We used numerically determined Coulomb plus linear plus hyperfine quark potential model wavefunctions, and evaluated the Born-order meson-meson scattering amplitude, which is the matrix element of the interquark Hamiltonian between scattering states with quark interchange (Fig.6). We included smeared OGE spin-spin, OGE color Cou- lomb and linear confining interactions, with parameters chosen to give a good fit to the $`q\overline{q}`$ and $`c\overline{c}`$ meson spectra. We find a somewhat smaller cross section (peaking at about 1 mb at $`\sqrt{s}4.0`$ GeV) for the sum of these inelastic $`J/\psi +\pi `$ reactions. The difference between our work and that of Blaschke et al. lies mainly in the treatment of the confining interaction; for simplicity, Blaschke et al. treated confinement as a color-independent Gaussian potential that acts only between quark and antiquark (hence they include only diagrams C1 and C2), whereas we used the standard $`\lambda ^a\lambda ^a`$ linear confining potential between all pairs of constituents.

We find destructive interference between the C and T diagrams, leading to a much reduced total cross section. In  we also treated $`\psi ^{}+\pi `$ scattering, which involved a simple change to a 2S $`c\overline{c}`$ wavefunction and a change of phase space, and found a rather large, ca. 5 mb maximum cross section (see Fig.7). Our $`c\overline{c}`$ \+ $`q\overline{q}`$ cross sections have their strongest support just a few hundred MeV in $`\sqrt{s}`$ above threshold, since the overlap integrals are damped by the tails of the wavefunctions at higher energies.

We also consider $`J/\psi +\rho `$ (Fig.8) and $`\psi ^{}+\rho `$, which are predicted to be much larger near threshold for the simple reason that they are exothermic; there is a $`1/v_{AB}`$ divergence in these cross sections as we approach threshold.

Although our scattering amplitudes and cross sections are evaluated numerically, it is interesting that a simple two-parameter function gives a useful approximation to our (single channel) cross sections. This function is

$$\sigma (s)=\sigma _{max}(ϵ/ϵ_{max})^pe^{p(1ϵ/ϵ_{max})},$$
(3)
where $`ϵ=\sqrt{s}M_CM_D`$, and $`\sigma _{max}`$ (mb) is the maximum value of the cross section, at $`ϵ_{max}`$ (MeV). The threshold exponent $`p`$ is fixed by the angular quantum numbers of the hadrons, and is $`\pm 1/2+L_{min}^{CD}`$ (for endothermic/exothermic), where $`L_{min}^{CD}`$ is the lowest angular momentum allowed for the final meson pair $`CD`$. As an example, in Fig.9 we show our numerical results for the reaction $`\eta _c+\pi ^+D^+\overline{D}^o`$ and a fit using the function (5.1) with $`p=1/2`$, as appropriate for an S-wave-allowed final state. The masses and parameters assumed for this example were $`M_{\pi ^+}=0.140`$ GeV, $`M_{\eta _c}=2.98`$ GeV, $`M_{D^+}=1.869`$ GeV, $`M_{D^o}=1.865`$ GeV, $`\alpha _s=0.6`$, $`b`$ (string tension) $`=0.18`$ GeV2, $`m_{u,d}=0.33`$ GeV, $`m_c=1.6`$ GeV, and the OGE contact hyperfine smearing distance was $`0.25`$ fm. These are all reasonably well established nonrelativistic quark model parameters. The Schrödinger equation was solved with these parameters to generate numerical wavefunctions, which were then used in a nine-dimensional Monte Carlo integration to evaluate the scattering amplitudes in the CM frame. The Monte Carlo amplitude evaluation in Fig.9 used $`4M`$ points for each diagram at each energy and each final $`D^+(\widehat{\mathrm{\Omega }})`$ direction; amplitudes along three final directions were evaluated, which were then projected into S-, P- and D-moments to separate the different partial waves. Cancellation of diagram sums in certain channels as well as evaluation of known exact results with SHO wavefunctions provided nontrivial checks of the numerical work.

## 6 Conculsions

We have reviewed the recent theoretical predictions and experimental status of the cross sections of $`c\overline{c}`$ on light hadrons, which is of great interest for the interpretation of heavy ion collisions. There are three scattering mechanisms currently being investigated, which are color-di- pole pQCD, $`t`$-channel meson exchange, and constituent interchange. The pQCD approach gives very small cross sections at low energies, whereas the meson exchange and constituent interchange models both predict peak cross sections near thres- hold of $`1`$-$`10`$ mb. If it is possible to establish these cross sections experimentally, we may achieve a much better understanding of the mechanisms of low energy hadron-hadron scattering.

## 7 Acknowledgements

We would like to acknowledge the kind invitation and support of the organisers of “Heavy Quark Physics 5”, which made it possible to attend the meeting and present these results. This research was supported in part by the DOE Division of Nuclear Physics, at ORNL, managed by UT-Bat- telle, LLC, under Contract No. DE-AC05-00OR 22725, Research Corp, and by the Deutsche For- schungsgemeinschaft DFG under contract Bo 56/ 153-1. We would also like to thank D.B.Blaschke, B.Kopeliovich, M.Peskin, A.Sibirtsev and S.Sor- ensen for useful discussions and information.