text:
# A study of randomness, correlations and collectivity in the nuclear shell model

## I Introduction

Nuclei display very regular spectral patterns. Low energy states in medium- and heavy-mass even-even nuclei allow their classification in terms of seniority, anharmonic vibrator and rotor nuclei zamfi94 , according to the ratio of the excitation energies of the states $`4_1`$ and $`2_1`$. While this regular behavior has been usually related with specific forces, the investigation of the energy spectra with random interactions joh98 ; joh99  has shown that many body states are strongly correlated even in the presence of random two-body interactions. Random interactions in bosonic Hilbert spaces, like those used in the IBM and the vibron model, exhibit a large predominance of vibrational and rotational spectra, strongly suggesting that in boson spaces collectivity is an intrinsic property of the space of nuclear states bij00 ; BF00\_1 ; BF00\_2 . Recently several studies were performed with random (tbre) and displaced random ensembles in order to simulate realistic systems hor01 ; vz02 , in particular to investigate the dominance of $`0+`$ states zhao02\_2 ; zhao02\_3

In the nuclear shell model the transition between random and collective behavior in the energy spectra of 20Ne generated by two-body forces was addressed in cor82 . Collectivity was generated with a quadrupole-quadrupole force, while a residual random interaction was included in the Hamiltonian in order to study its consequences on the system’s spectroscopic properties. Both the eigenvalue distribution and the overlap between the SU(3) and calculated wave functions exhibit the smooth path from a Hamiltonian dominated by the collective force to a random, non-collective, one.

The present work aims to extend these ideas by studying the transition from a realistic parametrization of the two-body force to a purely random one, in complex systems like 24Mg, 44Ti and 48Cr. The realistic interactions we have chosen are the universal Wildenthal Wil  interaction for the sd-shell, and the KB3 kb3  for the fp-shell. The probability for each state in the yrast band to follow a sequence where the higher energies correspond to the states with the larger angular momentum is studied by varying the mixing between realistic and random forces in the Hamiltonian. The evolution of the average energies for each member of the band, as well as its B(E2) values, is also reported.

## II Band structure

The combination of random and realistic interactions is taken as cor82

$`H=aH_C+bH_R\text{,   with }a+b=1,`$ (1)
where $`H_C`$ is a realistic Hamiltonian and $`H_R`$ is a two body random ensemble. Both $`H_C`$ and $`H_R`$ are written as DZ

$$H=\underset{rs,tu}{}V_{rstu}^{JT}Z_{rs,JT}^{}Z_{tu,JT},$$
(2)
in terms of scalar products of the normalized pair creation operators $`Z_{ij,JT}^{}`$ and its Hermitian conjugate $`Z_{kl,JT}`$, where $`r,s,\mathrm{}`$ specify sub-shells associated with individual orbits, and $`J,T`$ the coupled angular momentum and isospin. For the realistic Hamiltonian $`V_{rstu}^{JT}`$ is the Wildenthal or $`KB3`$ interaction. For the random case the $`V_{rstu}^{JT}`$ matrix elements are taken from a two body random ensemble (TBRE), i.e., to be real and normally distributed with mean zero and width $`\sigma `$ for the off-diagonals and $`\sqrt{2}\sigma `$ for the diagonals. The values of the width $`\sigma `$ are taken from the realistic interactions: 1.34 MeV and 0.60 MeV, for the sd- and fp-shells, respectively. The parameter $`b`$ is varied from 0 to 1, to cover the different mixing from the realistic interaction to a pure random force.

### II.1 The fp shell

As a first example we take 48Cr, a rotational nucleus which has been widely studied in full shell model calculations. In the present case the diagonalization is performed in the $`f_{\frac{7}{2}}p_{\frac{3}{2}}`$ shell with the code Antoine cau89 , using a KB3 interaction pz90  without single particle energies. Fig. 1 shows the average energies $`\overline{E_J}`$ of the lowest energy state for each angular momentum J= 0, 2, …, 16 $`\mathrm{}`$, the yrast band, calculated for 960 samples of the random interaction.

The vertical axis displays the average energy of each state, as a function of the mixing parameter $`b`$. At the left hand side the realistic energies are shown, with its distinctive rotor pattern. As $`b`$ increases to the right, the order between the different members of the bands is maintained, but their relative separation changes. The average energies in the right hand side, the pure random Hamiltonian, still exhibit a band structure but have lost their quadrupole collectivity, as discussed below in connection with their B(E2) values.

The evolution of the average ground state band of 44Ti, calculated in the full pf-shell, is shown in Fig. 2 as a function of the mixing parameter $`b`$. It can be seen that the relative ordering of the states with different angular momentum is maintained, but their relative separations vary significantly in the transition from the realistic to the random Hamiltonian. However, up to $`b=0.6`$ the change is mostly a scale variation, with the nearly equidistant structure of the band keeping its form. The energy spectra evolves from a vibrational equidistant form at $`b=0`$ to a mixed, yet ordered, spectrum for pure random forces. The average ground state energy increases with the mixing to a maximum value, with a slight decrease for $`b=1`$.

Similar patterns of evolution of the energy centroids for each angular momentum, as function of the mixing parameter, were found in several other nuclei like 46,48Ca and 46Ti. While for 46Ti, the centroid of the J=2 state is lower in energy than the J=0 state for a pure random Hamiltonian, in general for these nuclei and others in the sd-shell the average energy ordering is conserved, and there is a gradual change in the energy spectra.

The width $`\sigma `$ associated with the lowest average energy for each angular momentum, calculated as the square root of the variance, is shown in Fig. 3 for the gs-band of 48Cr, as a function of the angular momentum, for different values of the mixing parameter $`b`$. The energy width increases with the mixing, from 1.5 MeV for $`b=0.2`$ to 7-9 MeV for $`b=1.0`$. Given that, as $`b`$ increases, these widths become larger than the gaps between the average energies, they indicate a strong mixing between different bands, an effect closely associated with the lack of collectivity discussed bellow. As a function of the angular momentum all the widths are essentially flat, showing a moderate increase for the largest angular momenta.

The energy widths are larger than the relative energy differences between adjacent states in the g.s. band, implying that there are individual cases in which the spectrum is not ordered. However, the number of runs with a completely ordered spectrum is large even for pure random interactions, as is discussed in detail bellow.

### II.2 The sd-shell

Calculations in the sd-shell can be done for the full shell and an arbitrary number of active particles. The realistic interaction used in this case is the universal Wildenthal interaction (USD Wil ). The nucleus 24Mg offers a rich enough system, where three bands can be studied simultaneously. They correspond, for $`b=0`$, to the ground state (gs), $`\beta `$\- and $`\gamma `$\- bands. The first two start with J=0 and contain only even-J states, while the $`\gamma `$\- band starts with J = 2 and includes states with both even and odd angular momenta, up to J = 8. Following the evolution of the average energy levels as the mixing parameter $`b`$ increases, we reconstruct the three bands for each value of $`b`$. Nevertheless, for a pure random interaction one can show that the states have lost their quadrupole collectivity (see bellow), and for this reason they do not form bands in the usual sense.

In Fig. 4 the average energies for each angular momentum, for the ground, $`\beta `$\- and $`\gamma `$\- bands in 24Mg, are presented for different values of the mixing parameter $`b`$.

The most relevant features observed in Fig. 4 are that the three bands evolve in a similar way, keeping their internal order as well as the relative separation between the bands. The correlation between increasing energy and angular momentum is strictly followed in all the averaged bands. The absolute ground state energy becomes larger as the mixing parameter increases, following the same pattern found in the fp-shell.

## III Ordering probabilities

In order to further study these systems, we have analyzed the probability that each state in the band has the usual ordering i.e, the largest the angular momentum, the highest the energy for a given band. To do so we counted the number of runs where the first J=0 state is the ground state, or the first excited state in the band, the second excited state, and so on, in the 960 runs. We did the same for the states with angular momenta J = 2, 4, … In Table I the percentages for 48Cr with a mixing of $`b=0.4`$ are listed. All the states have a probability of at least 90% to occupy its physically expected place, while the dispersion is very small.

Table II lists the probability for states of each angular momentum in 48Cr to occupy the indicated place for a purely random Hamiltonian, i.e. $`b=1`$. The probabilities of being in its expected place run from 44% for J = 2 to 83% for J = 16. States with J = 0, 2 and 4 are those which more often fail to occupy their place, as they tend to exchange positions. These results represent an extension of previous studies concerning the probability of each state to be the ground state, listed in the first column Zhao01 .

Table III displays the probability for each angular momentum state to be in a given position, in the 44Ti ground state band with $`b=1.0`$. These probabilities have values between 44% for J=2 to 95% for J = 12. The fact that the states with J=2 and 4 have less than 50% probability of occupying in their places is strongly connected with the closeness of their average energies, shown in Fig. 2 for $`b=1`$.

Tables IV, V and VI display the probability that states with different angular momenta in 24Mg have to occupy a given place in each band for $`b=1.0`$, for the ground state, $`\beta `$\- and $`\gamma `$\- band, respectively. In most cases the diagonal probability, i.e. the probability that each state occupies its expected place, is larger than 50%. The exceptions are the states with J = 2 and 3 in the $`\gamma `$-band, and those with J= 0 and 2 in the $`\beta `$-band, whose probabilities lie between 38% and 48%.

## IV Collectivity

A sensitive measure of the quadrupole collectivity of the system is the B(E2,$`2_10_{gs}`$) transition strength. The distribution of B(E2) strengths in 48Cr, for four different values of the mixing parameter $`b`$, is shown in Fig. 5. Its experimental value is 230 $`e^2b^2`$. For $`b=0.2`$, shown in insert a), the distribution is concentrated around the experimental B(E2) value. For $`b=0.4`$, insert b), the increase in the random components of the Hamiltonian leads to an important fragmentation of the B(E2) intensity with four clusters, one near zero, the second one around 75 $`e^2b^2`$, the next one near 170 $`e^2b^2`$ and the last one close to the measured value. For $`b=0.6`$, insert c), most of the B(E2) values are very small, with some intensity at the collective B(E2) values. Finally, in insert d) the distribution of B(E2) values for a purely random Hamiltonian is shown. It is strongly concentrated at very small values, showing a complete lack of quadrupole collectivity, in consonance with the findings of Ref. cor82 .

Fig. 6 displays the distribution of B(E2,$`2`$ 0) values for 44Ti. Insert a) shows the results for $`b=0.2`$, with a narrow distribution around the measured value of 147 $`e^2b^2`$. In insert b) the distribution for $`b=0.4`$ is presented, which is concentrated around the same value. For $`b=0.6`$, Fig. 5 c), the distribution is still concentrated around the collective B(E2) values despite the dominance of the random component in the Hamiltonian. However, this collectivity is completely lost when pure random forces are employed, as shown in Fig. 5 d).

A useful indicator of collectivity is the energy ratio

$`R={\displaystyle \frac{E_{4_1^+}}{E_{2_1^+}}}.`$ (3)
The value R = 2 is associated with a harmonic oscillator spectrum, while R = 3.3 characterizes a rigid rotor structure.

Fig. 7 shows the distribution of energy ratios for 48Cr. The case $`b=0.2`$, shown in insert a), does exhibit the actual rotor behavior of this nucleus. This feature remains dominant for $`b=0.4`$, insert b), while for $`b=0.6`$ the distribution of energy ratios is wide and peaked at R=1. For $`b=1.0`$, shown in Fig 7 d), the distribution is very wide, with a clear dominance of the R=1 ratio, in correspondence with the near degeneracy of the states with J= 2 and 4 for $`b=1`$, as shown in Fig. 1.

44Ti has a structure closer to a harmonic oscillator than to a rotor. This feature can be seen in Fig. 8 a), which displays the distribution of energy ratios for $`b=0.2`$, narrowly concentrated around R = 2. For $`b=0.4`$ and $`0.6`$, inserts b) and c), the vibrational structure is wider but well defined. For a pure random interaction, Fig. 8 d), the distribution is peeked at R = 1, reflecting the near degeneracy of the average energies $`\overline{E}_2`$ and $`\overline{E}_4`$. The displacement of the most probable energy ratio R to 1 is accompanied by a lack of quadrupole coherence, in consistent with the previous analysis of the B(E2) transition strengths.

## V Correlation and coherence.

Having discussed the evolution of the energy centroids as a function of the mixing parameter $`b`$, and the probability for each state to occupy a certain position in the band, it is natural to study the probability that the ground state band has all its states properly ordered. This probability, given as the percentage of results from the 960 runs which are properly ordered, is shown in Fig. 9 for 48Cr, 44Ti and 24Mg. In the three cases it is apparent that the spectrum is always ordered when the realistic interaction dominates the Hamiltonian, and that the probability for the ground state band to be ordered decreases to about 35% for a fully random Hamiltonian.

It is worth emphasizing that if these probabilities were independent of each other, the probability to find a completely order band would be the product of the probabilities for each state, i.e., the product of the diagonal elements of the matrices shown in Tables II, III and IV. However, these products are all smaller than 0.04, far smaller than the probabilities of 25 - 45 % found in these calculations. These results show that the many body states obtained with two-body random forces are strongly correlated. This correlation in their energies is not, however, associated with quadrupole coherence, which is lost for pure random forces cor82 . The subtle connection between correlation and coherence requires further analysis.

The correlation $`r`$ for each value of the mixing parameter $`b`$, between the states with energies $`E_J`$ and $`E_{J+2}`$ in the sample, can be quantitatively calculated as DP

$`r={\displaystyle \frac{cov(E_J,E_{J+2})}{\sigma _{E_J}\sigma _{E_{J+2}}}}`$ (4)
where

$`cov(E_J,E_{J+2})={\displaystyle \frac{_i^n(E_{J}^{}{}_{i}{}^{}\overline{E_J})(E_{(J+2)}^{}{}_{i}{}^{}\overline{E_{J+2}})}{n}}`$ (5)
Here, $`cov(E_J,E_{J+2})`$ is the covariance of the two distributions, $`\sigma `$ is the square root of the variance $`cov(E_J,E_J)`$, and $`n`$ is the size of the space. If $`r=1`$ the distributions are fully correlated, while if $`r0`$ there is no correlation between them. In table VII the correlation $`r`$ for pairs of energies in the ground state band of 48Cr are listed .

The presence of strong correlations between the different wave functions obtained with two-body random forces, suggests that the many-body states could be well approximated by a small number of configurations which may correspond to definite shapes, as was found for bosonic models. This would imply that the very large number of shell-model many-body states would be limited or constrained by the geometry imposed by the existence of a two-body Hamiltonian, even for the case that its components are randomly selected chau .

## VI  The DTBRE case

In vz02  it was shown that a displaced two body random ensemble (DTBRE) gives rise to coherent rotor patterns. To complement the present study we analyze in this section the transition from the realistic KB3 interaction to a DTBRE in the ground state band of 48Cr. The DTBRE correspond to matrix elements with a normal distribution centered at $`c=1.0`$ and width $`\sigma =0.6`$. Fig. 10 shows the evolution of the average energies in the ground-state band of 48Cr as a function of the mixing parameter $`b`$.

At variance from the results shown in Fig. 1, 2 and 4, when the random component of the interaction increases the absolute energies continuously decrease. As expected, the rotor structure survives the transition from realistic to random interactions with nearly no changes.

In Table VIII we show the probability for states belonging to the ground state band in 48Cr, to occupy its expected place. Notice that even for $`b=1.0`$, when the pure displaced random ensemble is employed, the probability for each state to be in its place is larger than $`88\%`$, thus exhibiting a high ordering. The percentage of fully ordered bands for each mixing is equally large: 100% for $`b=0.2`$, 89% for $`b=0.6`$ and 70% for $`b=1.0`$. If the centroid of the DTBRE is displaced to more negative values, like $`c=3.0`$, 100 % of the bands are ordered vz02 . The distribution of B(E2) probabilities exhibits a clear presence of quadrupole coherence for $`c=1`$. They are concentrated in a narrow peak around the collective B(E2) values for $`c=3`$ vz02 .

## VII Summary and conclusions

The average energies of states with different angular momentum preserve their ordering inside the band when the Hamiltonian is changed smoothly from a realistic to a random one. Ground state energies increase as a function of the mixing parameter. The quadrupole collectivity is lost when the Hamiltonian is dominated by random two-body forces, and the probability that the ground state band remains ordered diminishes to 25-45% in the random limit, which is anyway far larger than the product of the probabilities for each state to be in its place, thus exhibiting the strong correlations between the different wave functions. On the other hand, when displaced two-body random ensembles are employed, the average energies decrease as the random component increases, and the rotor pattern remains unchanged.

###### Acknowledgements.

The exact diagonalizations were performed with the ANTOINE code. This work was supported in part by Conacyt, México.


format: text

text:
# Stellar Atmosphere and Accretion Disk Models for the Hot Component in Symbiotic Stars

## 1. Introduction

Symbiotic stars are interacting binaries with a relatively large separation when compared to cataclysmic variables. Their orbital periods are of the order of months up to many decades. Symbiotic binaries consist of a cool giant star and a hot ionizing radiation source (see e.g. Kenyon 2001). The giant looses mass by a wind and not by Roche-lobe overflow (Mürset & Schmid 1999) at a rate of typically $`10^6`$ M$`{}_{}{}^{}\mathrm{yr}_{}^{1}`$. Only a few percent of this amount is accreted by the companion. The hot component (in quiescent symbiotics) is in most cases a hot white dwarf. Two known systems are hard X-ray sources and contain accreting neutron stars (Chakrabarty & Roche 1997, Masetti et al. 2002). It is not clear if an accretion disk is present in all systems. The possible formation of stable disks is supported by hydrodynamic simulations (Mastrodemos & Morris 1998, Dumm et al. 2000).

The observed high luminosity of the hot component is generated by mass accretion from the giant wind. An accretion rate of $`\dot{M}`$$`10^7`$M$`{}_{}{}^{}\mathrm{yr}_{}^{1}`$ onto a white dwarf yields a luminosity of 100 L. It is believed that hydrogen shell burning of the accreted material can be sustained above the white dwarf core and in this case $`\dot{M}`$$`10^8`$M$`{}_{}{}^{}\mathrm{yr}_{}^{1}`$ is sufficient to generate that amount of luminosity.

The white dwarf is dominating the hot component in those symbiotics where steep UV continua are observed. The derived white dwarf parameters are $`R0.1`$R, $`M`$ = 0.5–1 M, and the effective temperatures range between $`T_{\mathrm{𝖾𝖿𝖿}}`$ = 30 000 up to 200 000 K (see e.g. Mürset et al. 1991). The accretion disk is probably dominant in about 10% of the systems with flat UV continua (Kenyon 2001). The disk temperature reaches 100 000 K at the inner boundary layer and decreases outwards.

It is obvious that the analysis of emission line spectra from symbiotic stars requires a realistic modeling of the white dwarf and the disk. Occasionally, the collision of winds from the white dwarf and the disk as well as jets can contribute significantly to the ionizing radiation. In the following we present our approach to calculate synthetic spectra of the hot component. In Sect. 2. we describe the construction of white dwarf and neutron star NLTE model atmospheres, assuming planar geometry and hydrostatic equilibrium. Then we present our method to compute NLTE disk models, which assume a radial Keplerian $`\alpha `$-disk structure and a detailed vertical structure (Sect. 3.). Up to now we are neglecting stellar winds and disk winds. We also disregard jets which are observed in some symbiotics and which could result from outbursts, disk instabilities, or very high accretion rates.

## 2. Stellar Atmosphere Models

Our computer code calculates plane-parallel, hydrostatic NLTE models in radiative equilibrium. It uses an Accelerated Lambda Iteration (ALI) technique and it is described in detail by Werner & Dreizler (1999) and Werner et al. (2002). The code solves self-consistently the equations of radiation transfer, hydrostatic and radiative equilibrium, particle conservation, and the NLTE rate equations for the atomic levels populations. Full NLTE line blanketing is possible. Iron group elements are included by a statistical approach (superlevels and -lines) and an opacity sampling method. Basic input parameters are $`T_{\mathrm{𝖾𝖿𝖿}}`$, $`\mathrm{log}g`$ , and the chemical composition. The code computes the atmospheric structure and an emergent spectrum. It is routinely used for the spectral analysis in the X-ray to IR spectral regions of many objects: Central stars of planetary nebulae, white dwarfs, and subdwarfs. For white dwarfs element diffusion strongly affects the (E)UV flux, however, this is unimportant for accreting white dwarfs because of the short diffusion time scales. This process could only play a role when accretion were shut off over a period of many years. We have successfully computed and applied such models (Dreizler & Wolff 1999, Schuh et al. 2002). Some symbiotic systems harbor a magnetic white dwarf. In these cases as well as in the systems with neutron star components the magnetic field affects the emergent stellar spectrum. A polarized radiation transfer calculation is possible with an extension of our code (Deetjen et al. 2002). We already applied our code to the analysis of UV and X-ray data from symbiotic stars (Jordan et al. 1994, 1996).

An extensive grid of stellar fluxes, also relevant for white dwarfs in symbiotic systems, is available from http://astro.uni-tuebingen.de/~rauch/ (Rauch 2002). The data have been calculated for photoionization modeling, in particular with the CLOUDY code (Ferland et al. 1998).

## 3. NLTE Accretion Disk Models

A new computer code, that has been developed from the stellar atmosphere code described above, calculates disk spectra under the following assumptions (for more details see Nagel et al. (2002) and in this volume). The radial disk structure is computed assuming a stationary, Keplerian, geometrically thin $`\alpha `$-disk (Shakura & Sunyaev 1973). This model is fixed by four global input parameters: Stellar mass $`M_{}`$ and radius $`R_{}`$ of the accretor, mass accretion rate $`\dot{M}`$, and the viscosity parameter $`\alpha `$. For numerical treatment the disk is divided into a number of concentric rings. For each ring with radius $`R`$ our code calculates the detailed vertical structure, assuming a plane-parallel radiating slab (in analogy to a stellar atmosphere).

In contrast to a (planar) stellar atmosphere, which is characterized by $`T_{\mathrm{𝖾𝖿𝖿}}`$ and $`\mathrm{log}g`$ , a particular disk ring with radius $`R`$ is characterized by the following two parameters, which follow from the global disk parameters introduced above. The first parameter measures the dissipated and then radiated energy. It can be expressed in terms of an effective temperature $`T_{\mathrm{𝖾𝖿𝖿}}`$:

$$T_{\mathrm{eff}}^4(R)=3GM_{}\dot{M}/8\sigma \pi R^3(1R_{}/R)^{1/2}.$$

The second parameter is the surface mass density $`M`$ of the disk ring:

$$M(R)=[1(R_{}/R)^{1/2}]\dot{M}/3\pi \overline{w}.$$

$`\sigma `$ and $`G`$ are the Stefan-Boltzmann and gravitational constants, respectively. $`\overline{w}`$ is the depth mean of viscosity $`w(z)`$, where $`z`$ is the height above the disk mid-plane. The viscosity is given by the standard $`\alpha `$-parametrization as a function of the total (i.e. gas plus radiation) pressure, but numerous other modified versions are used in the literature. We use a formulation involving the Reynolds number, as proposed by Kriz & Hubeny (1986) and Hubeny & Hubeny (1998).

The radiation transfer equations plus vertical structure equations are then solved like in the stellar atmosphere case, but accounting for two basic differences. First, the gravitational acceleration (entering the hydrostatic equation) is not constant with depth, but increases with $`z`$. This is simply the vertical component of the gravitational acceleration exerted by the central object (self-gravitation of the disk is negligible):

$$g=zGM_{}/R^3.$$

Second, the energy equation for radiative equilibrium balances the dissipated mechanical energy and the net radiative losses:

$$9/4\rho wGM/R^3=4\pi _0^{\mathrm{}}(\eta _\nu \kappa _\nu J_\nu )𝑑\nu ,$$

where $`\rho `$, $`\eta `$, $`\kappa `$, $`J`$ denote mass density, opacity, emissivity and mean intensity, respectively. In the case of a stellar atmosphere the left-hand side of this equation vanishes and we get the usual radiative equilibrium equation.

The total observed disk spectrum, which depends on the inclination angle, is finally obtained by intensity integration over all rings accounting for rotational Doppler effects.

## 4. Summary and Outlook

We have modeling tools at hand which can be used to calculate ionizing spectra of the hot component in symbiotic stars. We can compute synthetic spectra emerging from the white dwarf (or neutron star) and from the accretion disk. Our immediate aim for the near future is the inclusion of irradiation effects of the stellar spectrum onto the disk. We also want to include effects of a disk wind onto the emergent spectrum. Wind models for the compact star are already available (see e.g. Jordan et al. 1996) and must be utilized at least for analyses of those symbiotic systems where hot stellar winds are observed (Schmid 2000).

### Acknowledgments.

K. Werner thanks H.M. Schmid for helpful discussions. This work is supported by the Deutsche Forschungsgemeinschaft (DFG) and the Deutsches Zentrum für Luft- und Raumfahrt (DLR).

## References

Chakrabarty, D., & Roche, P. 1997, ApJ, 489, 254

Deetjen, J. L., Dreizler, S., Jordan, S., & Werner, K. 2002, in Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, The ASP Conference Series, in press

Dreizler, S., & Wolff, B. 1999, A&A, 348, 189

Dumm, T., Folini, D., Nussbaumer, H., Schild, H., Schmutz, W., & Walder, R. 2000, A&A, 354, 1014

Ferland, G. J., Korista, K. T., Verner, D. A., & Ferguson, J. W., Kingdon, J. B., & Verner, E. M. 1998, PASP, 110, 761

Hubeny, I., & Hubeny, V. 1998, ApJ, 505, 558

Jordan, S., Mürset, U., & Werner, K. 1994, A&A, 283, 475

Jordan, S., Schmutz, W., Wolff, B., Werner, K., & Mürset, U. 1996, A&A, 312, 897

Kenyon, S. 2001, in Encyclopedia of Astronomy and Astrophysics, IOP Publishing Ltd. and Nature Publishing Group, 3262

Kriz, S., & Hubeny, I. 1986, Bull. Astron. Inst. Czechoslovakia, 37, 129

Masetti, N., Dal Fiume, D., Cusumano, G., et al. 2002, A&A, 382, 104

Mastrodemos, N., & Morris, M. 1998, ApJ, 497, 303

Mürset, U., & Schmid, H. M. 1999, A&AS, 137, 473

Mürset, U., Nussbaumer, H., Schmid, H. M., & Vogel, M. 1991, A&A, 248, 458

Nagel, T., Dreizler, S., & Werner, K., 2002 in Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, The ASP Conference Series, in press

Rauch, T. 2002, in Ionized Gaseous Nebulae, ed. W. J. Henney, J. Franco, M. Martos, & M. Peña, RevMexAA Conf. Series, 12, 150

Schmid, H.M. 2000, in Thermal and Ionization Aspects of Flows from Hot Stars, ed. H. Lamers & A. Sapar, ASP Conference Series, 204, 303

Schuh, S. L., Dreizler, S., & Wolff, B. 2002, A&A, 382, 164

Shakura N. I., & Sunyaev R. A. 1973, A&A, 24, 337

Werner, K., & Dreizler, S. 1999, in Computational Astrophysics, ed. H. Riffert, & K. Werner, Journal of Computational and Applied Mathematics, 109, 65

Werner, K., Deetjen, J. L., Dreizler, S., Nagel, T., Rauch, T., & Schuh, S.L. 2002, in Stellar Atmosphere Modeling, ed. I. Hubeny, D. Mihalas, & K. Werner, The ASP Conference Series, in press


text:
# UW-PT 02-19, SNUTP-02-027 Matching coefficients for improved staggered bilinears Presented by S. Sharpe. Research supported in part by US-DOE contract DE-FG03-96ER40956/A006 and by BK21/KRF-2002-003-C00033.

## Abstract

We calculate one-loop matching factors for bilinear operators composed of improved staggered fermions. We compare the results for different improvement schemes used in the recent literature, including the HYP action and an action close to the Asqtad action. We find that all improvement schemes substantially reduce the size of the one-loop contributions to matching factors. The resulting corrections are comparable to, or smaller than, those found with Wilson and domain-wall fermions.

Staggered fermions provide an attractive option for calculating weak matrix elements. The main advantages are (i) that unquenched simulations are significantly faster than with other fermion discretizations, and (ii) that the axial symmetry allows the calculation of weak matrix elements relevant to kaon mixing and decays. Progress has been slowed, however, by the large $`O(a^2)`$ errors observed with unimproved staggered fermions—in particular large staggered-flavor (hereafter “taste”) symmetry breaking—and the large size of some 1-loop contributions to matching factors.

Both problems are alleviated with improved staggered fermions. The key ingredient is the use of “smeared” links, which reduce the taste-breaking couplings of high momentum gluons to quarks . This is one part of the Symanzik improvement program applied at tree level . Smeared links are expected to substantially reduce the non-perturbative taste-symmetry breaking, and this is observed in the pion spectrum . It turns out that the largest reduction is for multiple APE smearing or HYP links . The smeared links are also expected to reduce the 1-loop contributions to matching factors, since the taste-changing couplings are responsible for large tadpole-like contributions . The purpose of our present work is to test the latter expectation with a variety of smearing choices.

We carry out this test using “hypercube bilinears”, i.e. bilinear operators in which the quark and antiquark fields reside on the same $`2^4`$ hypercube. These are the most local bilinears which match onto the complete set of continuum operators, having all spins and tastes. They are also the building blocks for the four-fermion operators we use in calculations of weak matrix elements.

We compare the following choices of improved actions and operators (and use the associated numbers to refer to them below):111For details of the definitions of the various smeared links we refer to the original references.

1. Unimproved staggered action and operators made gauge invariant using the original links.
2. Unimproved action and operators except that all links are replaced by the Fat-7 smeared links introduced in Ref. .
3. Unimproved action and operators except that all links are replaced by the $`O(a^2)`$ improved links introduced in Ref.  (these involve Fat-7 smearing and the additional double-staple “Lepage” term).
4. Unimproved action and operators except that all links are replaced by HYP links . We consider two choices of the HYP smearing parameters, denoted 4A and 4B below.
5. Asqtad-like action (includes the Naik term but uses unimproved Wilson glue instead of the improved gluon action used by the MILC collaboration ) and operators containing $`O(a^2)`$ improved links.

We tadpole-improve all links (using the fourth-root of the plaquette to define $`u_0`$) except for the HYP links, for which we discuss mean-field improvement below.

The new feature introduced into perturbative calculations by smearing is that a link in, say, the $`\mu `$-th direction couples to gluons polarized in all directions $`\nu `$, not only to $`\nu =\mu `$ as for the original links. This introduces a vertex factor

$$\delta _{\nu ,\mu }D_\mu (k)+(1\delta _{\nu ,\mu })G_{\nu ,\mu }(k)$$
(1)
where

$`D_\mu (k)`$ $`=`$ $`1d_1{\displaystyle \underset{\nu \mu }{}}\overline{s}_\nu ^2+d_2{\displaystyle \underset{\genfrac{}{}{0pt}{}{\nu <\rho }{\nu ,\rho \mu }}{}}\overline{s}_\nu ^2\overline{s}_\rho ^2`$ (2)
$`d_3\overline{s}_\nu ^2\overline{s}_\rho ^2\overline{s}_\sigma ^2d_4{\displaystyle \underset{\nu \mu }{}}\overline{s}_\nu ^4,`$
with $`\overline{s}_\nu =\mathrm{sin}(k_\nu /2)`$, etc., and

$`G_{\nu ,\mu }(k)`$ $`=`$ $`\overline{s}_\mu \overline{s}_\nu \stackrel{~}{G}_{\nu ,\mu }(k)`$ (3)
$`\stackrel{~}{G}_{\nu ,\mu }(k)`$ $`=`$ $`d_1d_2{\displaystyle \frac{(\overline{s}_\rho ^2+\overline{s}_\sigma ^2)}{2}}+d_3{\displaystyle \frac{\overline{s}_\rho ^2\overline{s}_\sigma ^2}{3}}+d_4\overline{s}_\nu ^2,`$ (4)
with all indices ($`\mu ,\nu ,\rho ,\sigma `$) different. This general form encompasses all our smearing choices and agrees with previous results .

The different smearing choices are distinguished by the values of the coefficients $`d_{14}`$:
1. Unimproved: $`d_1=d_2=d_3=d_4=0`$;
2. Fat-7 links: $`d_1=d_2=d_3=1,d_4=0`$;
3 & 5. $`O(a^2)`$ improved: $`d_1=0,d_2=d_3=d_4=1`$;
4. HYP links: We consider two choices for the smearing parameters which define the HYP link.
4A. Those from Ref. , which were determined using a non-perturbative optimization procedure: $`d_1=0.89,d_2=0.96,d_3=1.08,d_4=0`$
4B. The second choice leads to the same coefficients as for Fat-7 links (case 2 above) and thus removes $`O(a^2)`$ taste-symmetry breaking. Note, however, that the two-gluon vertices are not the same as for the Fat-7 links, so that one-loop tadpole diagrams differ.

Using these new vertices, we calculate the one-loop matching to continuum bilinears regularized in the NDR scheme:

$$𝒪_i^{\mathrm{cont}}=𝒪_i^{\mathrm{lat}}+\frac{C_Fg^2}{16\pi ^2}\underset{j}{}(\delta _{ij}2d_i\mathrm{ln}(\mu a)+c_{ij})𝒪_j^{\mathrm{lat}}$$
(5)
where $`d_i=(3,0,1)`$ for spins $`(S/P,V/A,T)`$, and the labels $`i,j`$ run over spin-tastes. For example, the continuum operators are

$$𝒪_{i=(\gamma _S\xi _F)}^{\mathrm{cont}}=\overline{Q}_{\alpha ,a}^1\gamma _S^{\alpha \beta }\xi _F^{ab}Q_{\beta ,b}^2,\xi _F=\gamma _F^{}.$$
(6)
The superscripts on $`Q`$ indicate the continuum flavor—we consider only flavor non-singlets.

The details of the calculation are presented in Ref. . To check our results, we do two independent calculations using different methods. We focus here on the diagonal matching coefficients $`c_{ii}`$—four of the off-diagonal coefficients are non-vanishing but all are small, irrespective of smearing. Results for the $`c_{ii}`$ are given in Table 1. These are the matching coefficients if the NDR scale is taken to be $`\mu =1/a`$. This is expected to be within a factor of 2 of the optimum value (“$`q^{}`$”), and in any case the coefficients depend only weakly on this scale. It is sufficient here to quote only one decimal place—more accurate results can be found in Ref. .

We see from the table that all the choices of smeared link that we consider significantly reduce the size of the largest matching coefficients. Note that, in present simulations with $`1/a2`$GeV, $`\alpha _{\overline{\mathrm{MS}}}0.19`$ so $`c_{ii}=5`$ corresponds to a $`10\%`$ correction. Thus the matching coefficient for the unimproved staggered fermion scalar bilinear ($`c=29.4`$) leads to $`60\%`$ corrections, so that perturbation theory is not reliable. For all smearing choices $`|c|<10`$, which is a considerable improvement, and makes the corrections comparable to those for Wilson ($`c_S=0.1`$, $`c_P=9.7`$, $`c_V=7.8`$, $`c_A=2.9`$, $`c_T=4.3`$ in the tadpole improvement scheme of Ref. ) and domain-wall fermions ($`c_{S/P}=11.2`$, $`c_{V/A}=5.3`$, $`c_T=2.0`$, with the domain-wall mass parameter $`M=1.7`$ ). The “best” choice is case 2, i.e. Fat-7 links.

The corrections for HYP fermions can be further reduced using mean-field improvement. This is essentially a copy of the tadpole improvement scheme applied to the smeared links. We define $`u_0^{\mathrm{SM}}`$ to be the fourth-root of the plaquette built out of smeared links. Each smeared link in the action and operators is then divided by $`u_0^{\mathrm{SM}}`$. This is straightforward to implement in simulations and in perturbation theory. We expect this to reduce the matching factors by largely removing the residual fluctuations in the smeared links, but that the reduction should be by a smaller factor than for the original tadpole improvement applied to unimproved, Fat-7 or $`O(a^2)`$ links. This is borne out by our results, which are quoted in Table 2. We have also applied this “second-level” of mean-field improvement to the Fat-7 and $`O(a^2)`$ improved links, for which the effect is to slightly increase the corrections.222It turns out that, after mean-field improvement, Fat-7 and HYP links with “Fat-7” coefficients (case 4B) have identical 1-loop matching factors.

After mean-field improvement, the HYP matching factors are nearly as small as those for Fat-7 links in Table 1. Both have corrections at the 10% level. We expect the corrections for four-fermion operators to be roughly twice this size, and this is borne out by our preliminary results for four-fermion operators with HYP links . Given that HYP links lead to much greater reduction in non-perturbative taste-symmetry breaking in the pion spectrum, we are encouraged to pursue our calculations of weak matrix elements using HYP fermions.


text:
# Ionization rates in a Bose-Einstein condensate of metastable Helium

## Abstract

We have studied ionizing collisions in a BEC of He\*. Measurements of the ion production rate combined with measurements of the density and number of atoms for the same sample allow us to estimate both the 2 and 3-body contributions to this rate. A comparison with the decay of the number of condensed atoms in our magnetic trap, in the presence of an rf-shield, indicates that ionizing collisions are largely or wholly responsible for the loss. Quantum depletion makes a substantial correction to the 3-body rate constant.

The observation of Bose-Einstein condensation of metastable helium (He in the $`2^3S_1`$ state, denoted He\*) BecHe ; BecENS  constituted a pleasant surprise for experimentalists although the possibility had been predicted theoretically Shlyapnikov:94a . Success hinged, among other things, on a strong suppression of Penning ionization in the spin-polarized, magnetically trapped gas. Too high a rate of ionization would have prevented the accumulation of a density of atoms sufficiently high to achieve evaporative cooling. The ionization rate is not completely suppressed however, and when the atomic density gets high enough, a cold, magnetically trapped sample of He\* does produce a detectable flux of ions. As shown in BecHe , this signal can even be used as a signature of BEC. The observation of ions from the condensate opens the possibility of monitoring in real time the growth kinetics of a condensate NaissanceBEC . This is an exciting prospect, but in order to quantitatively interpret the ion rate, it is necessary to know the relative contributions of 2 and 3-body collisions.

In this paper we use the unique features of metastable atoms to detect, in a single realization of a BEC, the ionization rate, the density and the number of atoms. This allows us to extract 2 and 3-body rate constants without relying on fits to non-exponential decay of the number of atoms, which require good experimental reproducibility Roberts ; 3corpsDalibar ; Burt  and are notoriously difficult to interpret quantitatively Roberts . After estimating the ionization rate constants, a comparison with the observed decay of the number of atoms reveals no evidence for collisional avalanche processes. Thus, by contrast with 87Rb Rempe , He seems to be a good candidate for studying “hydrodynamic” regime as well as the effects of quantum depletion. Indeed in our analysis of the 3-body ionization process, quantum depletion makes a substantial correction DepletionGora .

Much theoretical Shlyapnikov:94a ; venturi:99  and experimental Hill:72 ; Vassen ; BecHe ; BecENS  work has already been devoted to estimating inelastic decay rates in He\*. The dominant 2-body decay mechanisms, which we will refer to collectively as Penning ionization,

$$He^{}+He^{}\{\begin{array}{c}He^++He(1S)+e^{}\hfill \\ He_2^++e^{}\hfill \end{array}$$
(1)
are known to be suppressed by at least 3 orders of magnitude in a spin-polarized sample, but the total rate constant has not yet been measured. The theoretical estimate of the rate at 1 $`\mu `$K is $`2\times 10^{14}\text{ cm}^3\text{ s}^1`$ Shlyapnikov:94a ; venturi:99 . The 3-body reaction,

$$\begin{array}{cccc}\hfill He^{}+He^{}+He^{}& & He_2^{}& +He^{}(1mK)\hfill \\ & & & He^++He(1S)+e^{}\hfill \end{array}$$
(2)
proceeds via 3-body recombination followed by autoionization of the excited molecule. The rate has been estimated theoretically Shlyapnikov:96b  to have a value of order $`10^{26}`$ cm6 s-1. Both reactions yield one positive ion which can easily be detected in our apparatus.

We define collision rate constants according to the density loss in a thermal cloud notedefcoef : $`\frac{dn}{dt}=\frac{n}{\tau }\beta n^2Ln^3`$ with $`n`$ the local density, $`\tau `$ the (vacuum limited) lifetime of the sample, and $`\beta `$ and $`L`$ the 2-body and 3-body ionizing rate constants defined for a thermal cloud. We have assumed here that there are no other loss processes. One can calculate an expected ionization rate per trapped atom ($`\mathrm{\Gamma }`$):

$$\mathrm{\Gamma }=\frac{Ionrate}{N_0}=\frac{1}{\tau ^{}}+\frac{2}{7}\kappa _2\beta n_0+\frac{8}{63}\kappa _3Ln_0^2,$$
(3)
for a pure BEC in the Thomas-Fermi regime with a number of atoms $`N_0`$, and a peak density $`n_0`$. The numerical factors come from the integration over the parabolic spatial profile and the fact that although 2 or 3 atoms are lost in each type of collision, only 1 ion is produced. The effective lifetime $`\tau ^{}\tau `$ is due to ionizing collisions with the background gas. The factors $`\kappa _i`$ take into account the fact that the 2 and 3-particle local correlation functions are smaller than those of a thermal cloud. For a dilute BEC $`\kappa _2=1/2!`$ and $`\kappa _3=1/3!`$ DepletionGora ; Burt . Because the He\* scattering length is so large, quantum depletion lead to significant corrections DepletionGora  to the $`\kappa `$’s as we discuss below.

Much of our setup has been described previously BecHe ; Therma ; Nowak:00 . Briefly, we trap up to $`2\times 10^8`$ atoms at 1 mK in a Ioffe-Pritchard trap with a lifetime ($`\tau `$) of 90 s. We use a ”cloverleaf” configuration Mewes:96  with a bias field $`B_0=150`$ mG. The axial and radial oscillation frequencies in the harmonic trapping potential are $`\nu _{}=47\pm 3`$ Hz and $`\nu _{}=1800\pm 50`$ Hz respectively ($`\overline{\omega }/2\pi =(\nu _{}\nu _{}^2)^{1/3}=534`$ Hz). A crucial feature of our set up is the detection scheme, based on a 2 stage, single anode microchannel plate detector (MCP) placed below the trapping region. Two grids above the MCP allow us either to repel positive ions and detect only the He\* atoms, or to attract and detect positive ions produced in the trapped cloud.

To detect the ion flux, the MCP is used in counting mode: the anode pulses from each ion are amplified, discriminated with a 600 ns deadtime and processed by a counter which records the time delay between successive events. Typical count rates are between $`10^2`$ and $`10^4`$ s-1. We have checked that the correlation function of the count rate is flat, indicating that there is no double counting nor any significant time correlation in the ion production. The dark count rate is of order 1 s-1. By changing the sign of the grid voltage, we have checked that while counting ions, the neutral He\* detection rate is negligible compared to the ion rate (less than 5 %) even when the radio frequency (rf) shield is on. We estimate the ion detection efficiency by assuming that only ions which hit the open channels of the MCP (60% of the total area) are detected (with a 100% quantum efficiency). We then multiply by the transmission of the two grids $`(0.84)^2`$. Based on Refs. Deconihout ; MCP , we assume this (0.42) is an upper limit on our detection efficiency.

To find the values of $`N_0`$ and $`n_0`$ corresponding to the measured ion rate, we use the MCP to observe the time-of-flight signal (TOF) of the He\* atoms released from the rapidly switched off trap. The instantaneous count rate can be as high as $`10^6`$ s-1, and the MCP saturates when used in counting mode. To avoid this problem, we lower the MCP gain, and record the TOF signal in analog mode with a time constant of 400 microseconds. Several tests were performed to verify the linearity of the detector.

In a typical run, forced evaporative cooling takes place for 40 s, down to an rf-knife frequency of 500 KHz, about 50 kHz above the minimum of the trapping potential. Near the end of the ramp, the ion rate increases sharply, signaling the appearance of a BEC (Fig. 4 in BecHe ). After reaching the final value, the rf-knife is held on at that frequency (rf-shield). The above sequence results in a quasi pure BEC for delay times up to 15 s (see Fig. 3). By quasi-pure we mean that we see no evidence of any thermal wings in signals such as shown in the inset of Fig. 1. From tests of our fitting procedure, we estimate that the smallest thermal fraction we can distinguish is about 20 %, corresponding to a temperature on the order of the chemical potential. Runs in which thermal wings were visible were discarded.

To acquire the TOF signals corresponding to a given ion rate, we turn off the rf-shield, wait 50 ms, and then turn off the magnetic trap, and switch the MCP to analog mode. To be sure that the rf has no influence on the ion rate, we use only the number of ions observed during the 50 ms delay to get the rate. We fit the TOF signals to an inverted parabola squared as expected for a pure BEC in the Thomas-Fermi regime, and for a TOF width ($`5`$ ms) narrow compared to the mean arrival time ($`100`$ ms) BecHe . Under these assumptions, the chemical potential $`\mu `$ depends only on the TOF width, the atomic mass and the acceleration of gravity Castin , and thus can be measured quite accurately. Figure 1 shows that $`\mu `$ varies as $`N_0^{2/5}`$ as expected, over almost 2 decades in atom number. Residuals from the linear fit do not show any systematic variation which is a good indication of the detection linearity. A fit on a log-log plot gives a slope of $`0.39`$.

To determine the collision rate constants $`\beta `$ and $`L`$, we need an absolute calibration of the number of atoms and the density. As discussed in Ref. BecHe , all the atoms are not detected, and the direct calibration has a 50 % uncertainty which is responsible for the large uncertainty in the scattering length $`a`$. In fact the measurement of the chemical potential gives an accurate value for the product $`n_0a=\mu m/4\pi \mathrm{}^2`$, and with the value of $`\overline{\omega }`$ gives the product $`N_0a=(1/15)(\mathrm{}/m\overline{\omega })^{1/2}(2\mu /\mathrm{}\overline{\omega })^{5/2}`$ as well. Therefore, in the hopes that the He\* scattering length will be measured more accurately in the future, we shall express $`N_0`$ and $`n_0`$ in terms of $`a`$. In this paper, unless stated otherwise, we suppose that $`a=20\text{nm}`$, and in our conclusions we shall discuss how our results depend on $`a`$.

Figure 2 shows the ion rate per atom $`\mathrm{\Gamma }`$ versus the peak density. The densest sample corresponds to $`N_0=2\times 10^5`$ atoms and $`n_0=2.5\times 10^{13}`$ cm-3. The corresponding Thomas-Fermi radii are $`r_{}5\mu `$m and $`r_{}200\mu `$m. The vertical intercept in Fig. 2 corresponds to ionizing collisions with the background gas ($`1/\tau ^{}`$). We have independently estimated this rate using trapped thermal clouds at 1 mK and 5 $`\mu `$K, and found $`1/\tau ^{}5\times 10^3`$ s-1. This value is negligible at the scale of the figure.

The curvature seen in Fig. 2 shows that 3-body ionizing collisions are significant. Before fitting the data to get $`\beta `$ and $`L`$, we must take into account several effects. First, for 3-body collisions, quantum depletion is important. For $`T=0`$, reference DepletionGora  gives a multiplicative correction CorrectionGora  to the factor $`\kappa _3`$ of $`((1+ϵ)=(1+A\frac{72}{\sqrt{\pi }}\times \sqrt{n_0a^3}))`$, where $`A0.84`$ and comes from an integration over the spatial profile using a local density approximation. At our highest density $`ϵ0.5`$. Two-body collisions are subject to an analogous correction but approximately 3 times smaller. The fits in Fig. 2 include the density dependence of $`\kappa _{2,3}`$, associated with quantum depletion. The $`n_0^{3/2}`$ dependence introduced for 2-body collisions is far too small to explain the curvature in the data. The density dependence of $`\kappa _{2,3}`$ does not significantly improve the quality of the fit, but it significantly affects the value of the fitted value of $`L`$ (reduction of 40%).

In addition, the fact that the sample probably contains a small thermal component means that collisions between the condensed and the thermal parts must be taken into account DepletionGora ; 3corpsDalibar . Assuming a 10 % thermal population ($`\frac{\mu }{k_bT}1.1`$), we find $`\kappa _3=\frac{1}{6}(1+ϵ+ϵ^{})`$, with an additional correction $`ϵ^{}0.35`$ for the densest sample QuantumDepletionT .

Taking into account all these corrections, the fitted values of the collision rate constants notedefcoef  are: $`\beta _{20}=2.9(\pm 2.0)\times 10^{14}`$ cm3 sec-1 and $`L_{20}=8.5(\pm 5.3)\times 10^{27}`$ cm6 sec-1, where the subscripts refer to the assumed value of $`a`$. These values are in good agreement with the theoretical estimates. The error bars are estimated as follows. We fix either $`\beta `$ or $`L`$ and use the other as a fit parameter. We repeat this procedure for different values of the fixed parameter and take the range over which we can get a converging and physically reasonable fit (i.e. no negative rate constants) as the uncertainty in the fixed parameter. These error bars are highly correlated since if $`\beta `$ is increased, $`L`$ must be decreased and vice-versa. The error bars do not include the uncertainty in the absolute ion detection efficiency (see below).

Until now we have assumed $`a=20`$ nm, but current experiments give a range from 8 nm to 30 nm BecHe ; BecENS . Using Eq. 3 and our parameterization of $`n_0`$ and $`N_0`$ in terms of $`a`$, one can see that, in the absence of quantum depletion, the values of $`\beta `$ and $`L`$ extracted from our analysis would be proportional to $`a^2`$ and $`a^3`$ respectively. Taking quantum depletion into account, no simple analytical dependence exists, but one can numerically evaluate $`\beta `$ and $`L`$ vs. $`a`$ and fit the results to expansions with leading terms in $`a^2`$ and $`a^3`$ respectively. The effect of quantum depletion is negligible for $`\beta `$ ($`\beta _a\beta _{20}(\frac{a}{20})^2`$). For $`L`$, we find $`L_aL_{20}(\frac{a}{20})^3[10.23\frac{a20}{20}]`$ with $`a`$ in nm.

To test the consistency of our measurements, we have studied the decay of the number of atoms in the BEC (Fig. 3). To acquire these data, we held the BEC in the trap in the presence of the rf-shield for varying times. This study involves multiple realizations of a BEC, which typically exhibit large fluctuations in the initial atom number. We have been able to reduce this noise by using the ion signal to select only data corresponding to the same ion rate 500 ms after the end of the ramp. This time corresponds to $`t=0`$ in the figure. We have also plotted the predicted decay curve (solid line) corresponding to ionization only. This curve results from a numerical integration of the atom loss due to ionization processes, calculated from the fitted values $`\beta _{20}`$ and $`L_{20}`$. The fact that the error bars on $`\beta `$ and $`L`$ are correlated leads to a small uncertainty on the solid curve that happens to be of the same order of magnitude as the typical error bars on the data. The observed decay agrees fairly well with the solid curve, and ionization apparently accounts for most of the loss. If the ion detection efficiency were actually lower than we assume, the predicted decay due to ionization would be faster than the observed decay, an unphysical situation. From this, we conclude that our estimate of the detection efficiency is accurate and does not lead to an additional uncertainty in $`\beta `$ and $`L`$.

We have also plotted the curves obtained from the same analysis but with scattering lengths of 10 and 30 nm. The curve corresponding $`a=30`$ nm lies below the data points. Based on our analysis, this means that $`a=30`$ nm is excluded. A scattering length of $`a=25`$ nm is the largest one consistent with our data. In contrast, the decay predicted for an analysis with $`a=10\text{ nm}`$ is slower than the observed decay of the number of condensed atoms. This would mean that there are additional non-ionizing losses (contributing up to half of the total loss), and/or that we have overestimated the ion detection efficiency by a factor as large as 2. In the latter case, the rate constants $`\beta `$ and $`L`$ should be multiplied by the same factor. This results for $`a=10`$nm in a supplementary systematic uncertainty on $`\beta `$ and $`L`$ of a factor as large as 2.

Even though the peak densities of our BEC are small compared to those in alkalis, the elastic collision rate is high because of the large scattering length, and one must consider the possibility of collisional avalanches. For $`a=20`$ nm our densest cloud has a mean free path of $`(8\pi a^2\overline{n})^17\mu `$m and using the definition of Rempe  the collisional opacity is $`0.8`$. With Rb atoms this would result in avalanche processes increasing by a large amount the atom loss Rempe . In our case we have to consider secondary collisions leading to both ion production and atom loss. However, no product of an ionizing collision can produce many secondary ionizing collisions, since the corresponding mean free paths are at least two orders of magnitude larger than $`r_{}`$. Hence secondary ionization is unimportant and Eq. 3 correctly describes the ionizing processes. This conclusion is supported by our observation that no correlation exists in the time distribution of detected ions.

The good agreement between the data and the curve in Fig. 3 indicates that losses due to non-ionizing collisional avalanches are not taking place either. This is in agreement with data on elastic collisions with He+, $`He_2^+`$ and He$`(1S)`$, which have small cross sections Beijerinck . Collisions with hot He\* atoms from the reaction of Eq. 2 are more likely to play a role, but due to the higher velocity, the elastic cross section for these atoms is smaller. This is in contrast to what happens with Rb atoms Rempe  where the total cross section is enhanced due to d-wave scattering resonance.

The theoretical analysis shows that quantum depletion strongly affects the measured 3-body rate constant. One way to experimentally demonstrate this effect would be to do similar measurements with thermal clouds, and compare them with the results reported here. Absolute calibration of ion and atom detection efficiency should play no role in this comparison, if one could prove that they are the same for both situations.

We thank F. Gerbier for stimulating discussions. Supported by the European Union under grants IST-1999-11055 and HPRN-CT-2000-00125, and by the DGA grant 00.34.025.


text:
# On uncertainty relations in noncommutative quantum mechanics

## Abstract

We discuss the uncertainty relations in quantum mechanics on nocommutative plane. In particular, we show that, for a given state, at most one out of three basic nontrivial uncertainty relations can be saturated. We consider also in some detail the case of angular momentum eigenstates.

In recent years noncommutative field theories (NQFT) have become very popular, mainly due to their relation to string theory in nontrivial backgrounds  and M–theory compactifications . As in the commutative case, the low energy limit of NQFT single–particle sector can be taken which results in noncommutative counterpart of one–particle quantum mechanics. There appeared many papers devoted to the study of various aspects of noncommutative quantum mechanics -. Their authors were mainly studying the spectra of hamiltonians describing particles moving on various noncommutative backgrounds.

In the present letter we analyse some very simple aspects of noncommutative quantum mechanics on the plane. Namely, we consider the uncertainty relations following from the basic commutation rules which define the theory. In particular, we are interested in the states which saturate the uncertainty relations and show that, contrary to the commutative case, not all basic nontrivial relations can be saturated simultaneously. Then we concentrate on eigenstates of angular momentum and find that while $`x`$$`x`$ relation can be saturated exactly, $`x`$$`p`$ ones — only with arbitrary good but not ideal accuracy.

The noncommutative quantum mechanics (NCQM) on the plane is defined by the following set of commutation relations:

$`[\widehat{x}_i,\widehat{x}_j]`$ $`=i\theta ϵ_{ij}I`$ (1a)
$`[\widehat{x}_i,\widehat{p}_j]`$ $`=i\mathrm{}\delta _{ij}I`$ $`i,j=1,2`$ (1b)
$`[\widehat{p}_i,\widehat{p}_j]`$ $`=0;`$ (1c)
up to renumbering $`12`$ one can assume $`\theta 0`$.

Algebra (1) looks like a deformation of Heisenberg–Weyl algebra; however, both algebras are, in fact, equivalent: the substitution

$$\begin{array}{c}\widehat{x}_i\stackrel{~}{x}_i\frac{\theta }{2\mathrm{}}ϵ_{ij}\stackrel{~}{p}_j\hfill \\ \\ \widehat{p}_i\stackrel{~}{p}_i\hfill \end{array}$$
(2)
transforms (1) into Heisenberg–Weyl algebra. The physical contents of both theories are, in general, different and depend on the interpretation of the relevant operators.

The angular momentum operator $`\widehat{L}`$ is defined as follows

$$\widehat{L}ϵ_{ij}\widehat{x}_i\widehat{p}_j+\frac{\theta }{2\mathrm{}}\widehat{p}_i\widehat{p}_i$$
(3)
and obeys

$$\begin{array}{c}[\widehat{L},\widehat{x}_i]=i\mathrm{}ϵ_{ij}\widehat{x}_j\hfill \\ \\ [\widehat{L},\widehat{p}_i]=i\mathrm{}ϵ_{ij}\widehat{p}_j\hfill \end{array}$$
(4)
Finite rotations are represented by unitary operators

$$U(\alpha )=e^{\frac{i\alpha \widehat{L}}{\mathrm{}}}$$
(5)
Eqs. (4) imply

$$U(\alpha )\left(\begin{array}{c}\widehat{x}_1\\ \widehat{x}_2\end{array}\right)U^{}(\alpha )=\left(\begin{array}{cc}\hfill \mathrm{cos}\alpha & \hfill \mathrm{sin}\alpha \\ \hfill \mathrm{sin}\alpha & \hfill \mathrm{cos}\alpha \end{array}\right)\left(\begin{array}{c}\widehat{x}_1\\ \widehat{x}_2\end{array}\right)$$
(6)
together with the similar rule for $`p`$’s.

Commutation rules (1) imply some uncertainty relations. Let us first remind the general scheme for such relations  (see also ). Given two selfadjoint operators $`\widehat{A}`$, $`\widehat{B}`$ obeying

$$[\widehat{A},\widehat{B}]=i\widehat{C}$$
(7)
one can derive the following inequalities

$$\mathrm{\Delta }A_\psi \mathrm{\Delta }B_\psi \frac{1}{2}|\widehat{C}_\psi |$$
(8)
for any normalized state $`\psi `$; here

$$\mathrm{\Delta }A_\psi \sqrt{(\psi ,(\widehat{A}A_\psi I)^2\psi )},\text{ etc.}$$
(9)
Morover, (8) is saturated iff

$$(\widehat{A}A_\psi I)\psi =i\gamma (\widehat{B}B_\psi I)\psi $$
(10)
for some $`\gamma `$. Multiplying both sides of eq. (10) by $`\widehat{A}A_\psi I`$ and using again (10) as well as the saturated version of (8) we obtain, provided $`\gamma 0`$,

$$\begin{array}{c}(\mathrm{\Delta }A)_\psi ^2=\frac{\gamma }{2}\widehat{C}_\psi \hfill \\ \\ (\mathrm{\Delta }B)_\psi ^2=\frac{1}{2\gamma }\widehat{C}_\psi \hfill \end{array}$$
(11)
In general, the operators entering uncertainty realtions are unbounded so necessary assumptions concerning their domains should be made.

Below we shall need yet another property: assume two selfadjoint operators $`\widehat{A}`$ and $`\widehat{B}`$ commute to unit operator, $`[\widehat{A},\widehat{B}]=icI`$, $`c0`$; then neither $`\widehat{A}`$ nor $`\widehat{B}`$ have normalizable eigenvectors (in their common invariant domain).

Consider now the uncertainty relations implied by the commutationrules (1):

$`\mathrm{\Delta }x_1\mathrm{\Delta }x_2{\displaystyle \frac{1}{2}}\theta `$ (12a)
$`\mathrm{\Delta }x_1\mathrm{\Delta }p_1{\displaystyle \frac{1}{2}}\mathrm{}`$ (12b)
$`\mathrm{\Delta }x_2\mathrm{\Delta }p_2{\displaystyle \frac{1}{2}}\mathrm{}`$ (12c)
We would like to know whether there exist states saturating the above inequalities. In the ”classical” ($`\theta =0`$) case there exist states saturating simultaneously both (12b) and (12c); these are famous coherent states. For $`\theta 0`$ situation is different. In fact we have the following theorem:

###### Theorem 1

For a given state $`\psi `$ at most one of the uncertainty relations (12) can be saturated.

Define $`\stackrel{˘}{x}_i\widehat{x}_i\widehat{x}_i_\psi I`$, $`\stackrel{˘}{p}_i\widehat{p}_i\widehat{p}_i_\psi I`$; then $`\stackrel{˘}{x}_i`$, $`\stackrel{˘}{p}_i`$ satisfy the same algebra (1). Assume now that $`\psi `$ saturates both (12a) and (12b). Then

$$\begin{array}{c}\stackrel{˘}{x}_1\psi =i\gamma \stackrel{˘}{x}_2\psi \hfill \\ \\ \stackrel{˘}{x}_1\psi =i\delta \stackrel{˘}{p}_1\psi ;\hfill \end{array}$$
(13)
note that both $`\gamma 0`$, $`\delta 0`$; indeed, according to the remark made above $`\psi `$ cannot be an eigenvector (in particular — null eigenvector) of $`\stackrel{˘}{x}_1`$. Eqs. (13) imply

$$(\gamma \stackrel{˘}{x}_2\delta \stackrel{˘}{p}_1)\psi =0$$
(14)
Then, due to

$`[\gamma \stackrel{˘}{x}_2\delta \stackrel{˘}{p}_1,\stackrel{˘}{p}_2]=i\mathrm{}\gamma I,`$ $`\gamma 0,`$ (15)
$`\psi `$ cannot be normalizable. This result is not very surprising: for $`\theta =0`$ the coherent states for $`x_1`$, $`p_1`$ are dispersionless neither with respect to $`x_1`$ nor $`x_2`$; however, a priori it could happen that the states saturating (12a) and (12b) diverge in $`\theta 0`$ limit.

Assume in turn, that (12b) and (12c) are saturated for some $`\psi `$, i.e.

$$\begin{array}{ccc}\stackrel{˘}{x}_i\psi =i\gamma _i\stackrel{˘}{p}_i\psi ,\hfill & \gamma _i,\hfill & i=1,2\hfill \end{array}$$
(16)
Multiply both sides by $`ϵ_{ji}\stackrel{˘}{x}_j`$ and sum over $`i`$, $`j`$

$$\begin{array}{cc}\hfill \theta \psi & =i\underset{ij}{}ϵ_{ij}\stackrel{˘}{x}_j\stackrel{˘}{x}_i\psi =\underset{i,j}{}\gamma _iϵ_{ji}\stackrel{˘}{x}_j\stackrel{˘}{p}_i\psi \hfill \\ & =\underset{i,j}{}\gamma _iϵ_{ji}\stackrel{˘}{p}_i\stackrel{˘}{x}_j\psi =i\underset{i,j}{}\gamma _j\gamma _iϵ_{ji}\stackrel{˘}{p}_i\stackrel{˘}{p}_j\psi =0\hfill \end{array}$$
(17)
i.e. $`\theta =0`$ or $`\psi =0`$. This concludes the proof of our theorem 1.

We are often interested in eigenvectors of $`\widehat{L}`$; for example, the ground state of rotationally invariant hamiltonian is an eigenstate of $`\widehat{L}`$. Assume that $`\psi `$ is an eigenvector of $`\widehat{L}`$,

$$\widehat{L}\psi =\mathrm{}l\psi $$
(18)
i.e.

$$U(\alpha )\psi =e^{il\alpha }\psi $$
(19)
We have the following

###### Theorem 2

If $`\psi `$ is an eigenvector of $`\widehat{L}`$ then neither (12b) nor (12c) are saturated.

The proof is extremely simple. It follows from (6) and (18) that $`(\psi ,\widehat{x}_i\psi )=(U(\pi )\psi ,U(\pi )\widehat{x}_i\psi )=(U(\pi )\psi ,U(\pi )\widehat{x}_iU^{}(\pi )U(\pi )\psi )=(\psi ,U(\pi )\widehat{x}_iU^{}(\pi )\psi )=(\psi ,\widehat{x}_i\psi )`$; on the other hand, using $`U(\frac{\pi }{2})`$ we find $`(\psi ,\widehat{x}_1^2\psi )=(\psi ,\widehat{x}_2^2\psi )`$; the same result holds for $`p`$’s. Therefore $`(\mathrm{\Delta }x_1)_\psi =(\mathrm{\Delta }x_2)_\psi `$, $`(\mathrm{\Delta }p_1)_\psi =(\mathrm{\Delta }p_2)_\psi `$ and both (12b) and (12c) must be saturated simultaneously. This is, however, impossible by Theorem 1.

Although (12a) and (12b) cannot be saturated for eigenstates of $`\widehat{L}`$ they can be satisfied with arbitrary accuracy. First, as we mentioned above, for eigenstates of $`\widehat{L}`$, $`\widehat{x}_i=0`$, $`\widehat{p}_i=0`$, $`\widehat{x}_1^2=\widehat{x}_2^2`$, $`\widehat{p}_1^2=\widehat{p}_2^2`$. Using (2) we have

$$\widehat{x}_1^2=\frac{1}{2}\widehat{x}_1^2+\widehat{x}_2^2=\frac{1}{2}\frac{\theta ^2}{4\mathrm{}^2}(\stackrel{~}{p}_1^2+\stackrel{~}{p}_2^2)+(\stackrel{~}{x}_1^2+\stackrel{~}{x}_2^2)\frac{\theta }{\mathrm{}}\widehat{L}$$

Take null eigenvector of $`\widehat{L}`$ and assume it saturates both $`\mathrm{\Delta }\stackrel{~}{x}_1\mathrm{\Delta }\stackrel{~}{p}_1\frac{\mathrm{}}{2}`$, $`\mathrm{\Delta }\stackrel{~}{x}_2\mathrm{\Delta }\stackrel{~}{p}_2\frac{\mathrm{}}{2}`$ (this is possible because we are now dealing with standard Heisenberg operators). Then (again $`\stackrel{~}{x}_i=0`$, $`\stackrel{~}{x}_1^2=\stackrel{~}{x}_2^2`$, $`\stackrel{~}{p}_i=0`$, $`\stackrel{~}{p}_1^2=\stackrel{~}{p}_2^2`$)

$$\begin{array}{cc}\hfill \widehat{x}_1^2& =\frac{1}{2}\frac{\theta ^2}{4\mathrm{}^2}(\stackrel{~}{p}_1^2+\stackrel{~}{p}_2^2)+(\stackrel{~}{x}_1^2+\stackrel{~}{x}_2^2)=\hfill \\ & =\frac{\theta ^2}{4\mathrm{}^2}\stackrel{~}{p}_1^2+\stackrel{~}{x}_1^2=\frac{\theta ^2}{4\mathrm{}^2}\stackrel{~}{p}_1^2+\frac{\mathrm{}^2}{4\stackrel{~}{p}_1^2}\hfill \end{array}$$

and

$$\widehat{x}_1^2\widehat{p}_1^2=\frac{\theta ^2}{4\mathrm{}^2}(\stackrel{~}{p}_1^2)^2+\frac{\mathrm{}^2}{4}$$

Now we can take normalizable eigenstate with arbitrary small dispersion $`\stackrel{~}{p}_1^2`$; however, the state with $`\stackrel{~}{p}_1^2=0`$ is not normalizable.

On the other hand we shall see that (12a) can be saturated for eigenvectors of $`\widehat{L}`$.

We shall show now that any of the inequalities (12) can be saturated. Define formally

$$T(\theta )=e^{\frac{i\theta }{2\mathrm{}^2}\widehat{p}_1\widehat{p}_2}$$
(20)
Then $`T^1(\theta )=T^{}(\theta )=T(\theta )`$ and

$$\begin{array}{cc}\hfill \stackrel{~}{x}_1& =T(\theta )\widehat{x}_1T^{}(\theta )\hfill \\ \hfill \stackrel{~}{x}_2& =T^{}(\theta )\widehat{x}_2T(\theta )\hfill \end{array}$$
(21)
Let us remind again that $`\stackrel{~}{x}_i`$ and $`\stackrel{~}{p}_i=\widehat{p}_i`$ form standard Heisenberg–Weyl algebra. Therefore, there exist coherent $`|\varphi `$ states saturating

$$\begin{array}{cc}\mathrm{\Delta }\stackrel{~}{x}_i\mathrm{\Delta }\stackrel{~}{p}_i\frac{\mathrm{}}{2},\hfill & i=1,2\hfill \end{array}$$
(22)
Due to (21) we conclude then that $`T(\theta )|\varphi `$ ($`T(\theta )|\varphi `$) saturates (12b) ((12c)). It remains to show that the action of $`T(\pm \theta )`$ on $`|\varphi `$ is well defined. This is, however, straightforward, because the representation of the algebra (1) can be explicitly constructed using (2) and well known Fock–space technique.

Finally consider the states saturating (12a). Inspired by the standard Fock space construction and the form of the commutation rule (1a) we define creation/anihilation operators

$$\begin{array}{cc}b=\frac{1}{\sqrt{2\theta }}(\widehat{x}_1+i\widehat{x}_2),\hfill & b=\frac{1}{\sqrt{2\theta }}(\widehat{x}_1i\widehat{x}_2)\hfill \end{array}$$
(23)
or, by (2),

$$\begin{array}{cc}\hfill b& =\frac{1}{\sqrt{2\theta }}\left(\stackrel{~}{x}_1\frac{\theta }{2\mathrm{}}\stackrel{~}{p}_2+i\left(\stackrel{~}{x}_2+\frac{\theta }{2\mathrm{}}\stackrel{~}{p}_1\right)\right)=\hfill \\ & =\frac{1}{\sqrt{2\mathrm{}}}\left(\frac{1}{\sqrt{2}}\left(\sqrt{\frac{2\mathrm{}}{\theta }}\stackrel{~}{x}_1\sqrt{\frac{\theta }{2\mathrm{}}}\stackrel{~}{p}_2\right)+\frac{i}{\sqrt{2}}\left(\sqrt{\frac{2\mathrm{}}{\theta }}\stackrel{~}{x}_2+\sqrt{\frac{\theta }{2\mathrm{}}}\stackrel{~}{p}_1\right)\right)\hfill \end{array}$$
(24)
Define the dilatation operator

$$D=\frac{1}{2\mathrm{}}([\stackrel{~}{x}_1,\stackrel{~}{p}_1]_++[\stackrel{~}{x}_2,\stackrel{~}{p}_2]_+)$$
(25)
It is simple to check that

$$b=e^{\frac{i}{2}\mathrm{ln}(\frac{2\mathrm{}}{\theta })D}\left(\frac{1}{\sqrt{2}}\left(\frac{\stackrel{~}{x}_1+i\stackrel{~}{x}_2}{\sqrt{2\mathrm{}}}\right)+\frac{i}{\sqrt{2}}\left(\frac{\stackrel{~}{p}_1+i\stackrel{~}{p}_2}{\sqrt{2\mathrm{}}}\right)\right)e^{\frac{i}{2}\mathrm{ln}(\frac{2\mathrm{}}{\theta })D}$$
(26)
or

$$b=e^{\frac{i}{2}\mathrm{ln}(\frac{2\mathrm{}}{\theta })D}a_{}e^{\frac{i}{2}\mathrm{ln}(\frac{2\mathrm{}}{\theta })D},$$
(27)
where

$$a_{}=\frac{1}{\sqrt{2}}(a_1+ia_2),$$
(28)
is the anihilation operator carrying definite angular momentum. The complementary anihilation operator is

$$a_+=\frac{1}{\sqrt{2}}(a_1ia_2).$$
(29)
The following prescription can be now given for constructing the states saturating (12a). We construct the general state saturating, say, $`\mathrm{\Delta }\stackrel{~}{x}_1\mathrm{\Delta }\stackrel{~}{p}_1\frac{\mathrm{}}{2}`$ and then replace $`a_1a_{}`$, $`a_2a_+`$; finally, $`\mathrm{exp}(\frac{i}{2}\mathrm{ln}(\frac{2\mathrm{}}{\theta })D)`$ is applied. Again, one can check explicitly that the whole procedure is well defined.

Acknowledgment
Piotr Kosiński kindly acknowledges numerous discussions with Prof. Paweł Maślanka.


text:
# Dusty Rings: Signposts of Recent Planet Formation

## 1 INTRODUCTION

Every planetary system forms from a thin disk of gas and dust in orbit around a young star. In the planetesimal theory, planets grow from collisions and mergers of smaller bodies, planetesimals, embedded in the disk. Protoplanets with radii of 100 km or more stir up the remaining planetesimals along their orbits. A cascade of collisions among rapidly moving planetesimals produces a ring of dust grains, which slowly disappears as protoplanets grow into planets. This entire process can lead to a solar system similar to our own (Lissauer, 1993; Mannings, Boss, & Russell, 2000).

Recent observations support this general picture. The dusty disks around many nearby stars are as large or larger than our solar system (Koerner, 2001; Weinberger et al., 2002). The dusty ring around the nearby 10 Myr old star, HR 4796A, has a thickness of $``$ 15 AU and lies $``$ 70 AU from its central star (Jayawardhana et al., 1998; Koerner et al., 1998; Schneider et al., 1999; Greaves, Mannings, & Holland, 2000; Telesco et al., 2000). The rings or partial rings around $`ϵ`$ Eridani, Vega, and other older stars have similar dimensions (Dent et al., 2000; Koerner, Sargent, & Ostroff, 2002; Wilner et al., 2002).

There have been few detailed numerical calculations of planet formation for comparison with these modern observations. Some calculations explore the early stages of planetesimal growth in a small range of disk radii (Greenberg et al., 1984; Wetherill & Stewart, 1993). These models follow the development of a single large planet, typically the Earth or Jupiter. Others use $`n`$-body simulations to investigate the last stages of planet formation, when large bodies coalesce to form a few planets (Lissauer et al., 1996; Levison, Lissauer, & Duncan, 1998; Chambers, 2001).

To simulate the formation of an entire planetary system, the calculation must span a large range of disk radii. Extending the calculations over a large fraction of the disk allows simulation of the diversity of observable phenomena in extrasolar planetary systems. Supercomputers now allow such multiannulus planetesimal calculations covering a decade or more in disk radius (Spaute et al., 1991; Weidenschilling et al., 1997; Kortenkamp, Wetherill, & Inaba, 2001). The multiannulus calculations we discuss predict the behavior of a planet-forming disk during the early and intermediate stages of planet formation. The models yield images for comparison with observations of disks around nearby stars.

## 2 THE MODEL

To calculate planet growth, we adopt the Safronov (1969) formalism, which treats planetesimals as a statistical ensemble of bodies with a distribution of masses and velocities orbiting a central star. Because direct orbit integrations of 1020 or more planetesimals are not possible, this statistical approach is essential. As long as planetesimals are numerous and gravitational interactions are small, statistical calculations reproduce the results of direct orbit calculations (Kokubo & Ida, 2000).

To evolve the mass and velocity distributions in time, we solve the coagulation and Fokker-Planck equations for bodies undergoing inelastic collisions, drag forces, and long-range gravitational forces (Kenyon & Bromley, 2002; Weidenschilling et al., 1997). We approximate collision rates as $`n\sigma vf_g`$, where $`n`$ is the number density, $`\sigma `$ is the geometric cross-section, $`v`$ is the relative velocity, and $`f_g`$ is the gravitational focusing factor. For small planetesimals with radii $``$ 1–10 km, $`f_g`$ = 1. For larger planetesimals, the gravitational cross-section exceeds the geometric cross-section and $`f_g>`$ 1. We use the ratio of the center-of-mass collision energy $`E_c`$ to the sum of the tensile strength $`S_0`$ and the gravitational binding energy of the merged pair $`E_g`$ to assign collision outcomes. Collisions with $`E_c/(E_g+S_0)`$ 1 produce mergers with little dust; $`E_c/(E_g+S_0)`$ 1 yields mergers with dust; $`E_c/(E_g+S_0)`$ 1 yields only dust. This algorithm matches data for simulated laboratory collisions and observations of the asteroid belt in our solar system (Davis et al., 1985; Davis, Ryan, & Farinella, 1994). We then compute velocity changes from gas drag, dynamical friction and viscous stirring (Kenyon & Bromley, 2002; Ohtsuki, Stewart, & Ida, 2002). Dynamical friction transfers kinetic energy from large bodies to small bodies and drives a system to energy equipartition. Viscous stirring transfers angular momentum between bodies and increases the velocities of all planetesimals.

Our numerical calculations begin with 1–1000 m planetesimals in 64 concentric annuli at distances of 30–150 AU from a 3 $`M_{}`$ star. The total mass in planetesimals is $`M_0`$. We divide the initial continuous mass distribution of planetesimals into a differential mass distribution of 30 mass batches with spacing $`m_{i+1}/m_i`$ = 2 between successive mass batches. We add mass batches as planetesimals grow in mass. Each planetesimal batch begins with a nearly circular orbit with eccentricity $`e_0`$ and inclination $`i_0=e_0/2`$.

## 3 RESULTS

Figures 1–3 summarize results for a model with $`M_0`$ = 100 $`M_{}`$ (1 $`M_{}`$ = $`6\times 10^{27}`$ g is the mass of the Earth) and $`e_0`$ = $`10^5`$. This initial mass is appropriate for the ‘minimum mass solar nebula,’ the minimum amount of solid material needed for the planets in our solar system (Hayashi, 1981; Weidenschilling, 1977). The adopted $`e_0`$ is a reasonable equilibrium value for 1–1000 m objects. This model assumes icy bodies with a mass density, $`\rho _p`$ = 1.5 g cm-3, and a tensile strength comparable to terrestrial snow, $`S_0=10^6`$ erg g-1. These parameters are similar to those adopted for calculations of planet formation in the Kuiper Belt of our solar system (Kenyon, 2002).

We separate the growth of planetesimals into three stages. When planetesimals are small, they have small geometric cross-sections and slow growth rates. This slow growth erases the initial conditions, including $`e_0`$ and the initial mass distribution. Slow growth ends when the gravitational cross-sections of the largest objects exceed their geometric cross-sections. This ‘gravitational focusing’ enhances collision rates by a factor of 100–1000. A few of the largest objects begin ‘runaway growth’; their radii grow rapidly from $``$ 10 km to $``$ 300 km. During the runaway, dynamical friction and viscous stirring increase the eccentricities of the 1 m to 1 km radius bodies from $`e10^4`$ to $`e10^3`$. Collisions between these objects then begin to produce substantial amounts of dust. Gravitational focusing factors diminish and runaway growth ends. As the largest bodies grow slowly to 1000–3000 km sizes, they continue to stir up the smaller bodies and produce more and more dust. Eventually, nearly all of the small bodies are gone and dust production slows down. Radiation pressure from the central star removes 1 $`\mu `$m dust grains; Poynting-Robertson drag pulls 1–100 $`\mu `$m dust grains into the central star. These processes remove dust on short timescales, $``$ 10 Myr or less (Artymowicz, Burrows, & Paresce, 1989; Backman, Dasgupta, & Stencel, 1995).

Figure 1 shows the evolution of the eccentricities $`e`$ of planetesimals at the inner edge of the disk (30–36 AU). All bodies begin with the same eccentricity. As objects grow from $``$ 1 km to $``$ 1000 km, dynamical friction maintains a roughly power law eccentricity distribution for the largest bodies. Viscous stirring steadily increases the eccentricities of the smallest bodies. When $`e`$ 0.01, collisions between small bodies begin to produce dust. Continued stirring yields larger $`e`$ and more dust, until the small bodies are gone. The amount of dust then steadily decreases with time.

Figure 2 shows snapshots of the time evolution of the largest object in each annulus. For most of the evolution, the growth timescale for planetesimals is roughly $`P/\mathrm{\Sigma }A^3`$ where $`P`$ is the orbital period for heliocentric distance $`A`$ and $`\mathrm{\Sigma }`$ is the surface density of planetesimals (Lissauer, 1987). Planetesimals thus grow first at the inner edge of the disk. It takes less than 1 Myr to produce 10 km bodies in the inner disk (30 AU) and $``$ 100 Myr to produce 10 km bodies in the outer disk (150 AU). Runaway growth leads to rapid production of 100–300 km bodies on timescales of 5 Myr in the inner disk and 100–200 Myr in the outer disk. A second slow growth phase produces objects ranging in size from Pluto (1000 km radius) to Jupiter’s moon Ganymede ($``$ 2500–3000 km radius). It takes 20–30 Myr to produce these objects in the inner disk and close to 2 Gyr to make large objects in the outer disk.

To visualize dust formation in a planet-forming disk, we use the geometric optics limit to derive the optical depth $`\tau `$ of particles in the grid (Kenyon et al., 1999). For each object with radius $`r_j`$ and space density $`n_j`$ in a single annulus of width $`\delta A_i`$, the optical depth is $`\tau _{ij}=2\pi n_jr_j^2\delta A_i`$. We derive $`\tau _i`$ for each annulus as the sum of $`\tau _{ij}`$ over all objects. The relative brightness of each annulus is then proportional to the solid angle of each annulus on the sky as seen from the central star, $`L_i/L_{}`$ = $`\tau _iH_i/A_i`$, where $`H_i`$ is the scale height of the dust in annulus $`i`$. This result assumes $`\tau _i`$ 0.1 at all wavelengths, which we verify for each timestep. The brightness per unit surface area of the disk follows from $`L_i`$ and the surface area of each annulus, $`\mu _iL_i/2\pi A_i\delta A_i`$. The total luminosity is the sum of $`L_i`$ for all annuli, $`L_d/L_{}=`$ $`\tau _iH_i/A_i`$, where $`L_d/L_{}`$ is the brightness of the disk relative to the brightness of the central star.

Figure 3 shows nine snapshots of the disk viewed along the rotational axis. In each image, the bright point in the center is the star. The sequence begins with a faint disk in the upper left corner at $`t`$ = 0. From $`t`$ = 0 to $`t`$ = 1 Myr, the inner disk fades slightly as planetesimals grow from 1 km to $``$ 30–50 km. By $`t`$ = 3 Myr, the largest objects have radii of 200–300 km. Collisions between the smallest objects produce modest amounts of dust; the inner disk brightens. In the middle panels, the inner disk brightens dramatically as collisions between the smallest objects produce more and more dust. Outer regions of the disk fade and then begin to brighten as large objects begin to form and stir up smaller objects. The contrast in brightness between the bright ring (shown in blue) and the rest of the disk is a factor of 30–100, making it readily observable. In the lower panels, the bright ring moves outwards. Throughout the expansion, this ring often contains narrow dark gaps that are local minima in the dust production rate. Dust disappears in the inner disk; collisions produce more dust in the outer disk. By $`t`$ = 2.5 Gyr, the bright ring reaches the outer boundary of the disk and the entire disk starts to fade and becomes unobservable (Habing et al., 2001; Spangler et al., 2001).

In the animation of Figure 3 included in the electronic version of the paper, planet formation appears as a set of three waves propagating outward through the disk. 111This animation and an animation for Figure 4 are also available at http://cfa-www.harvard.edu/ kenyon/pf/sp/movies.html. Slow growth from 1 km to $``$ 100 km produces a dark wave that lowers the brightness of the disk. Dust formed during runaway growth produces a bright wave that appears as a series of bright, narrow, concentric rings in the disk. The disk is brightest at $`t`$ 100 Myr, when $`L_d/L_{}`$ $`10^3`$. Finally, the disappearance of the dust grains yields a second dark wave that signals the last phases of planet formation in the disk. During this last phase, large bodies remaining in the disk will coalesce and form planets.

To test the robustness of these results, we calculated a series of models with a variety of initial conditions. We changed the initial eccentricity $`e_0`$, mass density $`\rho _p`$ and tensile strength $`S_0`$ of the bodies, and the total mass $`M_0`$ in the planetesimals. Most of these parameters change the timescale but do not change the character of the evolution (Wetherill & Stewart, 1989; Kenyon & Luu, 1999b). Because gravitational focusing is less effective at larger planetesimal velocities, larger initial eccentricities, $`e_0>10^5`$, delay runaway growth and slow the progress of the bright wave through the disk. The delay is roughly $`(e_0/10^5)^{1/2}`$. Planetesimals with smaller mass densities have larger cross-sections and thus grow faster. A larger initial mass in planetesimals increases the collision rate and shortens the evolution time. The evolution time is inversely proportional to the initial mass, $`tM_0^1`$, and to the planetesimal mass density, $`t\rho _p^{2/3}`$. Finally, changes in the tensile strength $`S_0`$ affect the growth time and the final mass of planetesimals. Weaker bodies with $`S_0<10^6`$ erg g-1 produce more dust in each collision and disrupt at smaller collision velocities. More dust during linear growth enhances dynamical friction and shortens the timescale to runaway growth. More dust after runaway growth speeds the exhaustion of smaller bodies and limits the growth of the largest bodies. For $`S_0`$ = $`10^2`$ to $`10^6`$ erg g-1, the net change in evolution time is $``$ 10% to 30%. The variation in the size of the largest object is $`r_{max}S_0^{0.20\mathrm{to}0.25}`$.

The evolution is more sensitive to stochastic processes. During runaway growth, random fluctuations in the collision rate can produce a single large body which grows much more rapidly than anything around it. By robbing other bodies of material, this runaway body slows the growth of large objects nearby. Stirring by the runaway object leads to more dust on shorter timescales. These runaways always produce bright dust rings.

To illustrate this process, Figure 4 shows disk snapshots for a model with $`M_0`$ = 100 $`M_{}`$, $`e_0`$ = $`10^5`$, and $`S_0`$ = $`10^4`$ erg g-1. The electronic version of the paper contains an animation of the surface brightness evolution. For $`t<`$ 1 Myr, this model follows the evolution of the model in Figure 3. At $`t`$ 1 Myr, single runaway bodies in two adjacent annuli begin to grow much more rapidly than anything else. They grow to 1000 km in $``$ 3 Myr. Stirring by these runaway bodies produces copious amounts of dust, which we see as two bright rings at 36 AU and at 40 AU. For $`t`$ = 10–100 Myr, collisions in the inner disk exhaust the supply of small bodies. At $`t`$ 100 Myr, a fluctuation produces another runaway body at 100 AU. Its rapid growth produces another bright dust ring outside the growing inner ring. As collisions deplete small bodies in the inner disk and begin to fragment small bodies in the outer disk, both of these rings and several fainter and narrower rings propagate outward through the disk.

## 4 SUMMARY

In our calculations, multiple rings tend to form when planets grow rapidly. Rapid growth occurs when the mass density of colliding bodies $`\rho _p`$ is small and when the total mass in planetesimals $`M_0`$ is large. Bodies with low tensile strength also promote rapid growth and multiple ring production. Calculations with $`S_010^4`$ erg g-1 produce multiple rings more often than do calculations with $`S_010^6`$ erg g-1. Multiple ring production is insensitive to $`e_0`$, the initial mass distribution, and other initial conditions.

Our calculations demonstrate that planet formation produces copious amounts of dust. The dust production rate ranges from $`10^{18}`$ g yr-1 to $`10^{21}`$ g yr-1. This dust absorbs and reradiates stellar energy with a relative luminosity of $`L_{\mathrm{dust}}/L_{}`$ $`10^5`$ to $`10^3`$, comparable to observed luminosities for dusty disks surrounding nearby stars (Habing et al., 2001; Spangler et al., 2001). Our dust formation timescales of 10–100 Myr are comparable to the ages of nearby stars with dusty disks (Habing et al., 2001; Song et al., 2000; Spangler et al., 2001). Dust first forms in large quantities when the largest bodies reach sizes of $``$ 1000 km. Dust disappears when disruptive collisions exhaust the supply of $``$ 1 km bodies and radiative processes remove dust from the ring. Thus, dust is concentrated in concentric rings which propagate outward through the disk as a function of time. The outer edge of each ring marks the location where 1000 km objects are just starting to form; the inner edge marks the location where collisions have exhausted the supply of $``$ 1 km bodies and dust has disappeared. Thus, dusty rings are signposts for recent formation of 1000 km or larger planets surrounding a star.

We acknowledge a generous allotment, $``$ 500 cpu days, of computer time on the Silicon Graphics Origin-2000 ‘Alhena’ at the Jet Propulsion Laboratory through funding from the NASA Offices of Mission to Planet Earth, Aeronautics, and Space Science. Advice and comments from J. Brauman and M. Geller greatly improved our presentation. R. Mackey of the JPL supercomputing group assisted with the animations of Figure 3 and Figure 4.


text:
# Velocity Distribution of Inelastic Granular Gas in Homogeneous Cooling State

## Abstract

The velocity distribution of inelastic granular gas is examined numerically on two dimensional hard disk system in nearly elastic regime using molecular dynamical simulations. The system is prepared initially in the equilibrium state with the Maxwell-Boltzmann distribution, then after a several inelastic collisions per particle, the system falls in the state that the Boltzmann equation predicts with the stationary form of velocity distribution. It turns out, however, that due to the velocity correlation the form of the distribution function does not stay time-independent, but is gradually returning to the Maxwellian immediately after the initial transient till the clustering instability sets in. It shows that, even in the homogeneous cooling state, the velocity correlation in the inelastic system invalidates the assumption of molecular chaos and the prediction by the Boltzmann equation fails.

Free cooling of granular gas under no gravity has been attracting much interest as a subject of statistical mechanics since people realized that the inelastic collisions between particles makes the system behave very different from the elastic system, that is the subject of the conventional statistical mechanics.

Besides the cooling, or loosing its kinetic energy due to the inelasticity, it has been well recognized by now the system shows series of instabilities however small the inelasticity may be as long as the system is large enough. If the system is prepared in a highly agitated state, initially it cools down uniformly as

$$T=\frac{T_0}{(1+t/t_0)^2},$$
(1)
which is called Haff state Haff83 . After a while, the vortex structure develop in the velocity field(the shearing instability)NE00 ; then the uniformity of the particle density is broken(the clustering instability)GZ93 .

After the clustering instability, the clusters of particles collide with each other, merge, and split in a complex way LH99 ; the system eventually develops high density region, where the inelastic collapse MY92 ; MY94  is likely to happen if one consider ideal hard sphere system with a constant restitution coefficient.

Regarding the velocity distribution, the Maxwell-Boltzmann distribution is an equilibrium velocity distribution for the elastic system, and the relaxation to the distribution is known to be very fast, i.e. within a several collisions per particles when it is spatially uniform. In the case of the inelastic system, obviously any velocity distribution cannot be stationary because the system loses kinetic energy at every collision, but it is plausible that the form of distribution stays stationary after a short transient if the velocity is scaled by the average speed $`v_0(t)`$:

$$f(𝒗,t)=\frac{1}{v_0(t)^d}\widehat{f}\left(\frac{𝒗}{v_0(t)}\right).$$
(2)
In fact, kinetic theories based on the Boltzmann equation predicts that, after a several collisions per particle, the velocity distribution for the inelastic system falls into a stationary form that is different from Gaussian NE98 ; BP00 ; HOB00 .

In this report, I present results of large scale two dimensional MD simulations and shows that the form of the velocity distribution does not stay stationary in the inelastic gas, but after a short initial transient the distribution gradually getting back to the Gaussian till the clustering instability sets in. This gradual change starts at very early stage where the inhomogeneity in the system is hardly visible.

The system we examine is the two dimensional system of hard disks that undergo inelastic collisions with a constant normal restitution $`r`$. The rotational motion is ignored. Then, the collision rule is given by

$`𝒗_i^{}`$ $`=`$ $`𝒗_i{\displaystyle \frac{1+r}{2}}\left(𝒏(𝒗_i𝒗_j)\right)𝒏`$
$`𝒗_j^{}`$ $`=`$ $`𝒗_j+{\displaystyle \frac{1+r}{2}}\left(𝒏(𝒗_i𝒗_j)\right)𝒏,`$
where $`𝒗_i`$ and $`𝒗_i^{}`$ are the velocity of the $`i`$’th disk before and after the collision with the $`j`$’th particle, respectively, and $`𝒏`$ is the unit vector parallel with the relative position of the colliding particles at the time of contact.

The average speed $`v_0(t)`$ for $`d`$-dimensional system, defined by

$$\frac{d}{2}v_0(t)^2=𝑑𝒗f(𝒗,t)𝒗^2,$$
(3)
decreases as the system looses energy; with this speed we scale the velocity distribution through eq.(2). In order to see the time-dependence of the scaled velocity distribution $`\widehat{f}`$, it is convenient to expand it using the Sonine polynomial as

$$\widehat{f}(𝒄,t)=\frac{1}{\sqrt{\pi }^d}e^{c^2}\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}a_{\mathrm{}}(t)S_{\mathrm{}}(c^2)$$
(4)
when the distribution is not very different from Gaussian; the $`\mathrm{}`$’th order Sonine polynomial is the $`\mathrm{}`$’th order polynomial orthogonalized with the $`d`$-dimensional Gaussian weight function:

$$S_0(x)=1,S_1(x)=x+\frac{1}{2}d,$$

$$S_2(x)=\frac{1}{2}x^2\frac{1}{2}(d+2)x+\frac{1}{8}d(d+2),\text{ etc.}$$

Due to the normalization and scaling of $`\widehat{f}`$, we have $`a_0=1`$ and $`a_1=0`$, thus any deviation from Gaussian distribution is seen in the non-zero values of $`a_{\mathrm{}}`$ for $`\mathrm{}2`$.

Simulations were performed by the event-driven method using the fast algorithm developed by Isobe I99 . Most of the simulations were done with particle number $`N=`$250,000, number density $`n=0.25`$ (area fraction $`\varphi \pi n/4=0.196`$). I employed the periodic boundary condition and the initial state is the equilibrium state that is prepared by running the system for long enough with the restitution constant $`r=1`$. We focus on the nearly elastic regime where the system stays uniform for a substantial length of time and the distribution does not deviate very much from the Gaussian when the system is in HCS.

In the following, the time is measured by the collision time $`\tau `$, which is defined as the number of collisions each particle experiences, i.e. $`\tau 2N_{\mathrm{coll}}/N`$ with $`N_{\mathrm{coll}}`$ being the total number of collisions (the factor 2 comes from the binary collision).

Figure 1 shows the energy decay as a function of time $`\tau `$ for $`r=0.9`$ and $`n=0.25`$ in the semi-logarithmic scale. In terms of $`\tau `$, the decay in the Haff state given by (1) is expressed as

$$E(\tau )=E(0)\mathrm{exp}[\gamma \tau ],$$
(5)
where $`\gamma `$ is a decay rate. The thin solid line in Fig.1 shows the exponential function with $`\gamma =0.093`$(the line is shifted vertically to avoid complete overlapping). The initial $`\tau `$-dependence fits to the exponential decay very well with the decay rate very close to the one obtained in the case of random collision: $`\gamma _0(1r^2)/d=0.095`$ for $`r=0.9`$. It eventually deviates from the exponential around $`\tau 70`$, when the clustering instability sets in.

The speed distribution for this system is plotted in Fig.2 for $`\tau =40`$ and 80. For both cases, the distribution is very close to Gaussian and the deviation from it is hardly seen.

The deviation, however, is clearly seen in $`a_2(\tau )`$ plotted in Fig.3, where the initial deviation from the Gaussian is shown for various values of restitution constant $`r`$. From this figure, it might seem that the scaled distribution becomes stationary after a several collisions per particle as is expected from the kinetic theories NE98 ; BP00 ; HOB00 . These “stationary” values of $`a_2`$ agree very well with the results of kinetic theory for the nearly elastic region $`r0.95`$ (Fig.4).

This form of distribution, however, is not really stationary as it may look in the initial stage data of Fig.3. The $`\tau `$ dependence of Sonine coefficients over longer time scale is shown in Fig.5 for $`r=0.9`$, 0.95, and 0.98. In this time scale, the plateau is hardly seen and absolute value of all the coefficient show gradual decrease towards zero till the time when the clustering instability sets in; after that time the distribution deviates from Gaussian drastically.

The Sonine coefficient $`a_2(\tau )`$ for various $`r`$ is plotted in Fig.6 where the time is scaled by the clustering time $`\tau ^{}`$ when the clustering instability sets in and $`a_2`$ is scaled by its maximum absolute value $`|a_2^{}|`$ for each $`r`$. The closer the value of $`r`$ is to 1, the smaller the slope in the $`\tau `$ dependence becomes, but for all cases, the gradual return to Gaussian starts almost immediately after the initial transient period finishes. It starts actually far before any instability becomes evident.

This behavior obviously contradicts to the results of the kinetic theories based on the Boltzmann equation NE98 ; BP00 ; HOB00 ; the theories predict the distribution shows the stationary form after the short initial transient. The stationary form should last till the Boltzmann equation becomes invalid due to correlations developed in the system. The fact that the distribution function starts to deviate from the stationary form quite early stage suggests that the correlation due to inelasticity becomes important much earlier than it is generally expected com1 .

The correlation that is responsible to the behavior of the velocity distribution is the velocity correlation. This can be seen by examining the behavior of the system where the particle velocity is artificially randomized by the operation that the velocity of each particle is shuffled by exchanging them between pairs chosen randomly. By doing this, we destroy the spatial correlation of velocity while preserving the velocity distribution. The dashed line in Fig.1 shows the energy decay with the velocity shuffling and indicates that the clustering is prevented by the velocity shuffle. From Fig.7, we can see clearly that this system shows the stationary form of the velocity distribution whose Sonine coefficients are closer to those predicted by the theories.

In Fig.7, the density dependence of the Sonine coefficient behavior is also examined for $`n=`$0.25, 0.111, and 0.0625. The general tendency that $`|a_2|`$ decreases toward zero after the initial deviation is the same and they are all clearly different from the case with the velocity shuffling although a certain density dependence exists in the present density regioncom2 .

The kinetic theories have succeeded in explaining many aspects of granular gas, but the Boltzmann-Enskog equation, on which most of the theories based, ignores the particle correlations except for the pair correlation factor of the position. This approximation is quite good in the equilibrium system thanks to the absence of the velocity correlation. In the inelastic systems, however, the system develops the velocity correlation, and the assumption of molecular chaos fails even at very early stage, where the system is still in HCS. This invalidates the prediction based on the Boltzmann equation that the functional form of the velocity distribution of the inelastic system in HCS becomes stationary.

In summary, using large scale MD simulations, I have demonstrated that in the inelastic system the velocity distribution does not stay in a stationary form, contrary to the expectation by the kinetic theories based on the Boltzmann equation. This is due to the velocity correlation developed through the inelastic collisions, and this effect manifests itself in the velocity distribution from very early stage where any instabilities caused by the inelasticity are still hardly visible.