text:
# Hexagon Formation as a Compromise Between Squares and Stripes in Vertically Vibrated Granular Layers

## I Introduction

Spatiotemporal patterns are found in numerous physical, chemical, and biological systems, and have drawn much interests in recent years . Like fluid systems (Faraday systems) , familiar but poorly understood granular systems  show various patterns including stripes, squares, hexagons and kinks when vertically vibrated . Novel two-dimensional localized structures, termed ’oscillons’, were also observed in these systems . Since the pattern formation phenomena in vertically vibrated granular systems were first reported, many models have been proposed to explain the mechanism . Molecular dynamics simulations  quantitatively reproduce every observed patterns except oscillons and kinks, but do not give the essentials of the mechanism. Phenomenological models  explain some of the phenomena by considering the grain transfer with a few qualitative features. Especially, the model in Ref.  reproduces stripes, squares, and hexagons by using two simple features, i.e., randomizing impacts and inelastic collisions. The continuum coupled map model  yields a phase diagram similar to that of the experiments with the particular choice of the map and the introduction of a second length scale and explains the universal bifurcation sequence related to the interaction between temporal period doubling and spatial pattern formation. The results obtained by using the order parameter equation  indicate that the localized structures are not specific to granular systems . Although these models give the key to the mechanism, hexagon formation has not been properly treated. The models of Refs.  exhibit hexagons, but do not elucidate the hexagon forming conditions. Other models reproduce only squares  or only stripes  besides oscillons.

In this paper, by incorporating the grain transfer  qualitatively but explicitly, we propose a model that can cover all patterns including squares, stripes, hexagons, and oscillons. Hexagon formation will be discussed in connection with squares and stripes.

Previously, we proposed a simple model  which gave some insights into the dynamics of stripes, oscillons, and kinks. The following two features were assumed to study stripes and oscillons :

(i) A local excitation organizes only through interaction with excitations within some distance;

(ii) The effective interaction or flow is approximately determined by the difference between the excitations.

Additionally, a period-doubling bifurcation nature of the local dynamics was needed to investigate the kinks. With these features and a Poincaré map-like idea , we simplified the dynamics of the layer and constructed a model which maps the pattern at the time of the layer-plate collision to the pattern at the time of the next collision.

However, room exists for more generalizations of the model, for it does not show squares or hexagons yet. In this regard, some authors  argue that mass conservation has an important role in pattern formation and grain mobility contributes to the selection of squares and stripes. In experiments , squares arise at low forcing frequencies, and at those frequencies, the horizontal mobility is large due to the large dilation of the layer and the small dissipation. However, stripes are observed at high frequencies, where the mobility is small due to the small dilation and the large dissipation.

In Sec. II, we propose a modified model by introducing mass conservation and grain mobility transitions. The results from numerical simulations are presented and discussed in Sec. III.

## II Model

In this modified model, we use densities instead of local heights to consider the grain transfer explicitly. The density $`\sigma _n(\stackrel{}{r})`$ denotes how many grains are piled up on a position $`\stackrel{}{r}`$ in a two-dimensional space at the $`n`$th time step. As mentioned above, the $`n`$th time step corresponds to the $`n`$th collision between the layer and the plate. Note that in the case of $`f/2`$ patterns, the densities behave as the local heights of the patterns we observe.

This model maps the densities at the $`n`$th time step to those at the $`(n+1)`$th time step through following two stages:

$`Thefirststage:\sigma _{}^{}{}_{n+1}{}^{}(\stackrel{}{r})`$ $`=`$ $`\sigma _n(\stackrel{}{r})+\mathrm{\Delta }\sigma (\stackrel{}{r}),`$ (1)
$`Thesecondstage:_t\sigma (\stackrel{}{r},t)`$ $`=`$ $`D_{}^{}{}_{}{}^{2}\sigma (\stackrel{}{r},t),`$ (3)
$`\mathrm{or}\mathrm{equivalently},`$ (4)
$`\overline{\sigma }_{n+1}(\stackrel{}{k})`$ $`=`$ $`e^{Dt_dk^2}\overline{\sigma ^{}}_{n+1}(\stackrel{}{k}),`$ (5)
where $`\mathrm{\Delta }\sigma (\stackrel{}{r})`$ is the increment or the decrement in the density at $`\stackrel{}{r}`$, $`D`$ is a diffusion coefficient, $`t_d`$ is the effective time duration of diffusion, and $`_{}^{}{}_{}{}^{2}=_x^2+_y^2`$. In Eq. (5), $`\overline{\sigma }_n(\stackrel{}{k})`$ denotes the spatial Fourier transform of $`\sigma _n(\stackrel{}{r})`$.

The two stages are similar to those of the model proposed in Ref. , but here we can control the grain mobility. Also, as our original model , this model shares a common spirit with the model of Ref. : these models are spatially continuous, but temporally discrete.

The First Stage.—— Since there is a limit on the distance a grain can travel during the time interval between the layer-plate collisions, we take the inside of a circle with radius $`R`$ centered at $`\stackrel{}{r}`$ in two-dimensional space as an interaction region. Some of the grains at this position $`\stackrel{}{r}`$ will flow into the interaction region $`𝒩(\stackrel{}{r})=\left\{\stackrel{}{r}^{}\right|0<|\stackrel{}{r}^{}\stackrel{}{r}|R\}`$. The radius $`R`$ can be determined by the horizontal velocity, the free flight-time, and the collisions with other grains. For simplicity, we set $`R`$ to the same constant for all $`\stackrel{}{r}`$’s.

The number of grains transferred from one position to another depends largely on the lateral velocity of the grains and the free-flight time during which the lateral transfer of grains occurs . The lateral velocity is approximately proportional to both the layer-plate collision velocity and the slope of the free surface of the layer ; the free flight-time is roughly proportional to $`1/f`$.

From these facts, we consider the quantities

$`\delta _n(\stackrel{}{r}\stackrel{}{r}^{})`$ (6)
$`={\displaystyle \frac{\alpha }{S(𝒩(\stackrel{}{r}))}}G\left(\sigma _n(\stackrel{}{r})\sigma _n(\stackrel{}{r}^{})\right)M_\epsilon \left({\displaystyle \frac{S(𝒩_<(\stackrel{}{r}))}{S(𝒩(\stackrel{}{r}))}}\right),`$ (7)
$$\mathrm{\Sigma }_n(\stackrel{}{r})=_{𝒩_<(\stackrel{}{r})}𝑑\stackrel{}{r}^{}\delta _n(\stackrel{}{r}\stackrel{}{r}^{}),\stackrel{}{r}^{}𝒩_<(\stackrel{}{r}),$$
(8)
where $`𝒩_<(\stackrel{}{r})=\{\stackrel{}{r}^{}𝒩(\stackrel{}{r})|\sigma _n(\stackrel{}{r}^{})<\sigma _n(\stackrel{}{r})\}`$. The quantity $`\delta _n(\stackrel{}{r}\stackrel{}{r}^{})`$ in Eq. (7) is the maximum number of grains that can be transferred from $`\stackrel{}{r}`$ to $`\stackrel{}{r}^{}`$ where the density is lower than $`\stackrel{}{r}`$ if there are enough grains at position $`\stackrel{}{r}`$. We call this maximum number the demand from $`\stackrel{}{r}^{}`$ to $`\stackrel{}{r}`$. The actual number transferred is limited by the density at $`\stackrel{}{r}`$. For this purpose, $`\mathrm{\Sigma }_n(\stackrel{}{r})`$, the total demand to $`\stackrel{}{r}`$, is given in Eq. (8). In this model, the parameter $`\alpha `$ acts like the layer-plate collision velocity (or the relative collision acceleration) which determines the onset of patterns . The dependence of the transferred number on the slope of the free surface is phenomenologically represented by the function $`G`$, which depends only on the difference between the densities at the two positions. Due to the avalanche-type flow , the number of grains transferred is naturally expected to increase faster than linearly as the difference increases, so we choose a slope-increasing function for $`G`$ \[$`G^{\prime \prime }>0`$; see Fig. 1 (a)\]. As in our previous model , this slope-increasing property of the function $`G`$ makes the hysteresis more evident and strong: As patterns grow from the small fluctuations in the densities, the density differences increase more, and the effect of this coupling function becomes more important. Therefore, this effect makes the quantities in Eqs. (7) and (8) large enough to support the excitations even if $`\alpha `$ is lowered below the critical values where the transitions from flat to patterns occur (see Fig. 4 to be discussed). The function $`M_\epsilon `$ controlling the mobility of the grains at a position $`\stackrel{}{r}`$ is a function of the status of the position $`\stackrel{}{r}`$ in the interaction region $`𝒩(\stackrel{}{r})`$. $`\epsilon `$ is a control parameter. The status of a position $`\stackrel{}{r}`$ is given by the quantity $`S(𝒩_<(\stackrel{}{r}))/S(𝒩(\stackrel{}{r}))`$, where $`S(𝒩_<(\stackrel{}{r}))`$ is the area with densities lower than $`\sigma _n(\stackrel{}{r})`$ in the interaction region $`𝒩(\stackrel{}{r})`$, and $`S(𝒩(\stackrel{}{r}))`$ is the total area of the interaction region. Thus, positions with relatively higher densities have higher status. For example, a peak in the layer corresponds to the highest status and a valley to the lowest one. Experimentally, there are qualitative changes in the overall mobility as the forcing frequency changes. The mobility is high at low frequencies while at high frequencies the grains are less mobile. Without $`M_\epsilon `$, in this model, the grains at positions with high status are more apt to move than those at positions with low status due to the large values of the quantities $`\delta _n(\stackrel{}{r}\stackrel{}{r}^{})`$ and $`\mathrm{\Sigma }_n(\stackrel{}{r})`$. We can control the relative mobilities between these positions by giving different $`M_\epsilon `$ values to them according to their status. If we give high $`M_\epsilon `$ values to positions with low status \[see Fig. 1 (b), $`\epsilon <0`$\], the mobilities at these positions increase. In such a case, most of the grains are apt to move. Contrary to this case, if we give low $`M_\epsilon `$ values to positions with low status \[see Fig. 1 (b), $`\epsilon >0`$\], the grains at these positions are even more immobile. The parameter $`\epsilon `$ in our model roughly plays a role analogous to that of the forcing frequency in the experiments. The denominator $`S(𝒩(\stackrel{}{r}))`$ in Eq. (7) is introduced to diminish the $`\alpha `$ value dependence on the radius of the interaction region $`𝒩(\stackrel{}{r})`$.

Since there are limited number of grains at the position $`\stackrel{}{r}`$ and we want to conserve the mass, satisfying equally the demands from $`\stackrel{}{r}^{}`$, we distribute the limited number of grains to $`𝒩_<(\stackrel{}{r})`$ in two possible ways:

Case 1: If the total demand $`\mathrm{\Sigma }_n(\stackrel{}{r})`$ is less than the density at $`\stackrel{}{r}`$, the outflow from $`\stackrel{}{r}`$ to $`\stackrel{}{r}^{}𝒩_<(\stackrel{}{r})`$ is set to the demand $`\delta _n(\stackrel{}{r}\stackrel{}{r}^{})`$.

Case 2: If the total demand $`\mathrm{\Sigma }_n(\stackrel{}{r})`$ is greater than the density at $`\stackrel{}{r}`$, the total outflow from $`\stackrel{}{r}`$ is set to the density at $`\stackrel{}{r}`$, and the outflow is distributed to $`𝒩_<(\stackrel{}{r})`$ in proportion to the demand $`\delta _n(\stackrel{}{r}\stackrel{}{r}^{})`$.

Note that grains are more actively transferred in Case 2 than in Case 1, for in such a case, all the grains at those positions leave their positions.

In this manner, we can calculate the inflow and the outflow of the grains for each position at this first stage.

The Second Stage.—— In the second stage, to include the diffusion effect due to the collisions between grains and between the layer and the plate , we apply Eq. (3) which automatically conserves mass. This process relaxes the high-density gradient occurring within the size of the characteristic length scale $`\lambda _0=2R`$. Since the net diffusion effect is determined by the parameter $`Dt_d`$, we choose $`Dt_d(\frac{2\pi }{\lambda _0})^2`$ such that the Fourier components of the density distribution with wave vectors greater than $`k_0=\frac{2\pi }{\lambda _0}`$ are sufficiently suppressed.

## III Numerical Results and Discussion

For the numerical simulations, we choose the functions (see Fig. 1)

$`G(x)`$ $`=`$ $`x+cx^2,c>0,x>0,`$ (9)
$`M_\epsilon (x)`$ $`=`$ $`x^\epsilon ,0<x1.`$ (10)
These functions are the simplest ones with the required properties. $`G(x)`$ in Eq. (9) is a monotonically increasing function with a slope-increasing property ($`G^{\prime \prime }>0`$). $`M_\epsilon `$ in Eq. (10) is monotonically increasing for $`\epsilon >0`$, but monotonically decreasing for $`\epsilon <0`$. Other functions with the same properties do not give qualitatively different results. We fix $`c=0.7`$ in the following simulations for the stability of oscillon structures. Higher (lower) $`c`$ values give larger (smaller) hysteresis.

We performed numerical simulations on a $`128\times 128`$ grid with spatial extension $`L=128`$ and periodic boundary conditions. Initially, the densities are distributed around the mean value $`\sigma _0=10`$ with small uniform random deviations within $`[10^3,10^3]`$.

Square and stripe patterns obtained using this model are shown in Figs. 2 and 3. These patterns oscillate with $`f/2`$ as in experiments ; that is, they return to the same configuration after 2 periods of the forcing (in this case, one time step corresponds to one period). Figure 4 shows the stability region for each pattern as a function of $`\alpha `$ and $`\epsilon `$. There are two transition regions between squares and stripes. At lower $`\alpha `$’s, the transition from squares to stripes occurs as the parameter $`\epsilon `$ increases from negative to positive values, and for positive $`\epsilon `$’s, the transition from stripes to squares occurs as $`\alpha `$ increases. Both of these transitions correspond to qualitative changes in the grain mobility. As expected, squares arise for $`\epsilon <0`$ or sufficiently large $`\alpha `$, where the quantities in Eqs. (7) and (8) are large enough that grains are actively transferred: All the grains at positions of more than $`85\%`$ of the entire space leave their positions according to Case 2 at the first stage of our model. On the contrary, when stripes appear, the percentage decreases below $`80\%`$. For the same reason, greater $`c`$ values for $`G(x)`$ in Eq. (9) lowers the transition points from stripes to squares for positive $`\epsilon `$ values.

Figure 3 also reflects the grain mobility transitions. As pointed out in Ref. , the densities are sharply clustered in the square patterns (the high-mobility situation) and more diffuse in the stripe patterns (the low-mobility situation); the peaks in the square pattern are higher than those in the stripe patterns (see Fig. 3).

In the hysteretic region, this model also exhibits period-2 oscillons and the oscillon structures such as dipoles, triangular tetramers, and chains (not shown), which are observed experimentally . When $`\alpha `$ is increased above the upper stability limit, the structures grow by nucleation of oscillons ($`\epsilon <0`$ in this model, low $`f`$ in experiments) or by elongation of oscillons and nucleation of waves ($`\epsilon >0`$ in this model, $`f>30`$ in experiments) .

Now, we consider the hexagon formation. Hexagons are also period-2 patterns like squares and stripes, but there is a difference. Two distinct phases – a set of isolated peaks on a triangular lattice and a hexagonal cell – appear in alternate cycles. Experimentally, hexagons arise when effective two-frequency forcing is applied internally or externally . In this situation, the layer has two different free-flight times and, consequently, collides with the plate alternately with two different colliding conditions. In the absence of any clear identification of the colliding conditions , from the two phases of hexagons, we conjecture that one is a high-mobility condition, and the other is a low-mobility condition. To determine the validity of our conjecture, we apply a high-mobility condition (a square-forming condition) in one step and in the next step a low-mobility condition (a stripe-forming condition) with the same radius. This yields hexagonal patterns as shown in Figs. 5 and 6.

As in experimental results , hexagonal patterns obtained in this model show two different phases - a peak phase and a cellular phase (Fig. 5). Note that the peak phase has localized density peaks like squares and always appears after applying the high-mobility condition \[see Fig. 3(a), Fig. 6(a), Fig. 5(a)\]. Similarly, note also that the cellular phase has ridges like stripes and always arises after applying the low-mobility condition \[see Fig. 3(b), Fig. 6(b), Fig. 5(b)\]. From these facts, we can regard the hexagon formation as a compromise between squares and stripes. Our results are consistent with the occurrence of hexagons in the transition region between squares and stripes in Fig. 2 of Ref. . However, the hexagon formation described by the order parameter equation  is different from ours. In our model, grain transfer is crucial for hexagon formation, but in the case of Ref.  grain transfer is treated irrelevant to hexagon formation.

Since the hexagons in this model are formed directly from the flat state and since the transitions are hysteretic, we have parameter regions with coexistence of hexagons and flat. In that region, we obtain stable localized hexagonal structures as shown in Fig. 7. These structures may be regarded as oscillon structures. Actually, if the parameter $`\alpha `$ is raised above the upper stability limit, these hexagonal structures grow into a global hexagonal pattern by nucleating oscillons of equal phase at the hexagonal lattice sites (see Fig. 8) as the oscillon structures at low forcing frequency grow through square-like structures into squares and at high frequency through stripe-like structures into stripes . Although such localized structures have not been reported in granular systems, similar structures have been found in other systems . Especially, the localized hexagonal oscillon structures found in the two-frequency forcing experiments with Newtonian fluids  seem to be closely related to the structures found in our model. The growing phenomena reported in Refs.  are also similar to the one described above. In contrast, the hexagonal patch in Ref.  grows by a six-faceted, moving front, and the three independent roll states in superposition give the hexagonal structure.

## IV conclusion

In summary, considering the grain transfer qualitatively but explicitly, we have constructed a model for pattern formation in vertically vibrated granular layers. Our model maps the density distribution at the time of the layer-plate collision to that at the time of the next collision in a continuous two-dimensional space. This model exhibits squares in high-mobility conditions and stripes in low-mobility conditions. To elucidate the hexagon forming situation, we simulated period doubling situation by applying a high-mobility condition (a square-forming condition) and a low-mobility condition (a stripe-forming condition) alternately and observed hexagons. Our simulations show that the peak phase of hexagons appears after the high-mobility condition and the cellular phase after the low-mobility condition. The peak phase has localized density peaks like squares, and the cellular phase has ridges like stripes. From these facts, we conclude that in vertically vibrated granular layers hexagon formation is a compromise between squares and stripes. Finally, the model presented here also reproduces the experimentally observed oscillon structures and their dynamics. In addition, we observed localized hexagonal structures and their growth by nucleation. Our model may be useful for the studies of multifrequency forcing experiments  and other pattern forming systems .

## acknowledgments

We thank Hwa-Kyun Park for useful discussions. This work was supported in part by the the Interdisciplinary Research Program (Grant No. 1999-2-112-002-5) of the Korea Science and Engineering Foundation and in part by the BK21 program of the Korea Ministry of Education. We also appreciate the support of the Korea Research Foundation (Grant KRF-2000-DP0097).

Figure Captions

Fig. 1. Functions in Eq. (7). (a) $`G(x)`$ is a function of the density difference between two positions. (b) $`M_\epsilon `$ is a function of the status (see text for definition). the parameter $`\epsilon `$ controls the overall grain mobility.

Fig. 2. Shadow graph of the patterns obtained numerically in our model. (White denotes high density.) (a) Period-2 square pattern at $`\alpha =3,\epsilon =3,R=6`$, and $`Dt_d=3.65`$. (b) Period-2 stripe pattern at $`\alpha =11,\epsilon =+3,R=6`$, and $`Dt_d=3.65`$.

Fig. 3. Three-dimensional perspective of the densities in (a) the square and (b) the stripe patterns of Fig. 2.

Fig. 4. Phase diagram of our model as a function of $`\alpha `$ and $`\epsilon `$ with $`R=6`$ and $`Dt_d=3.65`$. Closed (open) squares denote transitions from flats (patterns); closed triangles and circles indicate the transition points. The numbers near symbols are the corresponding $`\alpha `$ values.

Fig. 5. Hexagonal pattern obtained by alternately applying a high-mobility condition and a low-mobility condition. (a) Peak phase $`(n=odd)`$ after applying the high-mobility condition $`(\alpha _1=5,\epsilon _1=3)`$ and (b) cellular phase $`(n=even)`$ after applying the low-mobility condition $`(\alpha _2=5,\epsilon _2=+3)`$. $`R_1=R_2=6`$ and $`(Dt_d)_1=(Dt_d)_2=3.65`$.

Fig. 6. Three-dimensional perspective of the densities in the hexagonal pattern of Fig. 5 in (a) the peak phase and (b) the cellular phase.

Fig. 7. Localized hexagonal structure in the hysteresis region obtained by applying the same conditions as in Fig. 5, except for $`\alpha _1=\alpha _2=1.1`$.

Fig. 8. Growth of hexagonal structure by nucleation. The $`\alpha `$ values are raised from $`1.1`$ in Fig. 7 to $`1.5`$.


format: text

text:
# Growth of Bose-Einstein Condensates from Thermal Vapor

## Abstract

We report on a quantitative study of the growth process of 87Rb Bose-Einstein condensates. By continuous evaporative cooling we directly control the thermal cloud from which the condensate grows. We compare the experimental data with the results of a theoretical model based on quantum kinetic theory. We find quantitative agreement with theory for the situation of strong cooling, whereas in the weak cooling regime a distinctly different behaviour is found in the experiment.

The non-equilibrium path to the formation of a Bose-Einstein condensate  is one of the most intriguing physical processes in ultra cold quantum gases. When the critical temperature of this phase transition is crossed, the atoms accumulate in the ground state of the system, which leads to the characteristic dense core in the atomic distribution and to the long-range phase coherence. It is an experimental challenge to monitor this process in time, aiming for a precise understanding of the formation of Bose-Einstein condensates. Experimental studies of the growth of Bose-Einstein condensates have so far either focused on verifying the effect of bosonic stimulation  or on the case of attractive interactions between the atoms, which gives rise to special dynamics . However, the key role of the thermal cloud, from which the condensate grows, has not yet been addressed in experiments.

In a pioneering theoretical study Kagan et al.  proposed the appearance of a quasi-condensate with strong phase fluctuations prior to the formation of long-range phase coherence. Stoof  presented the idea of nucleation followed by a growth according to a kinetic equation. A simple quantitative theory has been developed by Gardiner et al.  on the basis of quantum kinetic theory. The key equation describing the time dependence of the number $`N_0`$ of atoms in the condensate reads

$$\dot{N}_0=2W^+(N_0)\left\{\left(1e^{\left[\mu _C(N_0)\mu \right]/k_BT}\right)N_0+1\right\},$$
(1)
where $`\mu _C(N_0)`$ is the chemical potential of the condensate, $`\mu `$ the chemical potential of the thermal vapor, and $`W^+(N_0)`$ an analytically known rate factor . Starting from $`N_0=0`$ the condensate starts to grow at a finite rate, which turns quickly into exponential growth, before it reaches its asymptotic equilibrium value when the condensate chemical potential approaches that of the thermal cloud. This equation describes only the growth of the population in the ground state of the system and it neglects certain scattering processes among atoms in higher lying states. Especially in the initial period the redistribution of population between higher lying states becomes important and modifies the growth . When the scattering terms are included one finds reasonable agreement with the experiment described in , but qualitative and quantitative discrepancies remain . Moreover, equation (1) neglects the dynamics of the thermal component by assuming its chemical potential to be stationary. The influence of the dynamics of the thermal cloud was included in the recent models . A purely numerical approach by Monte Carlo simulations is presented in .

We have identified and studied different stages of formation when the Bose-Einstein condensate grows from a thermal gas. In the experiment a cloud of thermal atoms is prepared above the critical temperature and then cooled into the quantum degenerate regime by continuous evaporation for a duration of up to 6 s. This strategy allows us to observe the formation process in slow motion. A distinct growth characteristic is revealed by directly controlling the properties of the thermal vapor. In comparison, previous experimental and most theoretical work addressed the situation in which a short radio-frequency pulse is applied to remove hot atoms, then turned off leaving the system to evolve toward thermal equilibrium without any further forced evaporation.

For a quantitative understanding it is crucial to prepare a well-defined initial state of the thermal cloud. As described previously , we load about $`10^9`$ atoms in the $`|F=1,m_F=1`$ hyperfine ground state into a magnetic Quadrupole-Ioffe configuration trap (QUIC-trap). The trapping frequencies are $`\omega _{}=2\pi \times 110`$ Hz in the radial and $`\omega _z=2\pi \times 14`$ Hz in the axial direction. In the magnetic trap we perform forced evaporative cooling for about 25 s by continuously lowering an RF frequency to a final value of $`\nu _{RF,0}`$=2120 kHz. This prepares a purely thermal cloud of $`N_i`$=(4.2$`\pm `$0.2)$`\times 10^6`$ atoms at a temperature of $`T_i`$=(640$`\pm `$30) nK. The minimum of the trap is determined by atom laser output coupling  and corresponds to a radio-frequency of 1955 kHz. The stability of the magnetic trap enables us to approach the formation of the condensate very slowly. The measured stability  and the good reproducibility of the experiment show that the formation of the condensate is not triggered by fluctuations of the magnetic field. When operating the trap at full current, we have measured a slow upwards drift of the bottom of the trapping potential of 5 kHz/s, which is due to thermally induced motion of the trap coils.

After this preparation of a purely thermal cloud, the atoms are cooled through the BEC phase transition. The radio-frequency field jumps abruptly to a lower frequency $`\nu _{RF,1}`$ and remains at this value for a period of up to 6 seconds. This limits the depth of the magnetic trap to $`ϵ_{cut}`$, which is expressed in units of $`T_i`$ by the evaporation parameter $`\eta =ϵ_{cut}/k_BT_i`$. Atoms with an energy larger than $`ϵ_{cut}`$ are lost from the trap due to spin flips into a non-trapped state. Equilibration of the cloud by elastic collisions leads to a decrease of temperature and the Bose-Einstein condensate forms during this process.

The time evolution of the condensation is probed by releasing the atoms from the magnetic trap after different evolution times, for otherwise identical repetitions of the experiment. The atom cloud is imaged by absorption imaging after 17 ms of ballistic expansion. Atom numbers are extracted from the absorption images by two-dimensional fits to the atomic density distribution. The thermal cloud is fitted by a Bose-Einstein distribution under the assumption of zero chemical potential, and the condensed fraction is fitted by a Thomas-Fermi distribution. The error in the determination of the atom number is estimated to be below 10%. Here, the ”condensate number” refers to the integral under the Thomas-Fermi profile. Temperatures are determined from a Gaussian fit to the far out wings of the thermal atom cloud. The statistical error on the temperature determination is 5%, whereas the systematic error is 10-15%.

The initially prepared atomic sample corresponds to $`\eta =6`$, where the gravitational sag of the atom cloud must be taken into account. Gravity pulls the atomic cloud away from the center of the magnetic trap and thereby reduces the dimensionality and the efficiency of evaporative cooling. We have performed experiments in the range $`0.75<\eta <4.8`$, corresponding to RF frequencies between 1990 kHz and 2095 kHz. For $`\eta >4.6`$ we could not obtain Bose-Einstein-condensates within the 6 s period of the measurement due to the low efficiency of the cooling. For $`\eta =0.75`$ we did not observe the formation of a condensate since the atom number of the initial thermal cloud was reduced too much for evaporative cooling to proceed.

Figure 1 shows the measured data for an evaporation parameter $`\eta =1.4`$. Every data point is an average over three measurements and the error bars show the statistical error. The final temperature of the cloud is (220$`\pm `$20) nK and the final condensate fraction is ($`55\pm 7`$)%. The solid line is a numerical calculation of the condensate growth according to a model based on quantum kinetic theory, where the influence of the thermal atoms is taken into account . The model was improved to include the continuous evaporation of atoms, the gravitational sag of the atomic cloud, the experimentally observed magnetic field drift, and three-body loss of atoms from the condensate . We find quantitative agreement between the model and the measured data without any free parameters entering into the calculation, which suggests that the model describes this regime of the formation of Bose-Einstein condensates sufficiently. The best agreement between data and theory has been obtained for the calculation with the starting conditions $`N_i=4.4\times 10^6`$ and $`T_i=610`$ nK or $`N_i=4.2\times 10^6`$ and $`T_i=610`$ nK, both of which lie within the range of the experimental uncertainties.

We have compared the time of the onset of Bose-Einstein condensation after the beginning of the final radio frequency cooling with the calculated curves and find excellent agreement for the situation of strong cooling (see Fig. 2). To quantify the initiation time $`t_i`$ for the measured data we have fitted the function

$$N_0(t)=N_{0,f}\left(1e^{\mathrm{\Gamma }(tt_i)}\right),$$
(2)
which describes the final stage of the condensate growth towards equilibrium (see Fig. 3). $`\mathrm{\Gamma }`$ describes a relaxation rate towards equilibrium. This functional form was suggested in Ref.  and subsequently found to be a reasonable approximation to more sophisticated calculations . To account for the range of validity of the equation we have restricted the fitting range to condensate numbers larger than 30% of $`N_{0,f}`$, where we have obtained good agreement with the data. For low values of the evaporation parameter, i.e. rapid cooling, the initiation time approaches a minimum of 200 ms. For the example shown in figure 1 we obtain an initiation time of $`t_i=(215\pm 10)`$ ms. This time corresponds to $`30\tau _0`$, where $`\tau _0=(n\sigma v)^1`$ is the classical collision time in a thermal gas ($`n`$: mean atomic density, $`\sigma `$: scattering cross section, $`v`$: average thermal velocity). For the homogeneous Bose gas the first stage of formation was predicted to be on the order of $`\tau _0`$ . Early numerical calculations have confirmed that this should also be valid for the case of a harmonic confining potential, provided that all scattering processes are included . Most theoretical treatments require sufficient ergodicity, meaning that ergodic mixing occurs on the order of a few collision times. At temperatures below a few $`\mu `$K this does not describe the experimental situation. Evaporation in the magnetic trap becomes one-dimensional due to the gravitational sag of the atom cloud. This leads to a significantly lower evaporation rate  and the ergodic mixing times can be on the order of several ten collision times .

The measured initiation times show a pronounced divergence at $`\eta 4.6`$ and become much longer than the ergodic mixing time of some ten $`\tau _0`$. Very slow cooling through the phase transition does not only affect the initiation time but it also changes the characteristics of the growth curve. We detected a slow, approximately linear growth of a condensate prior to the initiation time, determined from equation (2). The open circles in figure 2 indicate when a condensate is observed for the first time and the gray stars show the time the condensate population exceeds $`10^4`$ as obtained from theoretical simulations of the growth experiments .

In figure 3 the two distinct stages of condensate formation can be seen for the case of slow cooling ($`\eta =4.6`$): first, there is a approximately linear growth at constant rate of $`\dot{N}_0=(44\pm 8)\times 10^3`$ s-1, followed by a very rapid growth with an initial rate 8 times larger. The slow linear growth of the condensate indicates a less pronounced bosonic stimulation during this stage of the condensate formation. The numerical calculation does not reproduce the measured condensate growth data accurately in this regime. The calculation rather shows a growth rate of the same order of magnitude as the second stage of the observed growth, suggesting that the model works well if the system is sufficiently far below the transition temperature. For the thermal cloud the calculated temperatures agree with the experimentally determined values. We deduce that a possible heating rate of the atoms in the trap is below our sensitivity of 5 nK/s and thus probably negligible. The measured decrease in the total atom number for long time scales seems to be slightly larger than expected from the simulations. Whereas three body losses in the condensate are included in the calculations, we neglect three body recombination involving thermal atoms. The factor of 6 enhancement for these processes is negated by the much lower density of the thermal cloud cloud: for the given trap parameters and a temperature of 300 nK one obtains a lifetime of more than 40 s . Similarly, the background limited lifetime of the magnetic trap is at least 40 s.

We attribute the discrepancy between the theoretical model and the measured condensate numbers for slow cooling to approximations in the theoretical model which for example do not take into account coherence effects in the scattering process between the atoms. Possibly, the two-stage growth can be explained by the formation of quasi-condensates which have been predicted as an early stage in the formation process of a Bose-Einstein condensate . Quasi-condensates exhibit strong phase fluctuations that die out as the system grows. These phase fluctuations are due to a population of low-lying excited states in the trap. Density fluctuations are suppressed due to their comparatively large excitation energy and therefore the density exhibits a Thomas-Fermi distribution. The growth rate could then be significantly smaller than expected from the assumption of Bose-stimulation. Indeed, we find the growth rate to be almost one order of magnitude smaller than calculated from our model based on quantum kinetic theory. The time scale for both stages of the evolution is longer than the ergodic mixing time. The system is closer to equilibrium compared to the case of rapid cooling, i.e. small values of the evaporation parameter $`\eta `$. However, the kink in the growth curves clearly shows that the condensate is not formed adiabatically. In this case the atom number should increase according to $`N_0=N(1(T/T_c)^3)`$ when the temperature is lowered smoothly.

It was recently found that 3D Bose-Einstein condensates can exhibit strong phase fluctuations in elongated traps . The phase fluctuations are characterized by the parameter $`\delta _L^2=(T/T_c)(N/N_0)^{3/5}\delta _c^2`$, where $`\delta _c^2`$ is a factor that mainly depends on the trap geometry and atom number. Its value is about 0.15 for our experiment. The kink in the growth curves is observed at temperatures higher than 0.9 $`T_c`$, thus large phase fluctuations $`\delta _L^2>1`$ can be expected in the regime $`N_0/N<4\%`$. We find the transition from linear growth to relaxational growth for slow cooling always in the range $`1.5\%<N_0/N<3.5\%`$, which is compatible with the theoretical expectation. When the phase fluctuations have died out, i.e. $`\delta _L^21`$, a true Bose-Einstein condensate forms. This stage of the evolution then shows the same growth characteristics as observed for rapid cooling. So far calculation of the kinetics of the growth process have not taken into account the large aspect ratio of cigar shaped traps. Therefore axial phase fluctuations may be not included in these models. Evidence for axial phase fluctuations in a steady state of a very elongated condensate is reported in .

In conclusion, we have investigated the influence of the thermal atom cloud on the growth of a Bose-Einstein condensate. The measured data has been compared to a numerical model based on quantum kinetic theory. We find the minimum initiation time for condensate growth to be 30 times the classical collision time, in quantitative agreement with our calculations. This comparatively long timescale is presumably due to slow ergodic mixing of the trapped atomic cloud. Very slow cooling of the cloud enables us to observe a two stage process during the formation of a Bose-Einstein condensate, which is possibly beyond the approximations made in the model.

We would like to thank I. Bloch for help with the data analysis software, W. Eckardt for assistance with the laser system, Yu. Kagan, W. Ketterle, and W. Zwerger for helpful comments and DFG for financial support of the experimental work. MJD was supported by the UK EPSRC, CWG was supported by the Royal Society of New Zealand under the Marsden Fund Contract No. PVT-902.


text:
# 1 INTRODUCTION

## 1 INTRODUCTION

The axisymmetric structures, e.g., elliptical or bipolar, of most planetary nebulae (PNe), has led many people to suggest that stellar magnetic fields shape the winds from asymptotic giant branch (AGB) star progenitors of PNe (e.g., Pascoli 1997; Chevalier & Luo 1994; Garcia-Segura 1997; Garcia-Segura et al. 1999; Matt et al. 2000; Blackman et al. 2001). Common to all these models and scenarios, hereafter termed dynamic-magnetic models, is the dynamic role attributed to the magnetic field. The different models are not identical in all their ingredients. For example, in some models the shaping occurs close to the stellar surface, through magnetic tension and/or pressure (Pascoli 1997; Matt et al. 2000), while in others (Chevalier & Luo 1994; Garcia-Segura 1997) the shaping occurs in the nebula at large distances from the star. Despite these and other differences the models seem to suffer from the same basic problems, which, as we argue in the present paper, prevent any of them from being the correct model for the shaping of PNe. Following earlier papers (Soker & Harpaz 1992, 1999; Soker 1998) we argue in section 2 of the present paper that magnetic fields are not likely to play a dynamic role in shaping most PNe. Only in a minority of cases, where the progenitor AGB star was substantially spun-up via a common envelope or tidal interaction, it is possible, although not necessary, that the magnetic field had a dynamic role in shaping the descendant PN. Our objections to these models, some raised in earlier papers and some new, are summarized in section 2.

One of the problems we find in some of the papers cited above is that they scale properties from the solar model to AGB stars, e.g., they assume that magnetic field can shape the intensive AGB wind close to the stellar surface. To emphasize the differences between the Sun and AGB stars regarding magnetic activity and winds, we devote section 3 to a detailed comparison of the relevant physical parameters between the Sun and AGB stars. These differences cannot be ignored in a consideration of magnetic activity in AGB and post-AGB stars.

In section 4 we show that magnetic activity can indeed take place in AGB stars, but to a lesser degree than that required by the models cited above. We argue that the magnetic field is most likely amplified by a turbulent dynamo, an $`\alpha ^2\omega `$ dynamo where the main role of the rotation is fixing a symmetry axis, rather than by a solar type dynamo, the $`\alpha \omega `$ dynamo where the rotation plays a crucial role in amplifying the toroidal component of the magnetic field. Our estimate of the AGB magnetic activity is based on new papers by Brandenburg and collaborators (Brandenburg 2001, Brandenburg, Bigazzi, & Subramanian 2001, and Brandenburg & Dobler 2001, hereafter B2001, BBS, and BD01, respectively). As argued in earlier papers (e.g., Soker 1998; 2000) the magnetic field may become dynamically important at specific locations near the surface where it forms cool spots, which can regulate the mass loss process via the formation of dust above these cool spots. However, the average magnetic energy density is much below the thermal or kinetic energy density in the wind, hence the magnetic field has no direct influence on the mass loss process. The much weaker magnetic activity than that required in the dynamic-magnetic models is sufficient to explain the following observations of magnetic fields in cool giant stars. Kemball & Diamond (1997) detected a magnetic field in the extended atmosphere of the Mira variable TX Cam, with $`B<5G`$ at the locations of SiO maser emission at a radius of $`4.8\mathrm{AU}2\mathrm{R}_{}`$. Magnetic field of $`1\mathrm{m}\mathrm{G}`$ is detected in the OH maser emission regions of U Herculis (Palen & Fix 2000). The detection of X-ray emission from a few M giants (Hünsch et al. 1998; see also Schröder, Hünsch, & Schmitt 1998) also hints at the presence of magnetic fields in giant stars, but weak. We summarize our main results in section 5.

## 2 PROBLEMS WITH DYNAMIC-MAGNETIC MODELS

In the present section we list the fundamental problems we find in models that attribute a dynamic role to the magnetic field.

### 2.1 Too fast rotation and/or unrealistic angular momentum distribution

In earlier papers (Soker & Harpaz 1992; Soker 1998) it was shown that dynamic-magnetic models must incorporate a binary companion to spin-up the envelope, since single stars slow down markedly on the AGB (Soker 2001). The model of Pascoli (1997) was criticized by Soker (1998; 2001). The model proposed by Chevalier & Luo (1994) and extended by Garcia-Segura (1997; see also Garcia-Segura et al. 1999) was shown by Soker (1998) to require a binary companion. Basically, this model is based on the tension of the toroidal component of the magnetic field in the wind during the transition from the AGB to the PN phase at large distances from the star. In contrast to other dynamic-magnetic models, in this one close to the star the magnetic pressure and tension are negligible compared with the ram pressure and thermal pressure of the wind. Only when the wind hits the outer PN shell, which is the remnant of the slow wind, and goes through a shock and slows down, does the toroidal component becomes important and shape the nebula.

The recent paper by Blackman et al. (2001) seems to suffer the most from this problem, in that it assumes an unrealistic angular momentum distribution. Blackman et al. (2001) propose that the magnetic field is amplified via a solar type $`\alpha \omega `$ dynamo in the core-envelope interface (an amplification of the magnetic field close to the AGB core was already suggested by Pascoli 1997). For an $`\alpha \omega `$ dynamo a relatively large angular velocity gradient is required, which they assume is because each mass shell conserves its angular momentum from the main sequence up to the upper AGB. This assumption seems unrealistic, since a strong coupling is expected in the convective AGB envelope. Even in the radiative core the powerful weak-field MHD instability (Balbus & Hawley 1994) is likely to force a solid body rotation. The strong magnetic field obtained by Blackman et al. (2001) will force a very fast solid body rotation, as we now show, such that the uniformly rotating envelope will slow down very fast with mass loss (Soker 2001).

A curved magnetic flux tube embedded inside the envelope exerts a tension-force per unit volume given by

$$f=\frac{B^2}{4\pi R_c},$$
(1)
where $`B`$ is the magnitude of the magnetic field and $`R_c`$ the radius of curvature of the magnetic flux tube. At one point the flux tube rises to a height $`h`$, with an upward speed $`v_u`$, during a time $`t`$ given by $`th/v_u`$ (we ignore the unlikely case where the flux tube is exactly circular with the same radial positions around its entire circumference). The differential rotation results in an average (over the raising time $`t`$) relative azimuthal velocity between the upper and lower parts of the tube given by $`v_y\left(rd\omega /dr\right)h/2`$, where $`d\omega /dr`$ is the angular velocity gradient. The flux tube will be azimuthally bent along a distance given by $`yv_yt0.5\left(rd\omega /dr\right)ht`$, so that the radius of curvature is

$$R_c\frac{h^2}{2y}\frac{v_u}{r}\left(\frac{d\omega }{dr}\right)^1.$$
(2)
The azimuthal force resulting from the curved magnetic flux tube exerts a moment per unit volume of $`rf`$, where $`f`$ is given by equation (1), and tends to bring the envelope to a uniform rotation. Equating this moment to the rate of change of angular momentum per unit volume $`\rho r^2\left(d\omega /dt\right)`$, and assuming that the average magnetic flux tubes filling factor is $`\beta `$, gives the time required to bring the envelope to a uniform rotation

$$\tau w\left(\frac{d\omega }{dt}\right)^1\frac{v_ur}{v_A^2}\left(\frac{d\mathrm{ln}\omega }{d\mathrm{ln}r}\right)^1\beta ^1,$$
(3)
where $`v_A=\left[B^2/\left(4\pi \rho \right)\right]^{1/2}`$ is the Alfven velocity. The rising speed of a flux tube is some fraction of the sound speed, which Blackman et al. take to be the Alfven speed $`v_u=v_A`$. Substituting typical values used by Blackman et al. in their magnetic field amplification zone, we find

$$\tau 5\times 10^3\left(\frac{r}{R_{}}\right)\left(\frac{v_A}{4\mathrm{km}\mathrm{s}^1}\right)^1\left(\frac{d\mathrm{ln}\omega }{d\mathrm{ln}r}\right)^1\beta ^1\mathrm{yrs}.$$
(4)
For the magnetic field to play a significant dynamic role, as required by the model of Blackman et al., the filling factor $`\beta `$ cannot be too small, i.e., $`\beta 10^3`$, and since they have $`d\mathrm{ln}\omega /d\mathrm{ln}r1`$, we find the slowing-down time, which is about the time required to force a uniform rotation, to be $`\tau 10\mathrm{yrs}`$. As mentioned earlier, the powerful weak-field MHD instability which operates for much weaker magnetic fields (Balbus & Hawley 1994) is likely to force a solid body rotation as well. We therefore conclude that their assumed angular velocity profile is unrealistic.

In this regard it is not clear why Blackman et al. assume that the remnant WD can slow down by the interaction of its magnetic field with the wind, while it will not slow down (according to their assumption) by the interaction of the magnetic field with the much denser envelope during the AGB phase.

### 2.2 A Too Low Density Contrast and the Transition to Aspherical Mass Loss

In some papers there is no clear distinction between bipolar and elliptical PNe. One of the reasons is that the equatorial to polar density contrast achieved in these models is very low. To achieve even a moderate density contrast the models have to assume an unrealistically strong magnetic field (e.g., Matt et al. 2000) or extremely fast rotation (Garcia-Segura et al. 1999; see criticism by Soker & Harpaz 1999). In both these models, as well as that of Blackman et al. (2001), there is no satisfactory explanation for the observations that many PNe show a transition from an almost spherical to highly non-spherical mass loss during the late stages of the AGB and/or post-AGB. Each model has to assume ad hoc that the relevant mechanism starts to operate only near the termination of the AGB, but no satisfactory physical mechanism is proposed for the switch-on of the mechanism.

### 2.3 No Radiative Dust-Acceleration

The common view, supported both by observation and theory, is that the mechanism behind the intensive mass loss from AGB stars is radiation pressure on dust coupled with strong stellar pulsations (e.g., Wood 1979; Jura 1986; Knapp 1986; Fleischer, Gauger & Sedlmayr 1992; Habing 1996; Andersen, Loidl, & Höfner 1999). However, the papers on dynamic-magnetic activity (Matt et al. 2000; Blackman et al. 2001) omit radiation pressure altogether. For these models to work, the mass loss rate apparently has to be determined by direct magnetic activity, rather than by pulsation and radiation pressure on dust. This is in contradiction with observations.

### 2.4 A Too Strong X-Luminosity Is Predicted

In some of the dynamic-magnetic models the surface magnetic pressure is comparable to the thermal pressure of the gas (Matt et al. 2000; Blackman et al. 2001; Garcia-Segura, Lopez, & Franco 2001). This will lead to magnetic field reconnection, e.g., flares, which will cause a very strong X-ray emission. For typical values of surface magnetic fields in these models, $`B1\mathrm{G}`$, the expected X-ray luminosity is $`10^4`$ times stronger than that of the Sun, if the reconnection rate per unit surface area is similar to that in the Sun. If the X-ray luminosity is proportional to the optical luminosity the same factor holds. Note that in the Sun, where the mass loss rate is determined by magnetic activity, the average X-ray luminosity is of the same order of magnitude as the rate of kinetic energy carried by the wind. The X-ray luminosity in the ROSAT/PSPC band is in the range of $`3\times 10^{26}`$ to $`5\times 10^{27}\mathrm{erg}\mathrm{s}^1`$, at minimum and maximum, respectively (Peres et al. 2000). The solar wind’s kinetic energy falls between these values. If this is the case for AGB stars in the dynamic-magnetic models, the X-ray luminosity will be a factor of $`10^610^8`$ stronger than in the Sun. This expectation, of $`L_x10^{30}10^{35}\mathrm{erg}\mathrm{s}^1`$, in dynamic-magnetic models is in sharp contradiction with observations. From observation, the maximum X-ray luminosities of red giant stars are marginally larger than the solar X-ray luminosity (Schröder et al. 1998). In most cases $`L_x<10^{30}\mathrm{erg}\mathrm{s}^1`$, and further decreases in late giants’ evolution (Hünsch & Schröder 1996).

As mentioned in section 1, there are strong indication of magnetic fields around AGB stars. But we will argue in section 4 that these local fields (e.g., Palen & Fix 2000 for U Herculis) result from a much weaker magnetic activity.

A final comment to the entire section is our view that while the mechanism proposed by Chevalier & Luo (1994; also Garcia-Segura 1997), by Pascoli (1997), and by Garcia-Segura et al. (1999) may work, but only if the progenitor AGB star is spun-up by a stellar companion (although we still do not think these are the mechanisms for shaping most PNe), we think that the scenarios proposed by Matt et al. (2000) and Blackman et al. (2001) cannot work at all.

## 3 DIFFERENCES BETWEEN AGB STARS AND THE SUN

Soker (2000) discusses in detail a few major differences between dynamos in main-sequence stars, e.g., the Sun, and any dynamo model for AGB stars. To emphasis, clarify, and extend the list of these differences we present them in Tab1e 1. The compared variables are listed in the first column, and their symbols are given in the second column. The third column of Table 1 gives the units used, and the fourth and fifth columns give the typical values in the Sun and in upper AGB stars, respectively. The lower section of the Table gives variables relevant directly to the magnetic activity. For AGB stars these variables have to be scaled with the angular velocity $`\omega `$ and magnetic field intensity $`B`$. As representative values we take the orbital velocity to be $`0.001`$ times the Keplerian angular velocity on the equator, i.e., $`\omega =0.001\mathrm{\Omega }_3\omega _{\mathrm{kep}}`$, and for the magnetic field we take $`B=0.01B_2\mathrm{G}`$. Angular velocity larger than the above scaling requires the AGB star to be spun-up by a companion more massive than Jupiter (to an order of magnitude; see Soker 2001 for the exact values). The magnetic field is scaled according to the results of the next section.

The main relevant differences, and their implications, are as follows.
(1) In the Sun the mass loss rate is determined mainly by the magnetic activity. Therefore it is not surprising that the wind ram pressure near the solar surface is smaller than the magnetic pressure (row number 18). In AGB stars, mainly during the final intensive wind (FIW; also called superwind), the wind’s ram pressure is larger by several orders of magnitude than the magnetic field pressure, even if we take $`B=1\mathrm{G}`$. This implies that magnetic field has no dynamic role in determining the mass loss rate, and that mass leaving the star drags the magnetic field lines, rather than being dragged by the magnetic field lines.

(2) In the Sun the dynamo number is $`N_D>1`$ and the Rossby number is $`\mathrm{Ro}<1`$ (to an order of magnitude $`N_D\mathrm{Ro}^2`$) as is required by standard $`\alpha \omega `$ dynamo models. The dynamo number is the square of the ratio of the magnetic field amplification rate in the $`\alpha \omega `$ dynamo model, to the ohmic decay rate. In AGB stars, even if rotating close to the break-up velocity, the opposite inequalities hold (rows 19 and 23). As noted by Soker (2000), the low value of the dynamo number in AGB stars implies that the convective motion amplifies both the poloidal and toroidal magnetic components, i.e., $`\alpha ^2`$ dynamo. This is the subject of the next section.

(3) The fast convection motion in upper AGB stars (row 9; see Soker & Harpaz 1999 for details) means that the energy density in the convective motion is comparable to the thermal energy. In the Sun this ratio is very small. The strong convection also hints at the possible operation of a turbulent dynamo in AGB stars.

(4) As the AGB envelope gets depleted the mass loss time, defined as $`M_{\mathrm{env}}/\dot{M}`$, where $`M_{\mathrm{env}}`$ is the envelope mass, becomes shorter than the rotation period. In the Sun the mass loss time is longer by many orders of magnitude than the rotation period. The mass loss time is not given in the table since it changes over several orders of magnitude as the mass loss rate increases and envelope mass decreases toward the termination of the AGB.

(5) Some structural differences between AGB envelopes and the solar envelope are not presented here (see Soker & Harpaz 1999; note that the density scale given in their figs. 1-6 is off by a factor of 10; the correct density scaling is in their figure 6). An example is the convective region, which in AGB stars is very thick, whereas in the Sun its width is only $`0.3R_{}`$.

The main conclusion from this section, which was already mentioned in the previous section, is that the processes related to rotation, convection, mass loss, and magnetic activity cannot simply be scaled from the Sun to AGB stars. In particular, it is wrong to assume that the mass loss rate and/or geometry are dynamically dictated by the magnetic activity.

## 4 THE $`\alpha ^2\omega `$ DYNAMO IN AGB STARS

The amplification of a large-scale magnetic field in AGB envelopes via turbulent dynamo was mentioned before (Soker 2000), where simple arguments in favor of an $`\alpha ^2\omega `$ dynamo were given. The main role of rotation in the $`\alpha ^2\omega `$ dynamo is fixing a symmetry axis, rather than amplifying the azimuthal component of the magnetic field. We use recent results by Brandenburg and collaborators (B2001, BBS, and BD01), the most relevant for us being as follows. (1) The short mass loss time (see previous section) means that the dynamo has open boundaries, with significant implications for the amplification process (Blackman & Field 2000; BBS). (2) The turbulent dynamo ($`\alpha ^2`$ dynamo) can amplify a large-scale magnetic field (B2001, BD01). (3) In a homogeneous and isotropic model, even a small perturbation can lead to a preferred direction (B2001). In the present case the slow rotation of the AGB envelope fixes the symmetry axis. (4) With periodic boundary conditions the magnetic field energy reaches equipartition with the turbulent energy, or even exceeds it (B2001). (5) With open boundaries the magnetic energy density is lower by a factor of $`\eta /\eta _{\mathrm{eff}}`$ (BD01), where $`\eta `$ is the magnetic diffusion coefficient (due to processes on atomic scales), and $`\eta _{\mathrm{eff}}`$ is the effective magnetic diffusion coefficient due to convection. We now evaluate this ratio, and show that the expected magnetic field is strong enough to form magnetic cool spots on the surface of AGB stars.

The derivation of the magnetic energy density follows the results of BD01 (their $`\mathrm{\S }5.1`$) and B2001 (his $`\mathrm{\S }3.6`$). As stated above, with periodic boundary conditions the dynamo brings the magnetic energy density to an equipartition with the turbulent energy density, and even exceeds it by a factor of $`R_{}/l_T`$ (B2001 eq. 46), where $`R_{}`$ is the size of the periodic dynamo, taken here to be the stellar radius, and $`l_T`$ is the forcing scale, taken here as the mixing length, i.e., about the pressure scale height. We find that in AGB stars $`R_{}/l_T5`$, as in the calculations of B2001. Comparing eq. (23) of BD01 with equation (45) of B2001 gives the magnetic energy density for a dynamo with open boundaries $`E_{\alpha ^2}`$, in terms of the equipartition magnetic energy density $`E_e=E_T`$, where $`E_T=\rho v_T^2`$ is the turbulent kinetic energy density,

$$E_{\alpha ^2}=E_e\frac{\eta }{ϵ_H\eta _{\mathrm{eff}}}\frac{R_{}}{l_T},$$
(5)
where $`ϵ_H`$ ($`ϵ_H3`$ in the calculations of BD01) is a nondimensional quantity defined in equation (17) of BD01. The effective magnetic diffusivity is given by equation (22) of BD01. Since in AGB stars the turbulent magnetic diffusivity is much larger than the magnetic diffusivity coefficient (see Table 1), we can write to an order of magnitude

$$\eta _{\mathrm{eff}}\eta _T\frac{ϵ_Q}{2\pi ϵ_H},$$
(6)
where $`v_T=v_Tl_T`$, $`v_T`$ is turbulent velocity, and $`ϵ_Q`$ ($`\left|ϵ_Q\right|0.010.3`$ in BD01) is defined in equation (19) of BD01 as the ratio of the magnetic helicity flux to the quantity $`B^2v_T`$. The magnetic energy diffusivity coefficient due to processes on atomic scales in AGB envelopes is $`\eta 10^8\left(T_{\mathrm{env}}/10^4\mathrm{K}\right)^{3/2}\mathrm{cm}^2\mathrm{s}^1`$, where $`T_{\mathrm{env}}`$ is the envelope temperature. Substituting other typical values for upper AGB stars (e.g., Soker & Harpaz 1999), $`v_T10\mathrm{km}\mathrm{s}^1`$ and $`l_T50R_{}`$, we find

$$\frac{E_{\alpha ^2}}{E_e}10^8\left(\frac{ϵ_Q}{0.1}\right)^1\left(\frac{R_{}}{10l_T}\right)\left(\frac{T_{\mathrm{env}}}{10^4\mathrm{K}}\right)^{3/2}.$$
(7)
With the density in the outer convective region of AGB stars taken to be $`\rho =10^9\mathrm{g}\mathrm{cm}^3`$, the equipartition magnetic pressure is $`B_e=\left(8\pi E_e\right)^{1/2}=\left(8\pi \rho v_T^2\right)^{1/2}=160\mathrm{G}`$. Taking the square root of equation (7) and using this value of the equipartition magnetic field, we find the expected magnetic field intensity from an $`\alpha ^2`$ dynamo in upper AGB stars to be

$$B_{\alpha ^2}0.01G.$$
(8)
This is two orders of magnitude lower than the required magnetic field in most dynamic-magnetic models, but it is of the same order of magnitude as the magnetic field required for the formation of magnetic cool spots in the magnetic cool spots model (Soker 1998, eq. 10; note that in that equation $`\eta `$ has a different meaning).

## 5 SUMMARY

We analyzed recently proposed models which attribute the nonspherical mass loss process from upper AGB stars to strong magnetic fields, (dynamic-magnetic models) and used recent results (B2001, BBS, BD01) on the operation of a turbulent dynamo, i.e., an $`\alpha ^2`$ dynamo, and obtained the following results.
(1) We find ($`\mathrm{\S }2`$) problems in models where the magnetic field plays a dynamic role in shaping the AGB winds close to the stellar surface, i.e., the magnetic energy flux is of the same order of, or stronger than, the wind’s kinetic flux. From theoretical considerations we found that these models have to assume that the AGB envelopes rotate at very high speeds, and/or to assume unrealistic angular momentum distribution in the envelope. From observational considerations, the X-luminosity expected from these models, due to magnetic field reconnection (as in the Sun), is much higher than limits set by X-ray observations.
(2) We argue that most likely the amplification of magnetic fields in AGB stars is due to a turbulent ($`\alpha ^2`$) dynamo, where the azimuthal magnetic field component is amplified by convection, and not by differential rotation as in the solar $`\alpha \omega `$ dynamo. Only in a minority of AGB stars that have been substantially spun-up by stellar companions (Soker 2000) can an $`\alpha \omega `$ model be effective. We found by applying the recent results of B2001 and BD01 to AGB stars that the expected average magnetic field is $`B_{\alpha ^2}10^4B_e0.01\mathrm{G}`$, where $`B_e`$ is the equipartition magnetic field (where the magnetic energy density is equal to the turbulent energy density). This field intensity is indeed much below the magnetic field required by dynamic-magnetic models.
(3) B2001 found that in a homogeneous and isotropic $`\alpha ^2`$ dynamo even a small perturbation can lead to a preferred direction. We proposed that a slow rotation of the AGB envelope can fix the symmetry direction. For this role of the rotation, the dynamo in AGB stars is called an $`\alpha ^2\omega `$ (Soker 2000).
(4) Although the large-scale magnetic field is much weaker than the equipartition magnetic field, and cannot influence directly the wind from the AGB star, it is strong enough for the formation of magnetic cool spots on the AGB stellar surface (Soker 1998). As argued in earlier papers (e.g., Soker 1998; 2000) the magnetic field may become dynamically important in specific locations near the surface where it forms cool spots. The magnetic cool spots can regulate the mass loss process via the formation of dust above these cool spots, leading to axisymmetric mass loss and the formation of elliptical PNe (Soker 1998).

ACKNOWLEDGMENTS: This research was supported in part by grants from the the US-Israel Binational Science Foundation.


text:
# The Structure of Dark Matter Halos in an Annihilating Dark Matter Model

## 1 Introduction

Discrepancies between the Cold Dark Matter (CDM) model and observations on galactic scales (Salucci, 2001; Moore, 1994; Flores & Primack, 1994; Burkert, 1995; Salucci & Burkert, 2000; Moore, 2001; Sellwood & Kosowsky, 2000; Weiner et al., 2001b) have led to recent revival of the suggestion that the dark matter is self-interacting (Spergel & Steinhardt, 2000).

Simulations of halo formation in a CDM universe predict that dark matter halos have a cuspy inner profile. Most simulations find that, at small radii, $`\rho r^\alpha `$, with $`\alpha 1`$ to $`1.5`$ (Navarro et al., 1996b; Ghigna et al., 2000; Jing & Suto, 2000; Klypin et al., 2000). An inner density profile this steep may be consistent with the dark matter density profile inferred from many high surface brightness galaxies and clusters of galaxies (Eke et al. (2000), but see Borriello & Salucci (2001) ). However, low surface brightness galaxies, and dwarf irregular galaxies in particular, require an inner density profile which is less steep (Moore (1994); Burkert (1995); Flores & Primack (1994); for a contradictory view see van den Bosch & Swaters (2000)). In addition, Debattista & Sellwood (2000) have emphasized that too much dark matter in the centers of galaxies is incompatible with rotating bars, which are common in massive galaxies such as the Milky Way.

A separate problem for standard CDM model is that it predicts many more dark matter satellites of a Milky Way-sized galaxy than are observed (Moore et al., 1999; Klypin et al., 1999). This substructure cannot simply be dark satellites in which stars have not formed, because they would thicken the disk of the Milky Way as they passed through the plane of the disk.

While it is possible that the inclusion of baryonic physics would lead to a less cuspy dark matter halo, perhaps due to supernova driven gas outflow in dwarf galaxies (Navarro et al., 1996a), this seems unlikely, as cooling baryons would lead to further concentration of the dark matter halo through adiabatic compression. Exploding gaseous disks will drive a strong wind off the plane, but most simulations show that the cold, dense component of a multi-phase medium would be difficult to remove (Mac Low, 2001). Furthermore, these solutions to the cusp problem does not address the problem of too much substructure on small scales.

Recent attention has focused instead on the possibility that the dark matter is self-interacting. A remarkable number of different models have been proposed. For example, several authors have considered a model in which the dark matter is a self-interacting scalar field with quartic potential (Peebles, 2000; Peebles & Vilenkin, 1999; Goodman, 2000; Riotto & Tkachev, 2000). Hu et al. (2000) have suggested a model in which cores are produced in dwarf galaxies due to the wave nature of dark matter, if the matter is an ultra-light scalar particle. Another interacting dark matter candidate is Q-balls, non-topological solitons with properties that allow it to scatter, possibly with a time-dependent cross-section, or to merge (Kusenko & Steinhardt, 2001). Several authors have simulated halo evolution and formation in scenarios in which the dark matter scatters. All authors agree that on short time scales, scattering reduces the logarithmic slope of the inner density profile substantially; some authors have found that the inner density profile subsequently steepens as the halo undergoes core collapse (Kochanek & White, 2000; Yoshida et al., 2000), while others disagree (Burkert, 2000; Davé et al., 2001).

Self-annihilating dark matter is another possible solution to the problem (Kaplinghat et al., 2000), although the particle physics must be carefully contrived so that the dark matter does not completely self-annihilate in the early universe. An analytic calculation of the density of a self-annihilating dark matter halo demonstrates the formation of a core. Physically, the core develops because annihilation becomes very efficient at removing mass from the halo above a density that depends on the annihilation cross section. Kaplinghat et al. (2000) point out that this annihilation also causes an adiabatic expansion of the halo, and numerically calculate the effects of adiabatic expansion; the core density including expansion is smaller than that predicted from annihilation alone by a factor of ten (for dwarf galaxies) to three (for clusters). In addition, Hui (2001) points out that unitarity bounds require that annihilation be accompanied by scattering.

In this paper, we simulate the behavior of an isolated dark matter halo which undergoes both scattering and annihilation. In section 2, we describe our dark matter model and our simulation methods. In section 3 we present the results of our simulations, and in section 4 we compare our results with those of previous authors.

## 2 Model

### 2.1 Initial halo and units

We simulated the evolution of an isolated halo; the initial density profile is the Hernquist (1990) profile,

$$\rho =\frac{M}{2\pi r(r+r_s)^3},$$
(1)
where $`M`$ is the total mass of the halo, and $`r_s`$ is the scale radius. We take the dynamical time for this model to be the circular orbit time at the scale radius $`r_s`$,

$$t_{dyn}=4\pi \sqrt{\frac{r_s^3}{GM}}$$
(2)
In the simulations below, we work in units in which $`G=1`$, $`M=1`$, and $`r_s=1`$.

The profile has been generated to guarantee that the center of mass of the halo is also the origin used in calculating the density of the Hernquist profile. We first choose the energy for each particle, then select a radius for the particle; these determine the magnitude of the velocity of the particle. We adopt an isotropic velocity distribution. For each particle generated with position $`𝐫`$, a second particle is generated with position $`𝐫`$; the second particle has the same energy as the first. The distribution is truncated at a radius of 100$`r_s`$. The direction of the velocity of each particle is generated independently, but the net velocity of the halo is set to zero.

### 2.2 Model for scattering and annihilation.

We adopt a model similar to that used by Burkert (2000). The probability for interaction for a particle which moves a distance $`\mathrm{\Delta }x`$ in time step $`\mathrm{\Delta }t`$ is

$$P=\frac{\mathrm{\Delta }x}{\lambda },$$
(3)
where $`\lambda =1/(n\sigma )`$ is the mean free path for interactions; this expression is valid only when $`\mathrm{\Delta }x\lambda `$. Our simulations include both scattering and annihilation. We use the notation $`\sigma _a^{}`$ to denote the annihilation cross section per unit mass and $`\sigma _s^{}`$ to denote the scattering cross section per unit mass.

The probability that a *pair* of particles scatters in time step $`\mathrm{\Delta }t`$ is

$$P_{s,ij}=\sigma _s^{}v_{ij}\rho _l\mathrm{\Delta }t,$$
(4)
where $`\rho _l`$ is the local density, and $`v_{ij}`$ is the *relative* velocity of particles $`i`$ and $`j`$. Scattering is completely elastic, and is isotropic in the center of momentum frame of the pair of particles.

As pointed out by Kaplinghat et al. (2000), only $`s`$-wave annihilation is consistent with observations. For $`s`$-wave annihilation, $`\sigma |v|`$ is independent of velocity, so that the cross section for annihilation per unit mass, $`\sigma _a^{}`$, is given by

$$\sigma _a^{}=\frac{\stackrel{~}{\sigma }_a}{m}\frac{v_a}{v},$$
(5)
where $`\stackrel{~}{\sigma }_av_a`$ is determined by the scattering potential. The numerical value we assume for the scattering cross section is given below.

The probability for an annihilation interaction in time step $`\mathrm{\Delta }t`$ is:

$$P_a=\sigma _a^{}\rho _lv\mathrm{\Delta }t=\frac{\stackrel{~}{\sigma }_av_a}{m}\rho _l\mathrm{\Delta }t,$$
(6)
where $`\rho _l`$ is the local density. The annihilation probability is independent of velocity.

We model the decay of the dark matter particles represented by our simulation particles by assigning to each simulation particle a probability for annihilation given by equation (6). If the local density is large enough, $`P_a`$ is the probability that some of the dark matter particles represented by the simulation particle would interact and decay. If a simulation particle decays, its mass $`m_i`$ is reduced by fraction $`F`$,

$$m_{i,new}=(1F)m_{i,old}.$$
(7)
Thus, the effective probability for annihilation of individual dark matter particles in our model is the product $`FP_a`$. The mass lost in an annihilation is a pre-set fraction of the current mass of the particle undergoing annihilation. In other words, the *fractional* mass lost in an annihilation is held fixed. The absolute mass loss, $`m_{i,old}m_{i,new}`$, decreases as a particle undergoing repeated annihilation loses mass.

### 2.3 Numerical Parameters

Our simulations were run on a GRAPE3-AF special purpose hardware card (Okumura et al., 1993) using a direct summation code. The time step was chosen to be small enough to satisfy both the Courant condition and to ensure that $`v\mathrm{\Delta }t<\lambda `$. Each run contained 30,000 particles. Detailed numerical parameters for each run are shown in Table 1. Runs S0A0, S3A0, and S3A1 ran for roughly 10 dynamical times with approximately 2100 steps per dynamical time; a fourth run, S.3A.1, was done with a much larger Plummer softening and time step; it ran for 100 dynamical times with approximately 400 steps per dynamical time.

One run, S0A0, had no interactions, and was used to determine the stability of the integration. The run S3A0 had no annihilations, but non-zero scattering cross-section. Two runs with annihilation and scattering were performed. One, S3A1, was short (the final time was a tenth of a Hubble time assuming a 1010M halo), but had the same scattering cross-section as the scattering-only run. The other, S.3A.1, ran for ten times longer, with both annihilation and scattering cross section reduced by a factor of ten. In both annihilation runs the scattering and annihilation *probabilities* are equal at the scale radius; inside of that radius scattering will be enhanced relative to annihilation because the scattering probability increases with velocity. The fractional mass loss per annihilation in both was $`F=5\%`$.

The cross section per unit mass for scattering can be calculated in physical units from

$$\sigma _{phys,s}^{}=.48\sigma _s^{}\left(\frac{10^{10}M_{}}{M}\right)\left(\frac{R_s}{1\mathrm{k}\mathrm{p}\mathrm{c}}\right)^2\mathrm{cm}^2\mathrm{g}^1.$$
(8)
The cross-section for annihilation is given by

$$\sigma _{phys,a}^{}=.98\sigma _a^{}\left(\frac{10^{10}M_{}}{M}\right)\left(\frac{R_s}{1\mathrm{k}\mathrm{p}\mathrm{c}}\right)^2\left(\frac{100\mathrm{k}\mathrm{m}\mathrm{s}^1}{v}\right)\mathrm{cm}^2\mathrm{g}^1$$
(9)
For $`R_s=2`$ kpc and $`M=10^{10}`$M, the half mass dynamical time is .168 Gyr. Our cross-sections are $`\sigma _s^{}=5.76\mathrm{cm}^2\mathrm{g}^1`$ and $`\sigma _a^{}=5.88\mathrm{cm}^2\mathrm{g}^1(100\mathrm{k}\mathrm{m}\mathrm{s}^1/v)`$ for the run S3A1; in the run S.3A.1, both cross sections are a factor of ten smaller. Our values of $`\sigma _s^{}`$ span the range suggested by Wandelt et al. (2000). Kochanek & White (2000) considered a broader range of values for $`\sigma _s^{}`$; our choice for our runs with the largest scattering cross section is in the middle of their range.

In both annihilating runs, we chose our annihilation cross section to be comparable to the scattering cross section. However, this annihilation cross section in the run S3A1 is about 300 times larger than that suggested by Kaplinghat et al. (2000) based on fits of rotation curves to dwarf galaxies. The run S.3A.1, with an annihilation cross section ten times smaller, is closer to what they suggest. With a mass loss factor $`F`$ of 5%, the effective cross section $`F\sigma _a^{}`$ is comparable to the value required to match observations. In both cases, the scattering cross section is above the minimum required by unitarity bounds (Hui, 2001).

## 3 Results

The inclusion of annihilation in our simulations causes a reduction in the central density of the halo. Figure 1 shows the evolution in the density profile of a halo with both annihilation and scattering; the top panel shows a short run with high cross section (S3A1), and the bottom a run with smaller cross section run for longer time (S.3A.1). The central density rapidly decreases, and remains low for the duration of the simulations, about 100 (10) dynamical times for run S.3A.1 (S3A1).

The core is not due to numerical effects. The force softening is .05 (.01) in the run S.3A.1 (S3A1); the flattening of the inner density profile we observe occurs at radii 2–10 (10–50) times larger than that and is unlikely to have been caused by softening. To test the stability of our code and the effects of our relatively limited particle number, we ran a simulation with no interactions for ten dynamical times. The simulation with no interactions S0A0 shows no evidence of the development of a core. The density profile of the halo was essentially unchanged over the course of the simulation. Figure 1 shows the beginning profile and the final profile for this run.

We have examined the mechanism by which the core forms in our interaction simulations. There are three related mechanisms one might imagine contributing to the development of a core. First, mass lost due to annihilation and the disappearance of the relativistic decay products from the halo directly decreases the central density of the halo. Second, this direct mass loss reduces the gravitational potential of the halo, causing particles in the central region to become less bound. Third, energy is transferred from the central region of the halo outward by scattering, depressing the central density.

Although scattering at first suppresses the central density, core collapse begins fairly rapidly in halos with scattering only, as reported by Kochanek & White (2000). We have run one simulation with scattering alone to determine whether the core that initially forms persists. The results of this run are shown in Figure 2. The central density initially drops and then increases after a few dynamical times.111This run was terminated after 5.8 dynamical times. At that point, the central density became so large that the typical particle was moving more than one mean free path in a time step. In the absence of annihilation, scattering does not produce a persistent core in simulations of *isolated* halos. It is unclear whether or not cores persist in high resolution simulations of halo formation in a cosmological setting (Yoshida et al., 2000; Davé et al., 2001).

In order to differentiate between direct mass loss from annihilation and a change in the density profile due to a change in gravitational potential, we calculate the change in number density of simulation particles as a function of radius. If direct mass loss due to annihilation is causing most of the decline seen in the density profile, the number density profile ought to be unchanged, with particles in the central region of the halo simply becoming less massive. Instead, as can be seen in Figure 3, which shows the evolution in the number density of particles early in the simulation, we find that the number density rapidly diminishes in the center of the halo, indicating that much of the change in central density is due to the motion of particles away from the center, not from direct mass loss. The trend of particle number density decreasing in the center continues throughout the simulation, as can be seen from Figure 1. Density profiles are calculated with a fixed number of particles per bin, so the location of the first point in the density profile is a radius that encloses a fixed number of particles; this point moves farther from the origin throughout the simulation. Furthermore, the total mass lost by the halo over the entire simulation is only 12% of the initial halo mass.

The halo in the annihilating model forms a core very rapidly, and the size of the core continues to grow over the course of the simulation. We fit two different functional forms to the density profile of the halo. One is a modification of the initial Hernquist profile 1 to allow for a core,

$$\rho _C=\frac{\rho _0x_c}{(x+x_c)(1+x)^3},$$
(10)
where $`x=r/r_s`$, and $`x_c=r_c/r_s`$ is the core radius in units of the scale radius. The other density profile was derived by Kaplinghat et al. (2000); their profile takes into account modifications of the initial halo density profile due to direct annihilation only, without scattering. Adiabatic expansion of the core as mass is lost decreases the core density. We have derived the profile appropriate for an initial Hernquist model,

$$\rho =\frac{\rho _0x_c}{x(1+x)^3+x_c}$$
(11)
(note that $`x_c`$ is related to the annihilation cross section in Kaplinghat’s model (Eq. 13); in doing our fits we have treated $`x_c`$ as a free parameter to see how well the shape of the density profile is fit by the function). Fits to both functions are shown at four different epochs in Figure 4. The two profiles differ only in the innermost region of the halo. Burkert (1995) found that dwarf galaxies are well fit by the profile

$$\rho _B=\frac{\rho _0r_c^3}{(r+r_c)(r^2+r_c^2)},$$
(12)
where $`\rho _0`$ is the central density of the halo, and $`r_c`$ its core radius. We do not expect the Burkert profile to be a good fit to our halo because the initial halo profile has both a different radial dependence at large radii, and a length scale (the scale radius) in the initial conditions separate from the core radius.

The halo shows a clear evolution from a Kaplinghat profile to a modified Hernquist profile. At early times, the Kaplinghat model is a better fit than the modified Hernquist model. At later times, the modified Hernquist model is a better fit. This is good news, because the modified Hernquist profile is closer in form to the Burkert profile. The evolution in profile shape makes sense; early in the simulation the adiabatic effect of annihilation on the gravitational potential will be small, and the Kaplinghat profile should be a good fit. The core radius in the modified Hernquist profile continues to grow throughout the simulation, as shown in Figure 5. The final core radius is 3 kpc, in the range suggested by Burkert (1995) for a halo of the mass we studied.

## 4 Discussion

Our simulations demonstrate that annihilation produce cores of radius a few kiloparsecs in dark matter halos of dwarf galaxies; we now discuss whether this dark matter model preserves the successes of non-interacting CDM, the extent to which it addresses the problem of too much substructure in galactic dark matter halos, and the plausibility of the particle physics of the model.

There is little evidence for cores in clusters; for an annihilation model to be plausible it must modify small scale behavior of dark matter while leaving large scale behavior unchanged. In the analytic model of Kaplinghat et al. (2000), the core radius grows linearly in time, and is given by

$$x_c=\frac{r_c}{r_s}=\rho _s\sigma _at/m=\rho _sF\sigma _a^{}t$$
(13)
where $`\rho _s`$ is the characteristic density of the NFW profile and $`t`$ is the age of the halo. Since lower mass halos have higher density, and formed earlier, than high mass halos, one would expect a decreasing core radius with increasing mass. Figure 6 shows the trend of core radius with mass, assuming the core radius is given by equation (13); the core radius has been normalized to that of a $`10^{10}h^1M_{}`$ halo (approximately 3 kpc) to eliminate the dependence of the core radius on cross section. Characteristic densities and halo ages were calculated using the Eke et al. (2000) model, assuming a $`\mathrm{\Omega }_0=.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=1\mathrm{\Omega }_0=.7`$ flat cosmology with $`h=.65`$. Most of the decrease in core radius is due to the decrease in characteristic density as mass increases; the age difference of the lowest and highest mass halo accounts for only 20% of the difference in core radius. Our model clearly creates large cores in low mass/high density halos without producing them in high mass/low density halos. Thus if the annihilations lead to a core radius of 3 kpc for dwarf galaxies, as appropriate, the predicted dark matter core radius in rich clusters would be $`120`$ pc, which is deep within the central galaxy.

The core structure of the dark matter halo of a high surface brightness galaxies is unclear. Eke et al. (2000) argue that the concentrations of dark matter halos predicted in the standard $`\mathrm{\Lambda }`$CDM model are consistent with both the Tully-Fisher relationship and the amount of dark matter interior to the solar radius. On the other hand, Weiner et al. (2001a) argues that detailed mass modeling of the barred spiral NGC 4123 requires either a modest core in the dark matter halo, with a core radius of 1–6 kpc, or an NFW profile with very low concentration. Arguing more generally, Debattista & Sellwood (2000, 1998) find that a high density dark matter halo in a barred galaxy leads to rapid slowing of bar rotation. For the same choice of annihilation cross-section that generates a 3 kpc core in dwarf galaxies, the core radius in a $`10^{12}`$ M galaxy such as the Milky Way would be $`1.2`$ kpc, and would resolve the difficulties of bar spirals in standard $`\mathrm{\Lambda }`$CDM discussed by Debattista & Sellwood (2000).

The standard $`\mathrm{\Lambda }`$CDM model also predicts more substructure on galactic scales than is observed (Moore et al., 1999; Klypin et al., 1999). Although the simulations we have done cannot directly address the question of substructure, Figure 6 can be extended to low mass halos. As halo mass shrinks, the core radius increases while the virial radius increases. The two are equal for halos of mass 3$`\times 10^7`$M; halos with mass $`10^8`$ solar masses have a core radius that is half the virial radius. This suggests that halos smaller than a few times $`10^8`$ solar masses would not be present in an annihilating dark matter model. Further simulation of structure formation in a cosmological context with annihilating dark matter would be helpful.

One of the problems with earlier self-interacting dark matter proposals (Carlson et al., 1992; Machacek, 1994) is that ram pressure stripping will destroy small halos too efficiently (De Laix et al., 1995). Newer models in which the dark matter behaves like a viscosity-free fluid (Peebles, 2000; Peebles & Vilenkin, 1999; Goodman, 2000; Riotto & Tkachev, 2000) avoid this problem, but may overly-suppress the growth of structure on small scales in the linear regime (Peebles, 2000), a problem shared by the earlier self-interacting dark matter models (Machacek, 1994).

The model proposed by Carlson et al. (1992) has some similarities to our model. They investigated the properties of dark matter that self-annihilated (that is, interactions occurred in which the annihilation products were the dark matter particles). Their model affects the structure of halos by heating dense regions, and affects fluctuations in the early universe through the introduction of a non-trivial Jean’s length for the dark matter. Our model differs in that we assume the dark matter annihilates to something relativistic that no longer contributes to the gravitational potential of a halo after annihilation. This should lead to a greater effect on the core properties of the halo, through adiabatic changes to the gravitational potential, than the model of Carlson et al. (1992). The interacting dark matter model proposed by Spergel & Steinhardt (2000) is constructed so as to avoid excessive suppression of the power spectrum on small scales. The dark matter is presumed to be cold, and the interaction cross sections are small enough that the probability of interaction is small at the time of recombination.

There are additional problems with an annihilating dark matter model. First, the amount of energy released in dark matter annihilation is tremendous; in our simulation of a single halo, the mass loss from annihilation (in whatever the decay products are) has an energy equivalent of $`4.7\times 10^{45}`$ erg s-1 averaged over the lifetime of the universe. This is far too large to be consistent with observations of the gamma ray background. If the dark matter annihilates, the decay products cannot include photons. Second, as discussed at length by Kaplinghat et al. (2000), the particle physics underlying an annihilating dark matter model must be carefully contrived so as to avoid complete annihilation of the dark matter in the early universe. The suppression of annihilation must be in place for $`z10000`$, which is a rather low energy for any reasonable phase transitions in the dark matter sector.

## 5 Conclusion

Annihilating dark matter leads to halo density profiles consistent with those observed in dwarf spiral galaxies. Annihilation leads to a rapid reduction in the central density of a dark matter halo, driven mainly by adiabatic expansion of the core region as matter is lost. An annihilation cross section sufficient to generate a core in a dwarf galaxy halo generates a moderate core in a Milky Way sized halo, and a negligible core in a rich cluster.

The core radius decreases rapidly with halo mass, primarily because larger objects form later in hierarchical clustering models, when the mean density of the universe is lower. Objects smaller than $`3\times 10^8`$M are unlikely to be present in an annihilating dark matter model, because the core radius will be large compared to the virial radius. A reduction in structure on this scale would be useful; CDM predicts many more halos of this size than are observed.

On the other hand, there is no compelling motivation for such a complex dark matter candidate. A special mechanism is needed to suppress annihilation of the dark matter in the early universe, and the dark matter must decay to particles that do not interact with ordinary matter. Although a particle physics model with the required properties could no doubt be constructed, such a model would be contrived.

Further study of the impact of annihilation in a cosmological setting would be useful in evaluating the promise of this model.

## 6 Acknowledgements

This work was supported in part by NSF grant AST00-71048.