text:
# Stability of Toomre-Hayashi Disk

## 1 INTRODUCTION

Toomre (1982) and Hayashi, Narita, & Miyama (1982) found an analytic model for a self-gravitating gas cloud. The model describes a disk in which the self-gravity balances with the centrifugal force and gas pressure. The rotation velocity, $`v_\phi `$, and the sound speed, $`c_\mathrm{s}`$, are assumed to be constant. The density distribution of the model is expressed as

$$\rho (r,\theta )=\frac{\gamma ^2c_\mathrm{s}^2}{2\pi Gr^2}\left[\frac{\mathrm{cosh}\zeta }{\mathrm{cosh}(\gamma \zeta )}\right]^2,$$
(1)
where

$$\zeta =\mathrm{ln}\left(\frac{1+\mathrm{cos}\theta }{\mathrm{sin}\theta }\right),$$
(2)
and

$$\gamma =\mathrm{\hspace{0.33em}1}+\frac{v_\phi ^2}{2c_\mathrm{s}^2},$$
(3)
in the spherical polar coordinates, ($`r`$, $`\theta `$, $`\phi `$).

This simple model is useful for studying dynamics of self-gravitating gas disk. Earlier researchers (see, e.g., Kalnajs, 1971; Lemos, Kalnajs, & Lynden-Bell, 1991; Syer & Tremaine, 1996; Evans & Read, 1998a, b; Goodman & Evans, 1999; Shu et al., 2000) studied further simplified model, i.e., the razor-thin power law disk model, which is equivalent to the vertically averaged Toomre-Hayashi model. Besides of its simplicity, it is a good approximation to a centrifugally supported disk formed by collapse of a rotating gas cloud (see, e.g., Matsumoto, Hanawa, & Nakamura, 1997; Saigo & Hanawa, 1998). The centrifugally supported disk will evolve into single or binary stars associated with gas disks. The evolution may be driven by instabilities of the disk. Since the Toomre-Hayashi model is an exact equilibrium solution, it is ideal for linear stability analysis; we need not worry about inaccuracy of the unperturbed state.

In this paper we concentrate on the simplest mode, the axisymmetric one. We will show that the Toomre-Hayashi model is unstable for any $`\gamma `$, i.e., any $`v_\phi /c_\mathrm{s}`$. This result is different from prediction by Hayashi et al. (1982) and earlier analyses by Lemos et al. (1991) and Shu et al. (2000). While two-dimensional perturbation is taken into account in our analysis, only radial perturbation was considered in earlier analysis. We believe that our analysis will provide insights for stability of self-gravitating disk.

The rest of our paper has the following structure. In §2 we derive energy principle for axisymmetric perturbation around the Toomre-Hayashi disk. In §3 we obtain the marginally stable mode for the Toomre-Hayashi disk. In §4 we show that it is unstable for any $`\gamma `$ and discuss implications of our stability analysis. Main conclusions are summarized in §5.

## 2 PERTURBATION EQUATION AND ENERGY PRINCIPLE

We consider a small perturbation around the Toomre-Hayashi disk in the polar coordinates. In this paper we use the displacement vector, $`𝝃=(\xi _r,\xi _\theta )`$, in the $`\phi =\mathrm{const}.`$ plane to denote a perturbation. From the equation of continuity we can evaluate the change in the density as

$$\rho _1=\frac{1}{r^2}\frac{}{r}(r^2\rho _0\xi _r)\frac{1}{r\mathrm{sin}\theta }\frac{}{\theta }(\mathrm{sin}\theta \rho _0\xi _\theta ),$$
(4)
where $`\rho _0`$ denotes the density in the equilibrium and is given by Equation (1). From the angular momentum conservation we can evaluate the change in the rotation velocity as

$$v_{\phi 1}=\frac{v_{\phi _0}}{r}\left(\xi _r+\frac{\xi _\theta }{\mathrm{tan}\theta }\right),$$
(5)
where $`v_{\phi 0}`$ denotes the rotation velocity in equilibrium. The perturbation equation for the equation of motion is expressed as

$$\frac{^2\xi _r}{t^2}\frac{2v_{\phi 0}v_{\phi 1}}{r}+c_\mathrm{s}^2\frac{}{r}\left(\frac{\rho _1}{\rho _0}\right)+\frac{\varphi _1}{r}=\mathrm{\hspace{0.33em}0},$$
(6)
and

$$\frac{^2\xi _\theta }{t^2}\frac{2v_\phi v_{\phi 1}}{r\mathrm{tan}\theta }+\frac{c_\mathrm{s}^2}{r}\frac{}{\theta }\left(\frac{\rho _1}{\rho _0}\right)+\frac{1}{r}\frac{\varphi _1}{\theta }=\mathrm{\hspace{0.33em}0},$$
(7)
for the $`r`$\- and $`\theta `$-directions, respectively. The change in the gravitational potential, $`\varphi _1`$, is related to that in the density by

$$\left[\frac{1}{r^2}\frac{}{r}\left(r^2\frac{}{r}\right)+\frac{1}{r^2\mathrm{sin}\theta }\frac{}{\theta }\left(\mathrm{sin}\theta \frac{}{\theta }\right)\right]\varphi _1=\mathrm{\hspace{0.33em}4}\pi G\rho _1.$$
(8)
Before detailed manipulation we consider the characteristic of the Tooomre-Hayashi model. Since the Toomre-Hayashi model has neither inner nor outer boundary, its global analysis is not straight forward (see Goodman & Evans, 1999, for mathematical discussions on the stability of scale-free disks).

The Tooomre-Hayashi disk is scale free and has no characteristic length scale. It has neither characteristic timescale; both the rotation period and sound crossing time are proportional to the radial distance, $`r`$. Since the disk has no characteristic radius, inner and outer parts of a disk are equally either stable or unstable. This means that the disk has no discrete eigenmode which grows exponentially in time (see also Goodman & Evans, 1999).

If the Toomre-Hayashi disk is unstable, the perturbation grows faster in an inner part and propagate outwards since the dynamical timescale is shorter at a smaller radius. To show this we derive the energy principle,

$$𝒯+𝒲+𝒮=\mathrm{\hspace{0.33em}0}$$
(9)
where

$`𝒯`$ $`=`$ $`{\displaystyle _a^b}{\displaystyle _0^\pi }\rho _0\left|\xi _r^{}{\displaystyle \frac{^2\xi _r}{t^2}}+\xi _\theta ^{}{\displaystyle \frac{^2\xi _\theta }{t^2}}\right|r^2\mathrm{sin}\theta drd\theta ,`$ (10)
$`𝒲`$ $`=`$ $`{\displaystyle _a^b}{\displaystyle _0^\pi }\left\{\rho _0\left[2|v_{\phi 1}|^2+c_\mathrm{s}^2\left|{\displaystyle \frac{\rho _1}{\rho _0}}\right|^2\right]+\rho _1^{}\varphi _1\right\}r^2\mathrm{sin}\theta drd\theta ,`$ (11)
$`𝒮`$ $`=`$ $`\left[{\displaystyle _0^\pi }r^2\mathrm{sin}\theta \xi _r^{}(c_\mathrm{s}^2\rho _1+\rho _0\varphi _1)𝑑\theta \right]_a^b.`$ (12)
Equation (9) is derived from equations (4) through (8) by multiplying the complex conjugate of $`\rho _1`$, $`\xi _r`$, $`\xi _\theta `$, and $`\varphi _1`$. We can set the inner and outer radii, $`a`$ and $`b`$, arbitrarily. Equation (9) denotes the energy conservation in the second order in the interval of $`arb`$.

The term, $`𝒯`$, denotes the product of acceleration and displacement. If it is positive, the perturbation grows and the system is unstable. The term, $`𝒲`$, denotes the excess of rotation energy, thermal energy and gravitational energy. We consider the condition that the term, $`𝒮`$, vanishes. This condition is equivalent to

$`\xi _r(r,\theta )`$ $`=`$ $`\xi _r(\theta )r^{ik+1/2},`$ (13)
$`\xi _\theta (r,\theta )`$ $`=`$ $`\xi _\theta (\theta )r^{ik+1/2},`$ (14)
$`\rho _1(r,\theta )`$ $`=`$ $`\rho _1(\theta )r^{ik5/2},`$ (15)
$`\varphi _1(r,\theta )`$ $`=`$ $`\varphi _1(\theta )r^{ik1/2},`$ (16)
where $`k`$ denotes the wavenumber in the logarithmic scale (= $`/\mathrm{ln}r`$). When $`k\mathrm{ln}(b/a)=\mathrm{\hspace{0.17em}2}\pi `$, the third integral vanishes in equation (9). We call this the periodic boundary condition in the following. Our periodic boundary condition is equivalent to the surface density perturbation assumed by Kalnajs (1971) and others.

When we apply the periodic boundary condition, the term $`𝒮`$ is proportional to $`\mathrm{ln}(b/a)`$ and scale free. Also the term $`𝒯`$ should be scale free and the growth timescale should be proportional to the radius, $`/t\mathrm{\hspace{0.17em}1}/r`$. This means that there exists no exponentially growing mode satisfying the periodic boundary condition. We seek, instead, a marginally stable mode satisfying the periodic boundary condition.

Substituting the boundary condition into the perturbation equation we obtain

$$\left(ik+\frac{1}{2}\right)\xi _r+\frac{1}{\rho _0}\frac{\rho _0}{\theta }\xi _\theta +\frac{1}{\mathrm{sin}\theta }\frac{}{\theta }(\mathrm{sin}\theta \xi _\theta )+\left(\frac{\rho _1}{\rho _0}\right)=\mathrm{\hspace{0.33em}0},$$
(18)
$$2v_{\phi 0}v_{\phi 1}+\left(ik\frac{1}{2}\right)c_\mathrm{s}^2\left(\frac{\rho _1}{\rho _0}\right)+\left(ik\frac{1}{2}\right)\varphi _1=\mathrm{\hspace{0.33em}0},$$
(19)
$$\frac{2v_{\phi 0}v_{\phi 1}}{\mathrm{tan}\theta }+c_\mathrm{s}^2\frac{}{\theta }\left(\frac{\rho _1}{\rho _0}\right)+\frac{\varphi _1}{\theta }=\mathrm{\hspace{0.33em}0},$$
(20)
$$v_{\phi 1}=v_{\phi 0}\left(\xi _r+\frac{\xi _\theta }{\mathrm{tan}\theta }\right),$$
(21)
$$\left[\left(k^2+\frac{1}{4}\right)+\frac{1}{\mathrm{sin}\theta }\frac{}{\theta }\left(\mathrm{sin}\theta \frac{}{\theta }\right)\right]\varphi _1=\mathrm{\hspace{0.33em}4}\pi G\rho _1,$$
(22)
where

$$\rho _0=\frac{\gamma ^2c_\mathrm{s}^2}{2\pi G}\left[\frac{\mathrm{cosh}\zeta }{\mathrm{cosh}(\gamma \zeta )}\right]^2$$
(23)
and

$$v_{\phi 0}=c_\mathrm{s}\sqrt{2(\gamma \mathrm{\hspace{0.17em}1})}.$$
(24)
We solve equations (18) through (22) as an eigenvalue problem for $`k`$.

## 3 MARGINALLY STABLE MODE

### 3.1 Singular Isothermal Sphere

When $`\gamma =\mathrm{\hspace{0.17em}1}`$, the Toomre-Hayashi disk is spherically symmetric and reduces to the singular isothermal sphere. The density distribution is expressed as

$$\rho _0=\frac{c_\mathrm{s}^2}{2\pi Gr^2}.$$
(25)
Then the perturbation should be expressed by the spherical harmonics, $`Y_{\mathrm{}}^m(\theta ,\phi )`$. We concentrate on the axisymmetric perturbation for simplicity. Then the density perturbation is expressed as

$$\rho _1(\theta )=AP_{\mathrm{}}(\mathrm{cos}\theta ),$$
(26)
where $`A`$ and $`P_{\mathrm{}}`$ denotes an arbitrary scale factor and the Legendre’s polynomial, respectively. Then the change in the potential is given by

$$\varphi _1=\mathrm{\hspace{0.17em}4}\pi G\rho _1\left[k^2+\frac{1}{4}+\mathrm{}(\mathrm{}+\mathrm{\hspace{0.17em}1})\right]^1.$$
(27)
Substituting equation (27) into equation (20) we obtain

$$k^2=\frac{7}{4}\mathrm{}(\mathrm{}+\mathrm{\hspace{0.17em}1}).$$
(28)
From equation (28) we conclude that the marginally stable mode exists only for $`\mathrm{}=\mathrm{\hspace{0.17em}0}`$ and the critical wavenumber is $`k=\sqrt{7}/2`$.

The critical wavenumber coincides with that of the equilibrium isothermal sphere. As noted in the textbook of Chandrasekhar (1939), the density distribution is approximately given by

$$\rho =\frac{c_\mathrm{s}^2}{2\pi Gr^2}\mathrm{exp}\left\{\frac{A}{\sqrt{r}}\mathrm{cos}\left[\frac{\sqrt{7}}{2}\mathrm{log}\left(\frac{r\sqrt{4\pi G\rho _\mathrm{c}}}{c_\mathrm{c}}\right)\delta \right]\right\}.$$
(29)
in the region of $`rc_\mathrm{s}/\sqrt{4\pi G\rho _\mathrm{c}}`$, where $`\rho _\mathrm{c}`$ denotes the central density. This equilibrium isothermal sphere has almost the same density distribution as that of the singular isothermal sphere in the region of $`rc_\mathrm{s}/\sqrt{4\pi G\rho _\mathrm{c}}`$. The small difference between them can be regarded as a neutral perturbation on the singular isothermal sphere. Note the wavenumber, $`\sqrt{7}/2`$, in Equation (29).

The singular isothermal sphere is stable against non-spherical perturbation of $`\mathrm{}\mathrm{\hspace{0.17em}0}`$. Equation (28) leaves formally a possibility that the singular isothermal disk is unstable for any $`k`$. This possibility is, however, very unlikely. When $`k`$ is very large, the wavelength is so short to destabilizes the gas sphere. Thus we can exclude the possibility that the singular isothermal sphere is unstable.

### 3.2 Toomre-Hayashi Disk

In this subsection we solve the perturbation equation to obtain the critical wavenumber $`k`$ as a function of $`\gamma `$. First we obtain

$$\frac{v_{\phi 1}}{\theta }+\mathrm{tan}\theta \left(ik\frac{1}{2}\right)v_{\phi 1}=\mathrm{\hspace{0.33em}0},$$
(30)
from equations (19) and (20). Equation (30) has the solution,

$$v_{\phi 1}(\mathrm{sin}\theta )^{ik\mathrm{\hspace{0.17em}1}/2}.$$
(31)
To avoid the divergence of $`v_{\phi 1}`$ at $`\theta =\mathrm{\hspace{0.17em}0}`$, we seek modes of $`v_{\phi 1}=\mathrm{\hspace{0.17em}0}`$.

Substituting $`v_{\phi 1}=\mathrm{\hspace{0.17em}0}`$ into equations (18) through (22), we obtain

$$\xi _r=\frac{\xi _\theta }{\mathrm{tan}\theta }$$
(32)
$$c_\mathrm{s}^2\frac{\rho _1}{\rho _0}=\varphi _1,$$
(33)
$$\left(ik+\frac{1}{2}\right)\frac{\xi _\theta }{\mathrm{tan}\theta }+\frac{1}{\rho _0}\frac{\rho _0}{\theta }\xi _\theta +\frac{1}{\mathrm{sin}\theta }\frac{}{\theta }(\mathrm{sin}\theta \xi _\theta )\left(\frac{\varphi _1}{c_\mathrm{s}^2}\right)=\mathrm{\hspace{0.33em}0},$$
(34)
and

$$\left[\left(k^2+\frac{1}{4}\right)+\frac{1}{\mathrm{sin}\theta }\frac{}{\theta }\left(\mathrm{sin}\theta \frac{}{\theta }\right)\right]\varphi _1=\frac{4\pi G\rho _0}{c_\mathrm{s}^2}\varphi _1.$$
(35)
We solve Equation (35) with the boundary conditions $`\varphi _1/\theta =\mathrm{\hspace{0.17em}0}`$ at $`\theta `$ = 0 and $`\pi `$ as an eigenvalue problem. Once $`\varphi _1`$ is given, other perturbations, $`\xi _r`$, $`\xi _\theta `$, and $`\rho _1`$ are derived from Equations (32) through (34). For a given set of $`\gamma `$ and $`k`$, we integrated Equation (35) with the Runge-Kutta method. By try and error, we found pairs of $`\gamma `$ and $`k`$ for which Equation (35) has a solution satisfying the boundary condition at $`\theta `$ = 0 and $`\pi `$.

Figure 1 denotes the eigenvalue, $`k`$, as a function of $`\gamma `$. The solid curve connects the crosses which denotes the numerical obtained eigen values. The dashed curve does $`k=\gamma `$ for comparison. It is a good approximation to the numerically obtained eigenvalues for $`\gamma \mathrm{\hspace{0.17em}1}`$.

Figure 2 shows $`\rho _0/\rho _1`$ as a function of $`\theta `$ for the marginally stable modes. Each curve is labeled with $`\gamma `$. The relative density perturbation, $`\rho _1/\rho _0`$, is symmetric with respect to $`\theta =\pi /2`$ and has no zero point. When $`\gamma `$ is very large, the density perturbation is restricted near $`\theta =\pi /2`$. Note that Figure 2 shows not $`\rho _1`$ but $`\rho _1/\rho _0`$. The density perturbation, $`\rho _1`$, is very small near the pole ($`\theta =\mathrm{\hspace{0.17em}0}`$ and $`\pi `$) even for $`\gamma =\mathrm{\hspace{0.17em}2}`$, since $`\rho _0`$ is very small there.

After having completing numerical computations, we found an asymptotic solution for Equation (35). To derive the asymptotic solution, we rewrite Equation (35) into

$$\left[\left(k^2+\frac{1}{4}\right)+\mathrm{cosh}^2\zeta \frac{^2}{\zeta ^2}\right]\varphi _1=\mathrm{\hspace{0.17em}2}\gamma ^2\left[\frac{\mathrm{cosh}\zeta }{\mathrm{cosh}(\gamma \zeta )}\right]^2\varphi _1.$$
(36)
When $`\gamma \mathrm{\hspace{0.17em}1}`$, $`\varphi _1`$ has an appreciable amplitude only near the disk plane (i.e., $`\theta \pi /2`$ or equivalently $`\zeta \mathrm{\hspace{0.17em}0}`$). Thus, Equation (36) can be approximated as

$$\left[\left(k^2+\frac{1}{4}\right)+\frac{^2}{\zeta ^2}\right]\varphi _1=\frac{2\gamma ^2}{\mathrm{cosh}^2(\gamma \zeta )}\varphi _1.$$
(37)
Equation (37) has an analytic solution

$$\varphi _1=\frac{1}{\mathrm{cosh}(\gamma \zeta )},$$
(38)
for $`k^2=\gamma ^2\mathrm{\hspace{0.17em}1}/4`$. This denotes an asymptotic solution in the limit of $`\gamma \mathrm{\hspace{0.17em}1}`$.

Figures 3 and 4 show the marginally stable perturbation by cross section. Figure 3 shows the density distribution in the Toomre-Hayashi model of $`\gamma =\mathrm{\hspace{0.17em}10}`$. The solid curves are the contours of $`\mathrm{log}\rho =n/4`$ where $`n`$ is an integer. Only the sector near the disk plane is shown in Figure 3. Figure 4 is the same as Figure 3 but the marginally stable mode is superimposed on the Toomre-Hayashi model. As shown in Figure 4 the neutrally stable mode deforms the Toomre-Hayashi model into rings.

## 4 DISCUSSION

### 4.1 Instability Strip

We have shown in the previous section that the Toomre-Hayashi model has a marginally stable mode for any $`\gamma `$. This means that the Toomre-Hayashi model of any $`\gamma `$ is unstable against axisymmetric perturbation, since a marginally stable mode is adjacent to an unstable mode.

As shown by Hayashi et al. (1982), the Toomre-Hayashi model is unstable for total collapse when $`\gamma <\mathrm{\hspace{0.17em}2.441}`$. In other words, the Toomre-Hayashi model of $`\gamma =\mathrm{\hspace{0.17em}2.441}`$ has a marginally stable mode of $`k`$ = 0. This marginally stable point is also shown by an asterisk in Figure 1. It is also shown in the previous section that the Toomre-Hayashi model of $`\gamma =\mathrm{\hspace{0.17em}0}`$ (i.e., the singular isothermal sphere) is unstable when $`k<\sqrt{7}/2`$. Accordingly, we conclude that the Toomre-Hayashi model is unstable in the lower right side of the marginally stable curve in Figure 1.

The above conclusion apparently contradicts with the stability analyses by Lemos et al. (1991) and Shu et al. (2000) and prospect by Hayashi et al. (1982). All of them concluded that the Toomre-Hayashi model of $`2.5\gamma \mathrm{\hspace{0.17em}6}`$ is stable against axisymmetric perturbation.

The dash-dotted curves of Figure 1 denote the points of marginally stability obtained by Shu et al. (2000). They analyzed the stability of the Toomre-Hayashi disk under the thin disk approximation while taking account of magnetic fields. The magnetic field reduces the effective gravity and increases the sound speed (see, e.g., Li & Shu, 1997; Shu & Li, 1997; Nakamura & Hanawa, 1997). The effects can be fully taken into account by replacing the gravitational constant ($`G`$) and sound speed ($`c_\mathrm{s}`$) with the effective ones. Thus the stability of the magnetized disk is essentially the same as that of non-magnetized one. The condition for the marginal stability of the axisymmetric mode is denoted by

$$\gamma =\mathrm{\hspace{0.33em}1}+\frac{\left(k^2+\mathrm{\hspace{0.17em}1}/4\right)[K(k)\mathrm{\hspace{0.17em}1}]}{2\left(k^2+\mathrm{\hspace{0.17em}1}/4\right)K(k)},$$
(39)
where

$$K(k)=\frac{1}{2}\frac{\mathrm{\Gamma }[(1/2+ik)/2]\mathrm{\Gamma }[(1/2ik)/2]}{\mathrm{\Gamma }[(3/2+ik)/2]\mathrm{\Gamma }[(3/2ik)/2]},$$
(40)
denotes the Kalnajs function for $`m=\mathrm{\hspace{0.17em}0}`$. Lemos et al. (1991) obtained also the same dispersion relation although Equation (39) was not explicitly printed.

As shown in Figure 1, the disturbance is oscillatory and stable when $`\gamma `$ is medium. This discrepancy is due to incompleteness in their search for unstable modes. Shu et al. (2000) considered only radial perturbations neglecting vertical structure. However, vertical perturbation is dominant in the marginally stable mode that we found in this paper. Thus they could not find this marginally stable mode and adjacent unstable mode.

Equation (39) has an asymptotic form of

$$\gamma =k+\mathrm{\hspace{0.17em}2}+\frac{15}{8k}+𝒪(k^2),$$
(41)
in the limit of $`k\mathrm{\hspace{0.17em}3}`$. Note that the solid and dash-dotted curves are parallel in the region of large $`k`$ and hence for large $`\gamma `$. Change in the vertical structure has a significant effect even for a disk of a large $`\gamma `$. Comparison of our analysis and that of Shu et al. (2000) tells us importance of vertical perturbation in the instability of the Toomre-Hayashi model. The importance will be also true for stability of self-gravitating disk in general.

Hayashi et al. (1982) prospected on the basis of their numerical simulations that the Toomre-Hayashi is stable when $`\gamma `$ is medium. Their search was also incomplete since the the dynamic range in the radial direction is restricted. They could not find an unstable mode having a small $`k`$. Note that the wavenumber of the marginally stable mode is 3.102 and 5.058 for $`\gamma `$ = 3 and 5, respectively. This means that the radius increases by a factor of 7.58 and 3.46, respectively, at every node of the marginally stable mode. Thus large dynamic range is necessary to find an unstable mode from numerical simulations. We think that we need dynamic range covering at least two wavelengths.

As shown in the previous section, the change in the rotation velocity ($`\delta v_\phi `$) vanishes in the marginally stable mode. This means that the displacement is only in the vertical direction. The marginally stable mode changes hight of the disk but not the surface density.

### 4.2 Eigenmode

Although we have tried, we could not find the lower limit on $`k`$ for a given $`\gamma `$. In other words, we failed in finding the other side of the instability strip. This failure is not due to lack in our numerical technic.

One might think that the other marginally stable mode would stem from the edge of the total collapse, $`(k,\gamma )=(0,\mathrm{\hspace{0.17em}0.2441})`$. It, however, cannot stem. The marginally stable total collapse has radial dependence different from those assumed in Equations (13) through (16). The change in the density distribution is expressed as

$`\rho _1`$ $`=`$ $`{\displaystyle \frac{\rho _0}{\gamma }}`$ (42)
$`=`$ $`\left[{\displaystyle \frac{2}{\gamma }}\mathrm{\hspace{0.17em}2}\zeta \mathrm{tanh}(\gamma \zeta )\right]\rho _0`$ (43)
$``$ $`r^2,`$ (44)
for the marginally stable total collapse. This total collapse was called the breathing mode in Lemos et al. (1991).

The other marginally stable mode might have different radial dependence which have been not yet known. Usually we assume some temporal and spatial dependences when we search unstable modes. Often we assume that the perturbation grows exponentially. An eigenmode may not grow exponentially in time in the Toomre-Hayashi disk. Remember our argument based on the energy principle in §2 (see also Goodman & Evans, 1999, for another discussion on the stability of scale-free disks).

The expansion wave solution of Shu (1977) gives us a hint. The eigenmode may have a form of

$$\rho _1=r^2f(\frac{r}{c_\mathrm{s}t},\theta )$$
(45)
The density distribution is expressed as

$$\rho =\frac{c_\mathrm{s}^2}{2\pi Gr^2}g\left(\frac{r}{c_\mathrm{s}t}\right),$$
(46)
in the expansion wave solution. The difference from the singular isothermal sphere can be regarded as a perturbation propageting outwards while growing. As pointed out by Hunter (1977), the singular isothermal sphere can be regarded as a similarity solution. Similarly the Toomre-Hayashi model can be regarded as a similarity solution. An unstable mode can grow in proportion to the power of time when the unperturbed state is a similarity solution (see, e.g., Hanawa & Nakayama, 1997).

## 5 CONCLUSION

We discussed the stability of the Toomre-Hayashi model against axisymmetric perturbations in this paper. Our approach was conservative since we focused on the marginally stable mode. We did not obtain an unstable mode. Nevertheless we had some new findings. This is because we used least approximations and took account of the vertical structure exactly.

Although our analysis is restricted in the axisymmetric mode, the result will be helpful for more general analysis on the stability of self-gravitating disk. For later use we summarize the main results in the following.

1. The Toomre-Hayashi model is unstable against axisymmetric perturbation irrespectively of its model parameter, $`\gamma `$.
2. The wavelength of the unstable perturbation is proportional to the radius. The ratio of the wavelength to the radius is larger for smaller $`\gamma `$.
3. Change in the vertical height can induce instability even when the surface density does not change.
4. It is prospected that an unstable perturbation may grow in proportion of time in the power of time.
5. The Toomre-Hayashi model of $`\gamma =\mathrm{\hspace{0.17em}0}`$, i.e., the singular isothermal sphere is proven to be stable against a non-spherical perturbation.

The authors thank Tetsuro Konishi for telling us the asymptotic solution shown in §3. This research is financially supported in part by the Grant-in-Aid for Scientific Research on Priority Areas and that for Encouragement of Young Scientists of the Ministry of Education, Science, Sports and Culture of Japan (Nos. 10147105, 11134209, 1202103, and 12740123).


format: text

text:
# Acknowledgments

## Acknowledgments

This work was partially supported by CNPq and by PRONEX (Programa de Apoio a Núcleos de Excelência), BRAZIL. MVTM acknowledges the DF-UFPel for their warm hospitality. MBGD acknowledges enlightening discussions with Drs. Carlos Pajares and Alphonse Capella and the hospitality of the Departamento de Física de Partículas, U. Santiago de Compostela, where part of this work was accomplished. VPBG thanks FAPERGS and CNPq for support.

## Figure Captions

Fig. 1: The $`x_{IP}`$ dependence of the $`Q^2`$-slope at some typical values of $`\beta `$ for: (a) the pure CKMT model ; (b) the pQCD-inspired model . Kinematical contraint $`Q^2=Q^2(x)`$ from the MRST group  was used.

Fig. 2: The $`\beta `$ dependence of the $`Q^2`$-slope at some typical values of $`x_{IP}`$ for: (a) the pure CKMT model ; (b) the pQCD-inspired model . The kinematical contraint $`Q^2=Q^2(x)`$ from the MRST group  was used.

Fig. 3: The $`x_{IP}`$-slope versus $`Q^2`$ for the pQCD approach (BW)  and the CKMT model . The Donnachie-Landshoff intercept  is also depicted.


text:
# Experimental status of pionium at CERN

The experiment DIRAC  aims to measure the lifetime of the $`\pi ^+\pi ^{}`$ atom (pionium) in the ground state with $`10\%`$ precision. DIRAC is a magnetic double arm spectrometer and uses the high intensity 24 GeV/c proton beam of the CERN PS. Since this lifetime $`\tau `$ of order $`10^{15}`$ s is dictated by strong interaction at low energy, a precise measurement of $`\tau `$ enables to determine strong S-wave pion scattering lengths to $`5\%`$ . In Chiral Perturbation Theory (ChPT) pion scattering lengths are studied since many years, resulting in better and better predictions, now at an accuracy around $`2\%`$ . Such a measurement would submit the predictions of ChPT and hence the understanding of chiral symmetry breaking of QCD to a crucial test.

The experimental method to measure the pionium lifetime ($`3`$ fs) takes advantage of the Lorentz boost of relativistic pionia produced in high energy proton nucleus reactions at 24 GeV/c. After production in the target (e.g. 0.1 mm thin Ni foil) these relativistic ($`\gamma 15`$) atoms may either decay into 2 $`\pi ^0`$ or get excited or ionized in the target material. In the case of ionization or breakup, characteristic charged pion pairs, called “atomic pairs”, will emerge, exhibiting low relative momentum in their pair system ($`Q<3`$ MeV/c), small pair opening angle ($`\theta _+<3`$ mrad) and nearly identical energies in the lab system ($`E_+E_{}`$ at the 0.3% level). The task of the spectrometer is to identify charged pion pairs and to measure their relative momenta $`Q`$ with high resolution of $`\delta Q1`$ MeV/c. By these means, it is possible to determine the number of “atomic pairs” above background, arising from pion pairs in a free state (“free pairs”). For a given target material and thickness, the ratio of the number of observed “atomic pairs” to the total amount of produced pionia depends on $`\tau `$ in a unique way.

The experiment consists of coordinate detectors, a spectrometer magnet (bending power of 2.3 Tm) and two telescope arms, each equipped with drift chambers, scintillation hodoscopes, gas Cherenkov counters, preshower and muon detectors.

After tuning the primary proton beam as well as all detectors, the experiment started to accumulate data in 1999. The setup performance has been studied by analysing data collected with a platinum target. In order to track the relative momentum resolution, the invariant mass spectrum of p and $`\pi ^{}`$$`\mathrm{\Lambda }`$ decay products — was investigated. The position of the mass peak in Fig. 1 corresponds to $`\mathrm{\Lambda }`$ particles reconstructed in the DIRAC spectrometer. The width of the $`\mathrm{\Lambda }`$ peak, mainly given by the momentum resolution, is found to be $`0.43`$ MeV.

A preliminary search for “atomic pairs” in the 1999 platinum data sample is presented in Fig. 2. As mentioned above, “atomic pairs” are characterized by low relative momentum $`Q<3`$ MeV/c or equivalently $`F<3`$ 111Definition: $`F=\sqrt{(Q_L/\sigma _{Q_L})^2+(Q_X/\sigma _{Q_X})^2+(Q_Y/\sigma _{Q_Y})^2}`$ with resolutions $`\sigma _Q`$.. The distribution of “reals” $`N_{\mathrm{real}}(F)`$ is the sum of two distributions: “atomic pairs” (for $`F<4`$) and “free pairs”. In the region $`F>4`$ the “free pair” distribution $`N(F)`$ has been approximated by a fit function $`A(F)`$, based on “accidental pairs”. In Fig. 2 the difference between the experimental distribution $`N_{\mathrm{real}}(F)`$ and $`A(F)`$ is shown. The excess of $`\pi ^+\pi ^{}`$ pairs with $`F<2`$ amounts to $`160\pm 45`$, which is compatible with the expected number of “atomic pairs” $`240`$.


text:
# Gas Dynamics in Galaxy Mergers

## 1. Background

It has long been recognized that interstellar material plays an important role in galaxy interactions. Spitzer & Baade (1951) described fast collisions of galaxies as essentially hydrodynamic affairs. While emphasizing that stellar dynamics governs the outcome of galaxy encounters, Toomre & Toomre (1972) anticipated more than a decade of numerical work in suggesting that interactions ‘bring deep into a galaxy a fairly sudden supply of fresh fuel’.

The mechanics of this fueling process seem fairly clear. Negroponte & White (1983) found central concentrations of gas in simulated mergers of gas-rich disk galaxies, but the low resolution of their models left the physics uncertain. Noguchi (1987, 1988) and Combes, Dupraz, & Gerin (1990) showed that tidal encounters trigger bar formation in stellar disks and that gravitational torques exerted by these bars extract angular momentum from the gas. Hernquist (1989) reported rapid gas inflows in minor mergers where a disk galaxy swallows a small companion; in these experiments the gravitational torques driving the gas inward were generated by both the stellar disk and the companion. Barnes & Hernquist (1991, 1996) modeled major mergers between disk galaxies. This work extended the gravitational inflow picture to show how orbital decay could deliver $`5\times 10^9\mathrm{M}_{}`$ of gas within the central $`100\mathrm{pc}`$ of a merger remnant; such massive gas clouds are common in the central regions of infrared-selected galaxies, many of which are merger remnants (Sanders & Mirabel 1996). Hernquist & Barnes (1991) simulated the formation of a counterrotating gas disk in the aftermath of a major merger. Misaligned or counterrotating disks are found in both merger remnants (Schweizer 1982) and elliptical galaxies (e.g. Bender 1990), strengthening the hypotheses that merging disk galaxies form elliptical galaxies (Toomre & Toomre 1972) and that the centralized starbursts which power luminous infrared galaxies represent the formation of the cores of elliptical galaxies (e.g. Kormendy & Sanders 1992).

Existing simulations are still incomplete in one key respect: they don’t really capture the thermodynamics of the gas. Early studies used ‘sticky particles’ to mimic effects of dissipation. Hernquist and co-workers adopted ‘smoothed particle hydrodynamics’ (SPH) techniques which can include radiative heating and cooling. But in practice, radiative cooling is usually cut off below $`10^4\mathrm{K}`$ to prevent catastrophic instabilities, and most of the gas remains near this cut-off temperature throughout the simulation. Consequently, simulations including radiative processes are all but indistinguishable from those using an isothermal equation of state (Barnes & Hernquist 1996).

At this point there are two options. One is to adopt the isothermal gas model and survey a range of encounters. The other is to try to develop models which represent the thermodynamics of a multi-phase ISM. Below I discuss both.

## 2. Isothermal Encounter Survey

Isothermal simulations are relatively cheap, so an encounter survey is possible with modest computers. I chose four encounter geometries, and ran each once with a galactic mass ratio of 1:1 and once with a mass ratio of 3:1. Table 1 lists inclinations $`i`$ and pericentric arguments $`\omega `$ for the four encounter geometries; in the 3:1 encounters $`i_1`$ and $`\omega _1`$ refer to the more massive galaxy. These experiments used the same bulge/disk/halo galaxy models and close parabolic orbits as Barnes (1998, Ch. 4), but here each disk included a gaseous component amounting to $`12.5`$% of the disk mass. In simulation units with $`G=1`$, each galaxy in the 1:1 encounters, and each large galaxy in the 3:1 encounters, has total mass $`M_{\mathrm{bulge}}+M_{\mathrm{disk}}+M_{\mathrm{halo}}=\frac{1}{16}+\frac{3}{16}+1=1.25`$, half-mass radius $`r_{\mathrm{half}}0.28`$, rotation period $`t(r_{\mathrm{half}})1.2`$, and binding energy $`E=1.07`$; the gas has specific internal energy $`u_{\mathrm{int}}=0.014`$. The simulations, each using a total of $`N_{\mathrm{gas}}+N_{\mathrm{stars}}+N_{\mathrm{halo}}=24576+29696+32768=87040`$ particles, were run with a new N-body/SPH code featuring adaptive smoothing and time-stepping.

Figure 1 summarizes the evolution of two of these encounters. Here $`E_\mathrm{r}`$ is the energy lost in radiative shocks, while $`\rho `$ is the gas density. Encounter RET 1:1 (left) is the most dissipative of the eight studied here; radiative losses throughout this calculation amount to $`20`$% of its initial binding energy. These losses occur in large-scale shocks as the two galaxies plow into each other at $`t=1`$, separate, and fall back together. Gas densities increase with each burst of dissipation, and by the end of the simulation some $`95`$% of the gas lies in a barely-resolved disk at the center of the merger remnant. In contrast, encounter POL 3:1 (right) is the least dissipative, losing only $`5`$% of its initial binding energy. The first passage at $`t=1`$, while close enough to produce definite tidal features, barely registers in the traces of $`E_\mathrm{r}`$ and $`\rho `$. Later passages are more dramatic, eventually driving about $`60`$% into the central regions; most of the remaining gas settles into a warped disk.

The first passage of encounter RET 1:1, shown in Figure 2, exemplifies some consequences of a violent hydrodynamic interaction. This encounter’s geometry insures that most of the gas suffers strong shocks as the galaxies intersect. By the middle frame shown here, much of the gas in the in-plane disk has been swept into the center. This inflow is driven not by gravitational torques but by hydrodynamic forces; the gas loses its spin angular momentum by colliding with gas in its companion. Gas which escapes being swept inward forms a plume connecting the two galaxies; such structures may be fairly common in interpenetrating encounters (e.g. Condon et al. 1993).

After a tidal encounter, bridges and tails develop from any disk which suffers a reasonably direct passage. While the ends of such features may escape, much of this material is bound to fall back onto its parent – or, if the galaxies have already merged, onto their combined hulks. An example of such re-accretion is illustrated in Figure 3. Here the upper row shows gas in the larger disk of encounter DIR 3:1 responding after the direct passage of its lighter companion; note the pronounced bar typically formed in such passages. The lower row shows only the gas which has fallen back from the tidal features. Returning on elongated trajectories, this material is forced onto more circular orbits by shocks, some of which are visible as narrow curvilinear features. By the last frame the re-accreted gas has built up a large disk surrounding the bar.

One consequence of this disk rebuilding is that the gas can be highly sensitive to tidal perturbations on subsequent passages. Stellar disks can’t cool down once tidally heated, so they are less responsive and generate broader structures. Figure 4 illustrates the different responses of gas and stars in the second passage of encounter POL 1:1. The first frame shows two fairly relaxed gas disks, each wreathed in stellar debris; at this point the stars and gas have roughly similar distributions. In the second frame the disks have moved past each other, and a pronounced tail extends to the right of the lower galaxy. By the third frame the gas tail extends far beyond the stellar tail, crossing several contours of the stellar distribution; this frame also shows the two galaxies in the process of merging.

In these eight encounters, between $`50`$% and $`95`$% of the gas falls into a compact cloud at the center of each merger remnant. This is consistent with the results of earlier calculations (Negroponte & White 1983; Noguchi 1988; Hernquist 1989; Barnes & Hernquist 1991).

Most of the gas which doesn’t wind up in the central cloud nonetheless falls back into the remnant and settles onto closed, non-intersecting orbits. Thus at later times these remnants develop extra-nuclear disks and rings of gas with rather complex morphologies, as shown by the three examples in Figure 5. The remnant of encounter DIR 1:1 contains two rings, both rather sparse; the outer one lies within $`10^{}`$ the equatorial plane, while the inner one is tilted by $`70^{}`$. The origin of this tilt is unclear; the tilted ring is largely made up of gas from the $`i_2=71^{}`$ disk, but it’s born in a plane roughly perpendicular to the initial plane of its progenitor. In the remnant of POL 1:1 the extra-nuclear disk contains $`30`$% of the gas and exhibits an integral-sign warp; this remnant also has a compact nuclear disk, to be discussed below. Finally, the remnant of POL 3:1 has several nested disks and an extended off-center ring composed of material recently accreted from tidal debris. While the outer contour of the stellar distribution roughly parallels the gas ring, the stellar density falls monotonically with distance from the center; unlike the gas, the stars don’t form a true ring.

Some of the remnants in this survey have counterrotating nuclear disks. The most striking example was produced by encounter POL 1:1; Figure 6 shows that the outer gas disk and the stellar component both rotate in one direction, while the inner gas disk rotates in the other direction. This encounter had the same initial disk angles $`i_1`$, $`\omega _1`$, $`i_2`$, and $`\omega _2`$ as did a previous example of counterrotation (Hernquist & Barnes 1991), but here the first passage was considerably closer. Counterrotation seems insensitive to details of the encounter and simulation code, though the physical mechanism involved is not obvious from a cursory inspection of the simulations. More work is needed to determine the roles of gravitational and hydrodynamic forces in creating counterrotating disks and the range of disk angles yielding outcomes like the one in Figure 6.

## 3. Towards Proper Thermodynamic Models

The forgoing simulations presume that a good part of the ISM behaves somewhat like an isothermal gas at a temperature of $`10^4\mathrm{K}`$. This is probably not a desperately bad approximation on scales of order a kpc, but it doesn’t address some of the most interesting gas-dynamical effects in merging galaxies – for example, starbursts (Larson & Tinsley 1978), infrared emission (Joseph & Wright 1985), molecular content (Young et al. 1984), and superwinds (Heckman et al. 1996). All but the last of these involve gas which has cooled to well below $`10^4\mathrm{K}`$.

It’s hard to model the cold component with any degree of realism; individual GMCs are about $`10^4`$ times the total mass of the ISM in a galaxy like ours, and key processes are incompletely understood. One brave attempt by Gerritsen & Icke (1997) allowed gas to cool below $`10^4\mathrm{K}`$, and defined sites of star formation via a simple Jeans criterion; the radiation field of the stars, evaluated in the optically thin limit, was used to heat the gas. This yielded a self-regulating medium with two phases at temperatures of $`10^2\mathrm{K}`$ and $`10^4\mathrm{K}`$; as a bonus, the simulations also produced a Schmidt law with index $`n1.3`$. These successes may inhere to almost any approach which balances radiative cooling with heating by star formation.

In infrared-selected galaxies the ISM reprocesses $`90`$% to $`99`$% of the total luminosity (Sanders & Mirabel 1996), and this reprocessing has thermodynamic consequences (e.g. Maloney 1999) missing in optically thin models. Radiative reprocessing could be approximated in SPH by a Monte-Carlo procedure or by solving the diffusion equations, but galactic-scale models will require dramatic improvements in both spatial and temporal resolution; moreover, turbulent, magnetic, and radiation pressure terms may be needed as these probably supplant thermal pressure in a starburst’s ISM. This is a tall order for simulators!

If models of cold gas in starbursts are beyond the reach of current simulations then star formation may be included phenomenologically; for example, by setting $`\tau _{\mathrm{sf}}\rho ^n`$, where $`\tau _{\mathrm{sf}}`$ is the star formation timescale, $`\rho `$ is the local gas density, and $`n=0.5`$ for a Schmidt law (Katz 1992; Mihos & Hernquist 1994). This approach seems to reproduce nuclear starbursts, but doesn’t readily explain extended star formation in systems like NGC 4038/9 (Mirabel et al. 1998) and NGC 3690/IC 694. A more general approach could use the local dissipation rate $`\dot{u}`$ as a second parameter; for example, setting $`\tau _{\mathrm{sf}}\rho ^n\dot{u}^m`$. This might approximate Jog & Solomon’s (1992) model for triggered star formation via compression of molecular clouds.

While cold gas is problematic, it’s more straightforward to include hot gas in simulations of galactic collisions. If cooling is simply turned off, gas can be shock-heated to $`10^5\mathrm{K}`$ in the early stages of an encounter, and finally to the virial temperature of $`10^6\mathrm{K}`$ in a merger (e.g. Barnes 1998, Ch. 8). But such simulations fail to reproduce the superwinds often seen in starburst galaxies; purely mechanical heating doesn’t provide enough energy to drive bulk outflows. Injection of mass and energy by stellar winds, supernovae, or AGN is probably required if the simulations are to resemble real galaxies.

## 4. Future Directions

While stellar-dynamical modeling of galaxy collisions seems a fairly mature business, the same isn’t true of simulations including gas dynamics. Existing codes may do an acceptable job of simulating the behavior of a smooth gas at $`10^4\mathrm{K}`$, but the relationship between this hypothetical stuff and the interstellar medium of real galaxies is unclear. I close by listing some possible directions for further work:

* True 3-D modeling of thin disks. In published SPH calculations the smoothing radius exceeds the vertical scale height of the gas. An order of magnitude more gas particles are needed to vertically resolve thin disks.
* Interesting temperature structures. As noted above, attempts to include cold gas face many problems. Hot gas is more easily implemented, but source terms representing stellar processes are probably required to match the observations of this component.
* Self-regulating star formation. Until simulations including cold gas become practical it may be difficult to include “feedback” from star formation. Nonetheless, it may be worth trying simple tricks; e.g., disallowing star formation in recently-shocked gas.
* Magnetic field diagnostics. Polarization maps of continuum emission can provide information on flows in galaxies (e.g. Beck et al. 1999). Simulations including magnetic fields may be able to predict such polarization patterns.

### Acknowledgments.

Support for this work was provided by NASA grant NAG 5-8393 and by Space Telescope Science Institute grant GO-06430.03-95A.

## References

Barnes, J.E. 1998, in Galaxies: Interactions and Induced Star Formation, eds. D. Friedli, L. Martinet, & D. Pfenniger (Springer-Verlag: Berlin), 275

Barnes, J.E. & Hernquist, L.E. 1991, ApJ, 370, L65

Barnes, J.E. & Hernquist, L. 1996, ApJ, 471, 115

Beck, R. et al. 1999, Nature, 397, 324

Bender, R. 1990, in Dynamics and Interactions of Galaxies, ed. R. Wielen (Springer, Berlin), 232

Combes, F., Dupraz, C., & Gerin, M. 1990, in Dynamics and Interactions of Galaxies, ed. R. Wielen (Springer, Berlin), 205

Condon, J.J. et al. 1993, AJ, 105, 1730

Gerritsen, J.P.E. & Icke, V. 1997, A&A, 325, 972

Heckman, T.M. et al. 1996, ApJ, 457, 616

Hernquist, L. 1989, Nature, 340, 687

Hernquist, L. & Barnes J.E. 1991, Nature, 354, 210

Jog, C.J. & Solomon, P.M. 1992, ApJ, 387, 152

Joseph, R.D. & Wright, G.S. 1985, MNRAS, 214, 87

Katz, N. 1992, ApJ, 391, 502

Kormendy, J. & Sanders, D. 1992, ApJ, 390, L53

Larson, R.B. & Tinsley, B.M. 1978, ApJ, 219, 46

Maloney, P.R. 1999, A&AS, 266, 207

Mihos, J.C. & Hernquist, L. 1994, ApJ, 437, 611

Mirabel, I.F. et al. 1998, A&A, 333, L1

Negroponte, J. & White, S.D.M. 1983, MNRAS, 205, 1009

Noguchi, M. 1987, MNRAS, 228, 635

Noguchi, M. 1988, A&A, 203, 259

Sanders, D.B. & Mirabel, I.F. 1996, ARA&A, 34, 749

Schweizer, F. 1982, ApJ, 252, 455

Spitzer, L. & Baade, W. 1951, ApJ, 113, 413

Toomre, A. & Toomre, J. 1972, ApJ, 405, 142

Young, J.S. et al. 1984, ApJ, 287, L65

van der Hulst: How long do the counterrotating disks and warps in the models survive?

Barnes: In DIR 1:1 the tilt of the inner ring decreases with time, while in POL 1:1 the counterrotating nuclear disk seems quite stable. Warps should persist at least as long as continued infall from tidal tails feeds in misaligned material.

Dwarakanath: Is there an initial parameter space which does not lead to counterrotating gas disks at the center?

Barnes: Only one of the 1:1 remnants has a substantial counter rotating disk. And only one of the 3:1 remnants has a counterrotating disk of any kind. Counterrotation remains the exception rather than the rule.

Mathews: To what extent could the extreme central concentration of gas in post-merger galaxies be due to numerical viscosity?

Barnes: Barnes & Hernquist (1996) varied the artificial viscosity by a factor of three in either direction without changing the resulting central concentration of the gas. Moreover, SPH and sticky particle codes give similar results. It’s unlikely that numerical effects play a big role, but proving this beyond a shadow of a doubt is not easy!


text:
# Measurement of the charged pion electromagnetic form factor

## Abstract

Separated longitudinal and transverse structure functions for the reaction 1H$`(e,e^{}\pi ^+)n`$ were measured in the momentum transfer region $`Q^2`$ = 0.6 - 1.6 $`\text{(GeV}/c\text{)}^2`$ at a value of the invariant mass $`W=1.95`$ GeV. New values for the pion charge form factor were extracted from the longitudinal cross section by using a recently developed Regge model. The results indicate that the pion form factor in this region is larger than previously assumed and is consistent with a monopole parameterization fitted to very low $`Q^2`$ elastic data.

The pion occupies an important place in the study of the quark-gluon structure of hadrons. This is exemplified by the many calculations that treat the pion as one of their prime examples . One of the reasons is that the valence structure of the pion, being $`q\overline{q}`$, is relatively simple. Hence it is expected that the value of the four-momentum transfer squared $`Q^2`$, down to which a pQCD approach to the pion structure can be applied, is lower than e.g. for the nucleon. Furthermore, the asymptotic normalization of the pion wave function, in contrast to that of the nucleon, is known from the pion decay.

The charge form factor of the pion, $`F_\pi (Q^2)`$, is an essential element of the structure of the pion. Its behaviour at very low values of $`Q^2`$, which is determined by the charge radius of the pion, has been determined up to $`Q^2`$=0.28 $`\text{(GeV}/c\text{)}^2`$ from scattering high-energy pions from atomic electrons . For the determination of the pion form factor at higher values of $`Q^2`$ one has to use high-energy electroproduction of pions on a nucleon, i.e., employ the 1H$`(e,e^{}\pi ^+)n`$ reaction. For selected kinematical conditions this process can be described as quasi-elastic scattering of the electron from a virtual pion in the proton. In the t-pole approximation the longitudinal cross section $`\sigma _\mathrm{L}`$ is proportional to the square of the pion form factor. In this way the pion form factor has been studied for $`Q^2`$ values from 0.4 to 9.8 $`\text{(GeV}/c\text{)}^2`$ at CEA/Cornell  and for $`Q^2`$ = 0.7 $`\text{(GeV}/c\text{)}^2`$ at DESY . In the DESY experiment a longitudinal/transverse (L/T) separation was performed by taking data at two values of the electron energy. In the experiments done at CEA/Cornell this was done in a few cases only, and even then the resulting uncertainties in $`\sigma _\mathrm{L}`$ were so large that the L/T separated data were not used. Instead for the actual determination of the pion form factor $`\sigma _\mathrm{L}`$ was calculated by subtracting from the measured (differential) cross section a $`\sigma _\mathrm{T}`$ that was assumed to be proportional to the total virtual photon cross section. No uncertainty in $`\sigma _\mathrm{T}`$ was included in this subtraction. This means that existing values of $`F_\pi `$ above $`Q^2`$ = 0.7 $`\text{(GeV}/c\text{)}^2`$ are not based on L/T separated cross sections. This, together with the already relatively large statistical (and systematic) uncertainties of those data, precludes a meaningful comparison with theoretical calculations in that region.

Because of the excellent properties of the electron beam and experimental setup at CEBAF it is now possible to determine L/T separated cross sections with high accuracy and thus to study the pion form factor in the regime of $`Q^2`$ = 0.5 - 3.0 $`\text{(GeV}/c\text{)}^2`$. Using the High Momentum Spectrometer and the Short Orbit Spectrometer of Hall C and electron energies between 2.4 and 4.0 GeV, data for the reaction 1H$`(e,e^{}\pi ^+)n`$ were taken for central values of $`Q^2`$ of 0.6, 0.75, 1.0 and 1.6 $`\text{(GeV}/c\text{)}^2`$, at a central value of the invariant mass $`W`$ of 1.95 GeV.

The cross section for pion electroproduction can be written as

$$\frac{d^3\sigma }{dE^{}d\mathrm{\Omega }_e^{}d\mathrm{\Omega }_\pi }=\mathrm{\Gamma }_V\frac{d^2\sigma }{dtd\varphi },$$
(1)
where $`\mathrm{\Gamma }_V`$ is the virtual photon flux factor, $`\varphi `$ is the azimuthal angle of the outgoing pion with respect to the electron scattering plane and $`t`$ is the Mandelstam variable $`t=(p_\pi q)^2`$. The two-fold differential cross section can be written as

$`2\pi {\displaystyle \frac{d^2\sigma }{dtd\varphi }}`$ $`=`$ $`ϵ{\displaystyle \frac{d\sigma _\mathrm{L}}{dt}}+{\displaystyle \frac{d\sigma _\mathrm{T}}{dt}}+\sqrt{2ϵ(ϵ+1)}{\displaystyle \frac{d\sigma _{\mathrm{LT}}}{dt}}\mathrm{cos}\varphi `$ (3)
$`+ϵ{\displaystyle \frac{d\sigma _{\mathrm{TT}}}{dt}}\mathrm{cos}2\varphi ,`$
where $`ϵ`$ is the virtual-photon polarization parameter. The cross sections $`\sigma _X\frac{d\sigma _\mathrm{X}}{dt}`$ depend on $`W`$, $`Q^2`$ and $`t`$. The longitudinal cross section $`\sigma _\mathrm{L}`$ is dominated by the t-pole term, which contains $`F_\pi `$. The $`\varphi `$ acceptance of the experiment allowed the interference terms $`\sigma _{\mathrm{LT}}`$ and $`\sigma _{\mathrm{TT}}`$ to be determined. Since data were taken at two energies at every $`Q^2`$, $`\sigma _\mathrm{L}`$ could be separated from $`\sigma _\mathrm{T}`$ by means of a Rosenbluth separation.

The analysis of the experimental data included the following . Electron identification in the Short Orbit Spectrometer was done by using the combination of lead glass calorimeter and gas Cerenkov containing Freon-12 at atmospheric pressure. Pion identification in the High Momentum Spectrometer was largely done using time of flight between two scintillating hodoscope arrays. A small contamination by real electron-proton coincidences at the highest $`Q^2`$ setting was removed by a single beam-burst cut on $`e\pi ^+`$ coincidence time. Then $`Q^2`$, $`W`$, $`t`$, and the mass of the undetected neutron were reconstructed. Cuts on the latter excluded additional pion production. Backgrounds from the aluminum target window and random coincidences were subtracted. Yields were determined after correcting for tracking efficiency, pion absorption, local target-density reduction due to beam heating, and dead times. Cross sections were obtained from the yields using a detailed Monte Carlo (MC) simulation of the experiment, which included the magnets, apertures, detector geometries, realistic wire chamber resolutions, multiple scattering in all materials, reconstruction matrix elements, pion decay, muon tracking, and internal and external radiative processes.

Calibrations with the overdetermined 1H$`(e,e^{}p)`$ reaction were critical in several applications. The beam momentum and the spectrometer central momenta were determined absolutely to 0.1%, while the incident beam angle and spectrometer central angles were absolutely determined to better than 1 mrad. The spectrometer acceptances were checked by comparison of data to MC simulations. Finally, the overall absolute cross section normalization was checked. The calculated yields for $`e+p`$ elastics agreed to better than 2% with predictions based on a parameterization of the world data .

In the pion production reaction the experimental acceptances in $`W`$, $`Q^2`$ and $`t`$ were correlated. In order to minimize errors resulting from averaging the measured yields when calculating cross sections at average values of $`W`$, $`Q^2`$ and $`t`$, a phenomenological cross section model  was used in the simulation program. In this cross section model the terms representing the $`\sigma _X`$ of Eqn. (3) were optimized in an iterative fitting procedure to globally follow the $`t`$ and $`Q^2`$-dependence of the data. The dependence of the cross section on $`W`$ was assumed to follow the phase space factor $`(W^2M_p^2)^2`$.

The experimental cross sections can then be calculated from the measured and simulated yields via the relation

$$\left(\frac{d\sigma (\overline{W},\overline{Q}^2,t)}{dt}\right)_{\mathrm{exp}}=\frac{Y_{\mathrm{exp}}}{Y_{\mathrm{MC}}}\left(\frac{d\sigma (\overline{W},\overline{Q}^2,t)}{dt}\right)_{\mathrm{MC}}.$$
(4)
This was done for five bins in $`t`$ at the four $`Q^2`$-values. Here, $`Y`$ indicates that the yields were averaged over the $`W`$ and $`Q^2`$ acceptance, $`\overline{W}`$ and $`\overline{Q}^2`$ being the acceptance weighted average values for that $`t`$-bin. Even while the average values of $`W`$ and $`Q^2`$ differed slightly at high and low $`ϵ`$, the use of Eq. (4) with a MC cross section that globally reproduces the data allows one to take a common average $`(\overline{W},\overline{Q}^2)`$ value.

A representative example of the cross section as function of $`\varphi `$ is given in Figure 1. The dependence on $`\varphi `$ was used to determine the interference terms $`\sigma _{\mathrm{LT}}`$ and $`\sigma _{\mathrm{TT}}`$ after which the combination $`\sigma _{\mathrm{uns}}=\sigma _\mathrm{T}+ϵ\sigma _\mathrm{L}`$ was obtained at both the high and low electron energy in each $`t`$ bin for each $`Q^2`$ point. The statistical uncertainty in these cross sections ranges from 2 to 5%. Furthermore, there is a total systematic uncertainty of about 3%, the

most important contributions being: simulation of the detection volume (2%), dependence of the extracted cross sections on the MC cross section model (typically less than 2%), target density reduction (1%), pion absorption (1%), pion decay (1%), and the simulation of radiative processes (1%) . Since the same acceptances in $`W`$ and $`Q^2`$ and the same average values $`\overline{W}`$ and $`\overline{Q}^2`$ were used at both energies, $`\sigma _\mathrm{L}`$ and $`\sigma _\mathrm{T}`$ could be extracted via a Rosenbluth separation.

These cross sections are displayed in Figure 2. The error bars represent the combined statistical and systematic uncertainties. Since the uncertainties that are uncorrelated in the measurements at high and low electron energies are enlarged by the factor 1/($`\mathrm{\Delta }ϵ`$) in the Rosenbluth separation, where $`\mathrm{\Delta }ϵ`$ is the difference (typically 0.3) in the photon polarization between the two measurements, the total error bars on $`\sigma _\mathrm{L}`$ are typically about 10%.

The experimental data were compared to the results of a Regge model by Vanderhaeghen, Guidal and Laget (VGL) . In this model the pion electroproduction process is described as the exchange of Regge trajectories for $`\pi `$ and $`\rho `$ like particles. The only free parameters are the pion form factor and the $`\pi \rho \gamma `$ transition form factor. The model globally agrees with existing pion photo- and electroproduction data at values of $`W`$ above 2 GeV. The VGL model is compared to the data in Figure 2. The value of $`F_\pi `$ was adjusted at every $`Q^2`$ to reproduce the $`\sigma _\mathrm{L}`$ data at the lowest value of $`t`$. The transverse cross section $`\sigma _\mathrm{T}`$ is underestimated, which can possibly be attributed to resonance contributions at $`W=1.95`$ GeV that are not included in the Regge model. Varying the $`\pi \rho \gamma `$ transition form factor within reasonable bounds changes $`\sigma _\mathrm{T}`$ by up to 30%, but has a negligible influence on $`\sigma _\mathrm{L}`$, which is completely determined by the $`\pi `$ trajectory. This t-pole dominance was checked by studying the reactions 2H$`(e,e^{}\pi ^+)nn`$ and 2H$`(e,e^{}\pi ^{})pp`$, which gave within the uncertainties a ratio of unity for the longitudinal cross

sections. Hence the VGL model is still considered to be a good starting point for determining $`F_\pi `$.

The comparison with the $`\sigma _\mathrm{L}`$ data shows that the $`t`$ dependence in the VGL model is less steep than that of the experimental data. As suggested by the analysis  of older data, where a similar behaviour was observed, we attributed the discrepancy between the data and VGL to the presence of a negative background contribution to the longitudinal cross section, presumably again due to resonances. Since virtually nothing is known about the effect of these resonances on $`\sigma _\mathrm{L}`$, we proceeded on two paths to determine a trustworthy value of $`F_\pi `$. First we fitted the VGL prediction for $`\sigma _\mathrm{L}`$ to the data by adjusting $`F_\pi `$ at the lowest $`|t|`$ bin, as shown in Fig. 2, where it is assumed to be most reliable, owing to the dominant $`t`$ pole behaviour. However, since there is no reason to believe that the (negative) background is zero at the lowest $`t`$, the result is an underestimate for $`F_\pi `$. Secondly, $`F_\pi `$ was determined adding a $`Q^2`$ dependent negative background to $`\sigma _\mathrm{L}`$(VGL) and fitting it together with the value of $`F_\pi `$. The background term was taken to be independent of $`t`$. This was suggested by looking at the ’missing background’ in $`\sigma _\mathrm{T}`$, i.e., the difference between the data and VGL for $`\sigma _\mathrm{T}`$. That background is almost constant or slightly rising with $`|t|`$. Then, assuming that the background in $`\sigma _\mathrm{L}`$ has a similar $`t`$-dependence, a constant background leads to an overestimate of $`F_\pi `$. Our best estimate for $`F_\pi `$ is taken as the average of the two results. The model uncertainty (in relative units) is taken to be the same for the four $`Q^2`$ points, and equal to one half of the average of the (relative) differences. The results are listed in Table I and shown in the form of $`Q^2`$$`F_\pi `$ in Fig. 3. The error bars were propagated from the statistical and systematic uncertainties on the cross section data. The model uncertainty is displayed as the gray bar. The fact that the value of $`F_\pi `$ at $`Q^2`$ = 0.6 $`\text{(GeV}/c\text{)}^2`$ is close to the extrapolation of the model independent data from , and that the value of the background term is lower at higher $`W`$ (see below), gives some confidence in the procedure used to determine $`F_\pi `$.

For consistency we have re-analysed the older L/T separated data at $`Q^2`$ = 0.7 $`\text{(GeV}/c\text{)}^2`$ and $`W=2.19`$ GeV from DESY . We took the published cross sections and treated them in the same way as ours. The background term in $`\sigma _\mathrm{L}`$ was found to be smaller than in the Jefferson Lab data, presumably because of the larger value of $`W`$ of the DESY data, and hence the model uncertainty is smaller, too. The resulting best value for $`F_\pi `$, also shown in Fig. 3, is larger by 12% than the original result, which was obtained by using the Born term model by Gutbrod and Kramer . Here it should be mentioned that those authors used a phenomenological $`t`$-dependent function, whereas the Regge model by itself gives a good description of the $`t`$-dependence of the (unseparated) data from Ref. .

The data for $`F_\pi `$ in the region of $`Q^2`$ up to 1.6 $`\text{(GeV}/c\text{)}^2`$ globally follow a monopole form obeying the pion charge radius . It should be mentioned that the older Bebek data in this region suggested lower $`F_\pi `$ values. However, as mentioned, they did not use L/T separated cross sections, but took a prescription for $`\sigma _\mathrm{T}`$. Our measured data for $`\sigma _\mathrm{T}`$ indicate that the values used were too high, so that the values for $`F_\pi `$ came out systematically low.

In Fig. 3 the data are also compared to theoretical calculations. The model by Maris and Tandy  provides a good description of the data. It is based on the Bethe-Salpeter equation with dressed quark and gluon propagators, and includes parameters that were determined without the use of $`F_\pi `$ data. The data are also well described by the QCD sum rule plus hard scattering estimate of Ref. . Other models  were fitted to the older $`F_\pi `$ data and therefore underestimate the present data. Figure 3 also includes the results from perturbative QCD calculations .

In summary, new accurate separated cross sections for the 1H$`(e,e^{}\pi ^+)n`$ reaction have been determined in a kine-

matical region where the t-pole process is dominant. Values for $`F_\pi `$ were extracted from the longitudinal cross section using a recently developed Regge model. Since the model does not give a perfect description of the $`t`$-dependence of the data, our results for $`F_\pi `$ contain a sizeable model uncertainty. Improvements in the theoretical description of the 1H$`(e,e^{}\pi ^+)n`$ reaction hopefully will reduce those. The data globally follow a monopole form obeying the pion charge radius, and are well above values predicted by pQCD calculations.

The authors would like to thank Drs. Guidal, Laget and Vanderhaeghen for stimulating discussions and for making their computer program available to us. This work is supported by DOE and NSF (USA), FOM (Netherlands), NSERC (Canada), KOSEF (South Korea), and NATO.


text:
# Geometry of 4d Simplicial Quantum Gravity with a U(1) Gauge Field presented by T.Yukawa

## 1 Introduction

In this talk I am going to discuss on the geometrical property of space generated by the Monte Carlo simulation of 4d quantum gravity mostly with one U(1) gauge field. This is the third of a series of our reports on this topic at the past Lattice conferences. Let me first review the present status in order to recollect the motivation and aim of our work.

In 1998 there appeared a paper  which influenced us very much. It is claimed that by adding gauge matter fields a new phase, which is called the crinkled phase, appears. In some sense this result has been already predicted by the analytic calculation . Taking the conformal mode seriously into account, while the transverse mode is treated perturbatively, the string susceptibility is shown to have the similar form as the Liouville theory of 2d quantum gravity. In contrast to the 2d case adding matter fields brings the space towards a stable configuration. Although the numerical results gave us a lot of hope for further investigation along this line, there still remain a few pathological problems.

At Lattice 98 we added a new data in order to clear some of the problems. Firstly, we have observed the new phase even in the case with one gauge field. Furthermore, we did not have any pathological problems previously encountered at the transition point. In the Lattice 99 we have presented a result of the finite size scaling analysis showing that the transition is continuous and the fractal dimension measured by two independent methods both gave about 4.6 at the transition point.

Although there still remains some points such as the phase diagram to be discussed among two groups, we feel that in practice we should not worry too much about the phase diagrams as far as the phase transition remains continuous and possibly of second order.

In the followings we take our numerical results seriously, and study the consequence mostly on the geometry of space for the case with one U(1) field.

## 2 Machinery of the numerical simulation

Let me spend an ample time to describe our method employed in the numerical simulation.

$`Z(\kappa _2,\overline{N}_4)`$ $`=`$ $`{\displaystyle \underset{T}{}}{\displaystyle \frac{1}{c(T)}}{\displaystyle \mathrm{\Pi }_l𝑑A(l)}`$
$`\mathrm{exp}(S_GS_M\mathrm{\Delta }S),`$
where $`c(T)`$ is the symmetric factor of a triangulation $`T`$, and sums are taken over all possible triangulations with the topology fixed at $`S^4`$. The action is given on the 4d simplicial manifold by

$$S_G=\kappa _4N_4\kappa _2N_2,$$
(2)
for the Einstein-Hilbert term, where $`N_i`$ is the number of $`i`$-th simplex. The two parameters $`\kappa _2`$ and $`\kappa _4`$ correspond to the inverse of the gravitational constant and the cosmological constant respectively. The action,

$$S_M=\underset{t_{ijk}}{}o(t_{ijk})(A_{ij}+A_{jk}+A_{ki})^2,$$
(3)
corresponds to gauge matter fields, where $`A_{ij}`$ is the non-compact gauge fields defined on the link between vertices $`i`$ and $`j`$ with the sign convention, $`A_{ij}=A_{ji}`$. Sums are taken over all triangles $`t_{ijk}`$ in the triangulation $`T`$ with the weight factor $`o(t_{ijk})`$, which is given by the number of 4-simplices sharing the triangle $`t_{ijk}`$. The term

$$\mathrm{\Delta }S=\frac{\delta }{2}(N_4\overline{N_4})^2$$
(4)
is added to the action artificially in order to control $`N_4`$ around $`\overline{N_4}`$ by choosing the parameter $`\delta `$ appropriately.

Measurements are made at every 100 sweeps and totally 5 to 10K configurations for each point in the graphs.

## 3 Results and discussions

Our primary concern to the present work is to make sure that we are dealing with the Einstein gravity. To this end one may ask oneself the following questions: 1) Does the average property of space reflect the classical solution of Einstein equation? 2) Is there massless graviton acting between masses? Unfortunately, we do not have simulation data with large space configuration enough to answer the second question, and in this report we consider only the first question.

The first measurement we make is the distribution of local scalar curvature. The local scalar curvature at a 4-simplex $`s`$ is defined as

$$R_s=2\underset{ts}{}\frac{\delta _t}{o(t)}\frac{V_t}{V_s},$$
(5)
where $`\delta _t`$ is the deficit angle of a triangle $`t`$,

$$\delta _t=2\pi \mathrm{cos}^1\left(\frac{1}{4}\right)o(t),$$
(6)
and where $`V_t`$ and $`V_s`$ are the area and the volume of the elementary triangle and the 4-simplex respectively. In this definition the local curvature is an average of curvatures attached to 10 triangles $`t`$ which belong to the 4-simplex $`s`$. The lattice constant is chosen so that $`2V_t/V_s=1`$. In contrast to the 2d simulation the scalar curvatures distribute sharply about average values.

In order to see how the curvature distributes locally in space, we measure the curvature $`R_s(D)`$ averaged over a strip between two rings with geodesic distances $`D`$ and $`D+1`$ from a fixed 4-simplex $`s`$ on the $`S^4`$ sphere (Fig.1). From the figure the space looks homogeneous. Measuring $`R_s(D)`$ for several starting simplices we check the rotational symmetry of the space. It indeed shows that the space is fairly symmetric.

The Einstain equation with the cosmological term has a homogeneous and symmetric solution with $`S^4`$ topology of the form

$$ds^2=\frac{dr^2}{1\frac{r^2}{r_0^2}}+r^2d\mathrm{\Omega }_3^2,$$
(7)
where the parameter $`r_0`$ is related to the scalar curvature as $`R=\frac{3}{r_0^2}`$ and to the cosmological constant as $`\mathrm{\Lambda }=\frac{12}{r_0^2}`$. The space simulated by the dynamical triangulation shows fractal property with the Hausdorff dimension a little bigger than 4 , while the space obtained as a classical solution of the Einstain equation is smooth and its dimension is apparently 4.

In order to draw the coordinate frame on such a fractal space we need to rescale the distance. By comparing the D-step volume, i.e. the number of simplices covered within D-steps from a fixed simplex, to the volume within a distance $`r`$ from a point of the classical space (eq.(7)),

$$V(r)=4\pi ^2_0^r\frac{r^3dr}{\sqrt{1\frac{r^2}{r_0^2}}}.$$
(8)
The result is shown in Fig.2 by rescaling $`r=const\times D^{\frac{d_H}{4}}`$, adjusting the parameters as $`const=0.095`$ and $`d_H=4.6`$.

The second measurement is for the effect of matter. We then check the effect of ’mass’ by fixing the scalar curvature at a simplex $`s`$ as $`R_s=1.0`$. Numerically we reject those geometrical moves which change $`R_s`$ without any constraint in matter field configurations. The effects are shown for the scalar curvature in Fig.1. The average property of curvature distribution is symmetric and it can be written as a symmetric solution of the Einstein equation of $`S^4`$ universe:

$$ds^2=\frac{dr^2}{1r^2f(r)}+r^2d\mathrm{\Omega }_3^2.$$
(9)
The scalar curvature of the space is given by

$$R(r)=12f(r)+3rf^{}(r).$$
(10)
From this equation it is possible to obtain the function $`f(r)`$ numerically by measuring $`R(r)`$,resulting as

$$f(r)=\frac{1}{3r^4}_0^rR(r^{})r^3𝑑r^{}.$$
(11)
From these numerical results we find that the space with and without mass can be considered as the solution of Einstein equation with the cosmological term, if we rescale the coordinate appropriately. However, in order to check whether it is the theoretically expected solution of quantum gravity is left for the analytic calculation for the gauge matter energy concentration.

## 4 Conclusion

Taking the numerical simulation of simplicial quantum gravity with U(1) gauge fields seriously, we have measured the geometry of space. Although the microscopic structure of space exhibits fractal, the Hausdorff dimension is $`4.6`$, which is rather close to the classical dimension 4. Furthermore, the scalar curvature distributes sharply about its average, and the global geometry seems to have a correspondence to the solution of Einstein equation for both cases with and without mass. These results encourage us to proceed further for the correlation measurement in order to detect graviton. We also feel that it is now worth considering more seriously to make effort to compare numerical results to the analytic calculation along the conformal quantum gravity.


text:
# Monitoring of the accelerator beam distributions for internal target facilities

## 1 Introduction

Internal cluster target facilities – such as COSY-11  installed at the cooler synchrotron COSY-Jülich  – permit the study of the production of mesons in the proton-proton interaction with high luminosity ($`10^{31}`$ cm-2 s-1) in spite of very low target densities ($`10^{14}`$ atoms cm-2). Such conditions minimize changes in the ejectiles’ momentum vectors due to secondary scattering in the target, and hence facilitate the precise study of reactions with cross section values at the nanobarn level.

An exact extraction of absolute cross sections from the measured data demands a reliable estimation of the acceptance of the detection system. This in turn crucially depends on the accuracy of the determination of the position and dimensions of both the beam and the target. Using the COSY–11 detection system  as illustration, we will describe a method for estimating the dimensions of the proton beam, based on the momentum distribution for elastically scattered protons, which can be measured simultaneously with the investigated reaction.

At the COSY–11 facility (see Fig. 1) the production of short-lived uncharged mesons ($`\eta `$, $`\omega `$, $`\eta ^{}`$, or $`\varphi `$) and hyperons ($`\mathrm{\Lambda }`$ and $`\mathrm{\Sigma }^0`$) is investigated by means of the missing mass technique via the $`ppppX`$ and $`pppK^+Y`$ reactions, respectively . In this way, the four-momenta of positively charged particles are fully determined experimentally. The accuracy of the extracted missing mass value determines whether the signal from the given meson is visible over a background and depends on the precision of the momentum reconstruction of the registered protons or kaons. Consequently, it varies with the detector resolution, the momentum spread and the geometrical spread of the accelerator proton beam interacting with the internal cluster target beam. The momentum reconstruction is performed by tracing back trajectories from drift chambers through the dipole magnetic field to the target (Fig. 1,), which is an infinitely thin vertical line in the ideal case. In reality, however, the reactions take place in a region of finite dimensions where the beam and the target overlap (see Fig. 2). Assuming an infinitesimal target, our analysis shows a smearing out of the momentum vectors and hence of the resolution of the missing mass signal. However, if the target and beam dimensions are known it is possible to determine the average smearing of the missing mass originating from this effect and hence to infer for example the natural width of the measured mesons, provided that the beam momentum spread is known as well.

## 2 Description of the method

Part of the COSY–11 detection setup used for the registration of elastically scattered protons is shown in Figure 1. Trajectories of protons scattered in the forward direction are measured by means of two drift chambers (D1 and D2) and a scintillator hodoscope (S1), whereas the recoil protons are registered in coincidence with the forward ones using a silicon pad detector arrangement (Si) and a scintillation detector (S4). The two–body kinematics gives an unambiguous relation between the scattering angles $`\mathrm{\Theta }_1`$ and $`\mathrm{\Theta }_2`$ of the recoil and forward protons (see Fig. 1). Therefore, as seen in Figure 3, events of elastically scattered protons can be identified from the correlation line formed between the position in the silicon pad detector Si, and the scintillator hodoscope S1, the latter measured by the two drift chamber stacks. For those protons which are elastically scattered in the forward direction and deflected in the magnetic field of the dipole, the momentum vector at the target point can be determined. According to two–body kinematics, momentum components parallel and perpendicular to the beam axis form an ellipse. A section of this ellipse is shown as a solid line in Figure 4a superimposed on the data which were selected according to the correlation criterion from Figure 3 for elastically scattered events.

In Figures 3 and Figure 4a, the elastically scattered events stand out clearly against a low level background. But more importantly the mean of the elastically scattered data is significantly shifted from the expected line, indicating that the reconstructed momenta are on average larger than expected. This discrepancy is difficult to explain in terms of the uncertainties of either the proton beam momentum or the proton beam angle. For instance, to account for this discrepancy, the beam momentum must be changed by more than 120 MeV/c (dashed line in Figure 4a), which is 40 times larger than the conservatively estimated error of the absolute beam momentum ($`\pm `$ 3 MeV/c) . Similarly, the effect could have been corrected by changing the beam angle by 40 mrad (see dotted line in Figure 4a), which exceeds the uncertainty of the beam angle ($`\pm `$ 1 mrad ) by at least a factor of 40.

A reasonable agreement is obtained, however, by shifting the assumed reaction point relative to the nominal value by $`0.2`$ cm perpendicular to the beam axis towards the center of the COSY–ring, along the X-axis defined in Figure 1. The experimentally extracted momentum components obtained in this manner are shown in Figure 4b, together with the expectation depicted by the solid line. There is good agreement and the data are now spread symmetrically around the ellipse. This spread is essentially due to the finite extensions of the cluster target and the proton beam overlap (about $`\pm 0.2`$ cm as inferred from the shift in the target position required). Other contributions, including the spread of the beam momentum and still some multiple scattering events, appear to be negligible. We now consider in more detail the influence of the effective target dimensions and the spread of the COSY beam.

We assume the target beam to be described by a cylindrical pipe, with a diameter of 9 mm, homogeneously filled with protons . Furthermore, the COSY proton beam density distribution is described by Gaussian functions with standard deviations $`\sigma _X`$ and $`\sigma _Y`$ for the horizontal and vertical directions (see Fig. 2), respectively . Monte–Carlo calculations are performed varying the distance between the target and the beam centres ($`\mathrm{\Delta }_X`$), and the horizontal proton beam extension ($`\sigma _X`$).

In order to account for the angular distribution of the elastically scattered protons, an appropriate weight $`w`$ was assigned to each generated event according to the differential distributions of the cross sections measured by the EDDA collaboration . The generated events were analyzed in the same way as the experimental data. A comparison of the data from Figure 4a (data analyzed using the nominal interaction point) with the corresponding simulated histograms allows both $`\sigma _X`$ and $`\mathrm{\Delta }_X`$ to be determined. For finding an estimate of the parameters $`\sigma _X`$ and $`\mathrm{\Delta }_X`$ we construct the $`\chi ^2`$ statistic according to the method of least squares:

$$\chi ^2=\underset{i}{}\frac{\left(\alpha N_i^s+b_iN_i^e\right)^2}{\alpha ^2{\displaystyle \underset{i^{th}bin}{}}w^2+N_i^e+b_i},$$
(1)
where $`N_i^e`$ and $`N_i^s={\displaystyle \underset{i^{th}bin}{}}w`$ denote the content of the $`i^{th}`$ bin of the P–versus–P spectrum determined from experiment (Fig. 4a) and simulations, respectively. The background events $`b_i`$ in the ith bin amount to less than one per cent of the data  and were estimated by linear interpolations between the inner and outer part of a broad distribution surrounding the expected ellipse originating from elastically scattered protons. The free parameter $`\alpha `$ allows adjustment of the overall scale of the fitted Monte–Carlo histograms. Thus, by varying the $`\alpha `$ parameter, the $`\chi _{min}^2`$ for each pair of $`\sigma _X`$ and $`\mathrm{\Delta }_X`$ was established as a minimum of the $`\chi ^2(\alpha )`$ distribution.

The results of the Monte–Carlo calculations are shown in Figure 5.
Figure 5a shows the logarithm of $`\chi _{min}^2`$ as a function of the values of $`\sigma _X`$ and $`\mathrm{\Delta }_X`$. The $`\chi _{min}^2`$ function has a valley identifying a minimum and allowing the unique determination of the varied parameters $`\sigma _X`$ and $`\mathrm{\Delta }_X`$. The overall minimum of the $`\chi _{min}^2(\sigma _X,\mathrm{\Delta }_X)`$–distribution, at $`\sigma _X`$ = 0.2 cm and $`\mathrm{\Delta }_X`$ = - 0.2 cm is more evident in Figures 5b and 5c, which show the projections of the valley line onto the respective axis.

The same results were obtained when employing the Poisson likelihood $`\chi ^2`$ derived from the maximum likelihood method

$$\chi ^2=2\underset{i}{}[\alpha N_i^s+b_iN_i^e+N_i^eln(\frac{N_i^e}{\alpha N_i^s+b_i})],$$
(2)
The vertical beam extension of $`\sigma _Y`$ = 0.51 cm was established directly from the distribution of the vertical component of the particle trajectories at the centre of the target. This is possible, since for the momentum reconstruction only the origin of the track in the horizontal plane is used, while the vertical component remains a free parameter. As shown in Figure 6, the reconstructed distribution of the vertical component of the reaction points indeed can be described well by a Gaussian distribution with the $`\sigma `$ = 0.53 cm. The width of this distribution is primarily due to the vertical spread of a proton beam. The spread caused by multiple scattering and drift chambers resolution was determined to be about 0.13 cm (standard deviation) . Therefore, $`\sigma _Y=\sqrt{0.53^20.13^2}`$ 0.51 cm.

The effect of a possible drift chamber misalignment, i. e. an inexactness of the angle $`\beta `$ in Figure 1, was estimated to cause a shift in the momentum plane by 0.15 cm for $`\mathrm{\Delta }_X`$. This gives a rather large systematical bias in the estimation of the absolute value of $`\mathrm{\Delta }_X`$, but it does not affect the parameter $`\sigma _X`$, and still allows for the determination of relative beam shifts $`\mathrm{\Delta }_X`$ during the measurement cycle.

## 3 Monitoring of the beam geometry during the measurement cycle

It is possible to use the method described here to control the size and position of the beam relative to the target. This was demonstrated during the experimental run performed in February 1998 , when for the first time longitudinal and horizontal stochastic cooling were used at the COSY accelerator.
The horizontal stochastic cooling  is used to squeeze the proton beam in the horizontal direction until the beam reaches the equilibrium between the cooling and the heating due to the target. The size of an uncooled beam increases during the cycle, as can be seen in Figure 7a, which shows the spreading of the beam in the vertical plane during the 60 minutes cycle. This was expected since stochastic cooling in the vertical plane was not used in this case. The influence of the applied cooling in the horizontal plane is clearly visible in Figure 7b. During the first five minutes of the cycle the horizontal size of the beam of about $`210^{10}`$ protons was reduced by a factor of 2, reaching the equilibrium conditions, and remaining constant for the rest of the COSY cycle.

Figure 8 depicts the accuracy of the momentum reconstruction, reflected in the spread of the data, which is mainly due to the finite horizontal size of the beam and target overlap. The upper and lower panel correspond to the first and the last minute of the measurement cycle, respectively. The data were analyzed correcting for the relative target and beam shifts $`\mathrm{\Delta }_X`$ as determined by the method discussed above. The movement of the beam relative to the target during the cycle is quantified in Figure 7c. The shift of the beam denotes also changes of the average beam momentum, due to the nonzero dispersion at the target position. The beam reaches stability after some minutes. This can be understood as the equilibration between the energy losses when crossing $`1.610^6`$ times per second through the $`H_2`$ cluster target, and the power of the longitudinal stochastic cooling which can not only diminish the spread of the beam momentum but also shifts it as a whole.

## 4 Conclusions

We presented a method which allows the proton beam dimensions and its position relative to the target to be controlled for internal beam experimental facilities. The technique is based on the measurements of the momentum distribution of elastically scattered protons and its comparison with the distributions simulated with different beam and target conditions. We demonstrated the applicability of the method in controlling the proton beam during measurement when stochastic cooling was used at the synchrotron COSY for the first time.

Acknowledgements
One of the authors (P.M.) acknowledges hospitality and financial support from the Forschungszentrum Jülich. We appreciate the careful reading of the manuscript by J.N. Tan.
This research project was supported in part by the BMBF (06MS881I), the Bilateral Cooperation between Germany and Poland represented by the Internationales Büro DLR for the BMBF (PL-N-108-95) and by the Polish State Committee for Scientific Research, and by FFE grants (41266606 and 41266654) from the Forschungszentrum Jülich.


text:
# Superconductivity in Ropes of Single-Walled Carbon Nanotubes

## Abstract

We report measurements on ropes of Single Walled Carbon Nanotubes (SWNT) in low-resistance contact to non-superconducting (normal) metallic pads, at low voltage and at temperatures down to 70 mK. In one sample, we find a two order of magnitude resistance drop below 0.55 K, which is destroyed by a magnetic field of the order of 1T, or by a d.c. current greater than 2.5 $`\mu A`$. These features strongly suggest the existence of superconductivity in ropes of SWNT.

Metallic carbon nanotubes are known to be model systems for the study of 1D electronic transport . Electronic correlations are expected to lead to a breakdown of the Fermi liquid state. Nanotubes should then be described by Luttinger Liquids (LL) theories , with collective low energy excitations and no long range order. Proof of the validity of LL description in ropes was given by the measurement of a resistance diverging as a power law with temperature down to 10 K . However, this measurement was done on nanotubes separated from measuring leads by tunnel junctions. Because of Coulomb blockade , the low temperature and voltage regime were not explored. In contrast, we have developed a technique in which measuring pads are connected through low contact resistance to suspended nanotubes . We previously showed that when the contact pads are superconducting, a large supercurrent can flow through nanotubes . In this letter, we report experimental evidence of intrinsinc superconductivity below 0.55 K in ropes of carbon nanotubes connected to normal contacts.

The samples are ropes of SWNT suspended between normal metal contacts (Pt/Au bilayers). The SWNT are prepared by an electrical arc method with a mixture of nickel and yttrium as a catalyst . SWNT with diameters of the order of 1.4 nm are obtained. They are purified by the cross-flow filtration method . The tubes are usually assembled in ropes of a few hundred parallel tubes. Isolation of an individual rope and connection to measuring pads are performed according to the procedure we previously used , where ropes are soldered to melted contacts. The contact resistance is low and the tubes can be structurally characterized with a transmission electron microscope (TEM). For the three samples presented here, the contacts were trilayers of sputtered $`Al_2O_3/Pt/Au`$ of respective thicknesses 5, 3 and 200 nm. This procedure insures that the tubes do not contain any chemical dopants such as alkalis or halogens. The contacts showed no sign of superconductivity down to 50 mK. The samples were measured in a dilution

refrigerator, at temperatures ranging from 1 K to 0.05 K, through filtered lines . Magnetic fields up to 5 T could be applied perpendicularly to the contacts and the tubes. The resistance was measured by applying a small (1 nA to 10 nA, 30 Hz) a.c. current though the sample and measuring the a.c. voltage using lock-in detection.

We select samples with a room temperature (RT) resistance less than 10 k$`\mathrm{\Omega }`$. As is generally observed, we find that the resistance increases as the temperature is lowered between 300 K and 1K . Things change however below 1K, as shown in Fig. 1 for the three samples Pt1, Pt2, and Pt3, measured in magnetic fields ranging from 0 to 2.5 T. At zero field, the zero-bias resistance of Pt3 increases as T is reduced, whereas the resistances of Pt1 and Pt2 decrease drastically below $`T_1^{}=140`$ mK for Pt1 and $`T_2^{}=550`$ mK for Pt2. The resistance of Pt1 is reduced by 30% at 70 mK. That of Pt2 decreases by more than two orders of magnitude, and saturates below 100 mK at a value $`R_r=74`$ $`\mathrm{\Omega }`$. We define a transition temperature $`T_{C_2}`$ by the inflexion point of R(T). $`T_{C_2}`$ is 370 mK at zero field, decreases at higher magnetic fields, and extrapolates to zero at 1.35 T (Fig 4c). At fields above 1.25 T, the resistance increases with decreasing temperature, similarly to Pt3, and becomes independent of magnetic field. The resistance of Pt1 follows qualitatively the same trend, but the full transition did not occur down to 70 mK. Figures 2 and 3 show that in the temperature and field range where the zero-bias resistance drops, the differential resistance is strongly bias-dependent, with lower resistance at low bias. These data suggest that the rope Pt2 (and, to a lesser extent, Pt1) is intrinsically superconducting. Although the experimental data of Pt2 seem similar to those of SWNT connected to superconducting contacts , there are major differences. In particular the $`V(I),dV/dI(I)`$ do not show any supercurrent because of the existence of a finite residual resistance.

We now analyse the superconductivity in these systems, taking into account several features: the large normal contacts, the coupling between tubes within the rope, the 1D character of each tube, and their finite length compared to relevant mesoscopic and superconducting scales. The resistance of any superconducting wire measured through normal contacts (an NSN junction) cannot be zero because the number of channels in the wire is much smaller than in the contacts: a metallic SWNT, with 2 conducting channels, has a contact resistance of half the resistance quantum, $`R_Q/2`$ (where $`R_Q=h/(2e^2)`$=12.9 k$`\mathrm{\Omega }`$), even if it is superconducting. A rope of $`N_m`$ parallel metallic SWNT will have a minimum resistance of $`R_Q/(2N)`$. Therefore we use the residual resistance $`R_r=74`$ $`\mathrm{\Omega }`$ of Pt2 to deduce that Pt2 has at least $`N_m=R_Q/2R_r90`$ metallic tubes. This is approximately one quarter of the number of tubes in the rope, measured by TEM (Fig 1c). Similarly, $`R_Q`$ is also the maximum resistance of any phase coherent metallic wire. As a consequence, the high value (9.2 k$`\mathrm{\Omega }`$) of the resistance at 1K (which corresponds to an average resistance per metallic tube of $`9.2`$ k$`\mathrm{\Omega }N_m=830`$ k$`\mathrm{\Omega }=130`$ $`R_Q`$) cannot be understood if the nanotubes are independent, unless considering a very short (unphysical) phase coherence length $`L_\phi (1K)=L/130=`$ 8 nm. On the other hand if the electrons are free to move from tube to tube, the resistance is simply explained by the presence of disorder. The mean free path is deduced from the RT resistance $`R_2=4.1`$ k$`\mathrm{\Omega }`$ through  $`l_{e2}\frac{L}{R_2}\frac{R_Q}{N_m}18`$ nm. We conclude that Pt2 is a diffusive conductor with a few hundred conduction channels.

With such a small number of channels, we expect the superconductivity to differ from 3D superconductivity. In particular, we expect to observe a broad resistance drop starting at the 3D transition temperature  $`T^{}`$ and going eventually to $`R_r`$ at zero temperature. This is what is observed in Pt2 (see figure 1.c). We estimate the gap through the BCS relation $`\mathrm{\Delta }=1.76`$ $`k_BT^{}`$ : $`\mathrm{\Delta }85`$ $`\mu eV`$ for Pt2. We can then deduce the superconducting coherence length along the tube in the diffusive limit $`\xi _2=\sqrt{\mathrm{}v_fl_e/\mathrm{\Delta }}0.3`$ $`\mu `$m where $`v_f`$ is the longitudinal Fermi velocity $`8\times 10^5`$ m/s . Consistent with 1D superconductivity, $`\xi _2`$ is ten times larger than the diameter of the rope.

Finally, reminiscent of measurements of narrow superconducting metal wires , we find jumps in the differential resistance as the current is increased (Figures 2 and 3). For Pt2 the differential resistance at low currents remains equal to $`R_r`$ up to 50 nA, where it strongly rises but does not recover its normal state value until 2.5 $`\mu `$A (fig 3a). The jump in resistance at the first step corresponds approximately to the normal state resistance of a length $`\xi _2`$ of Pt2. Each peak corresponds to a hysteretic feature in the V-I curve (fig 3b). Above 1 T the differential resistance is peaked at zero current. This is also the case for Pt3 (data not shown). The variations of the differential resistance of Pt1 are similar to those of Pt2 close to its transition temperature. These jumps are identified as phase slips , which are the occurrence of normal regions located around defects in the sample. Such phase slips can be thermally activated (TAPS), leading to an exponential decrease of the resistance instead of a sharp transition, in qualitative agreement with our experimental observation (fig 4a). At sufficiently low temperature, TAPS are replaced by quantum phase slips (QPS), which, when tunneling through the sample, contribute an additional resistance to the zero temperature resistance. Moreover, QPS are predicted to supress the transition when the normal state resistance of the sample on the phase coherence scale is larger than $`R_Q/2`$ (as confirmed by recent experiments ). Our data on Pt2 show no evidence of such an effect, even though the normal state resistance, measured above T\*, is 40% larger than $`R_Q/2`$. The current above which the jumps disappear, 2.5 $`\mu A`$, is close to the critical current $`I_C=\mathrm{\Delta }/R_re1`$ $`\mu `$A of a superconducting wire without disorder and with the same number of conducting channels . This large value of critical current would also be the maximum supercurrent in a structure with this same wire placed between superconducting contacts (with gap $`\mathrm{\Delta }_S`$), and is much larger than the Ambegaokar-Baratoff

prediction $`R_NI_C=\mathrm{\Delta }_S/e`$. This might explain the anomalously large supercurrent measured in a previous experiment , where nanotubes were connected to superconducting contacts.

We now discuss the effect of the magnetic field. The field at which the resistance saturates to its normal value and at which the critical current vanishes, 1.25 T, coincides with the field obtained by extrapolation of $`T_C(H)`$ to zero temperature (fig. 4b). It is difficult to say what causes the disappearance of superconductivity. The value of Hc(0) should be compared to the depairing field in a confined geometry , and corresponds to a flux quantum $`\mathrm{\Phi }_0`$ through a length $`\xi `$ of an individual SWNT of diameter d, $`\mu _0H_C=\mathrm{\Phi }_0/(2\sqrt{\pi }d\xi )=1.35`$ T. But $`H_C(0)`$is also close to the field $`\mu _0H_p=\mathrm{\Delta }/\mu _B=1.43`$ T at which a paramagnetic state becomes more favorable than the superconducting state . Note that this value is of the same order as the critical field that was measured on SWNT connected between superconducting contacts, i.e. much higher than the critical field of the contacts.

We now estimate the superconducting coherence length of the two other samples, to explain the extent or absence of observed transition. Indeed, investigation of the proximity effect at high-transparency NS interfaces has shown that superconductivity resists the presence of normal contacts only if the length of the superconductor is much greater than $`\xi `$ , i.e. if the wire contains a superconducting reservoir. This condition is nearly fulfilled in Pt2 ($`\xi _2L_2/3`$). Using the high temperature resistance values of Pt1 and Pt3, and assuming a gap $`\mathrm{\Delta }`$ and number of metallic tubes equal to those of Pt2 we find $`\xi _1L_1/2`$ and $`\xi _32L_3`$. These values explain qualitatively a reduced transition temperature for Pt1 and the absence of a transition for Pt3. Moreover we can argue that the superconducting transitions we see are not due to a hidden proximity effect : if the $`Al_2O_3/Pt/Au`$ contacts were made superconducting by the laser pulse, the shortest nanotube (Pt3) would become superconducting at temperatures higher than the longer tubes (Pt1 and Pt2). The main result, i.e. no visible transition with a short rope, and a visible transition in a long rope, are confirmed by measurements on two other samples which are not presented here.

We now consider the possible mechanism of superconductivity. It has been suggested that coupling with low energy phonons can turn repulsive interactions in a Luttinger liquid into attractive ones and drive the system towards a superconducting phase . Such low energy phonons have been experimentally observed in the form of mechanical bending modes of a suspended SWNT rope . It was also shown that superconducting fluctuations can dominate at low temperature in ladders such as tubes . In this case the system must be away from half-filling, a condition probably fulfilled in our experimental situation, due to hole doping from the contacts . Finally, the superconductivity reported here recalls that of graphite intercalated with alkalis (Cs,K), which also occurs between 0.2 and 0.5 K . Much higher temperatures were observed in alkali doped fullerenes  because of the coupling to higher energy phonons. This suggests the possibility of increasing the transition temperature by chemically doping the nanotubes.

We have shown that ropes of carbon nanotubes are intrinsically superconducting. This is the first observation of superconductivity in a system with such a small number of conduction channels. The understanding of this superconductivity calls for future experimental and theoretical work and motivates in particular a search of superconducting fluctuations in a single SWNT.

A.K. thanks the Russian foundation for basic research and solid state nanostructures for financial support, and thanks CNRS for a visitor’s position. We thank M. Devoret, N. Dupuis, T. Martin, D. Maslov, C. Pasquier for stimulating discussions.