text:
# Search for Decays of 𝐵⁰ Mesons Into Pairs of Leptons: 𝐵⁰→𝑒⁺𝑒⁻,𝐵⁰→𝜇⁺𝜇⁻ and 𝐵⁰→𝑒^±𝜇^∓.

## I Introduction

In the Standard Model the decays $`B^0e^+e^{}`$ and $`B^0\mu ^+\mu ^{}`$ are allowed via the diagrams shown in Fig 1. The Standard Model predictions for the branching fractions are $`(B^0e^+e^{})=1.9\times 10^{15}`$ and $`(B^0\mu ^+\mu ^{})=8.0\times 10^{11}`$, respectively. These rates put them beyond the reach of current experiments. Their observation at CLEO would provide clear evidence for physics beyond the Standard Model.

In two Higgs doublet models, the corresponding branching fractions can be strongly enhanced due to additional diagrams involving Higgs bosons . Supersymmetric particles can further modify the expectation, but calculations of SUSY contributions have not yet been made. In some models with $`Z`$-mediated flavor changing neutral currents  the corresponding branching fraction can be enhanced by a factor of 400 compared to Standard Model predictions.

The decay $`B^0e^\pm \mu ^{}`$ is forbidden in the Standard Model by lepton number conservation. However it can occur in Pati-Salam leptoquark models, due to the existence of particles that couple to both quarks and leptons.

Currently the best 90% confidence level limits on the branching fractions are $`(B^0e^+e^{})<59\times 10^7`$ (from CLEO), $`(B^0e^\pm \mu ^{})<35\times 10^7`$ (from CDF) and $`(B^0\mu ^+\mu ^{})<6.8\times 10^7`$ (from CDF).

## II Data sample

We use data taken with the CLEO II/II.V detector operating at the Cornell Electron Storage Ring (CESR). The sample consists of 9.1 $`fb^1`$ taken on the $`\mathrm{{\rm Y}}(4S)`$ resonance, corresponding to approximately $`9.7\times 10^6B\overline{B}`$ pairs. An additional sample of 4.5 $`fb^1`$ accumulated approximately 60 MeV below the $`\mathrm{{\rm Y}}(4S)`$ resonance is used to estimate the non-resonant background. Charge conjugation is implied throughout this paper and equal production of charged and neutral $`B`$ mesons at the $`\mathrm{{\rm Y}}(4S)`$ resonance is assumed. The CLEO II charged particle tracking system consists of a 6-layer straw tube chamber, a 10-layer precision tracking chamber and a 51 layer main drift chamber. In the CLEO II.V configuration the straw tube chamber was replaced by a 3-layer, double-sided silicon tracker. Beyond the tracking chambers, but inside the 1.5T solenoid magnet, are a time-of-flight system and an electromagnetic calorimeter consisting of 7800 CsI crystals. Electron identification is performed with a likelihood function combining calorimeter ($`E/p`$, the ratio of energy deposited by particle in calorimeter to its momentum) and specific ionization ($`dE/dx`$) information. Muons are identified using proportional counters placed at various depths in the steel return yoke of the magnet.

Since, at a symmetric $`e^+e^{}`$ collider, $`B`$ mesons at the $`\mathrm{{\rm Y}}(4S)`$ resonance are produced nearly at rest, the decays $`B^0e^+e^{}`$, $`B^0\mu ^+\mu ^{}`$ and $`B^0e^\pm \mu ^{}`$ will contain two nearly back-to-back leptons with momenta $`p2.6\mathrm{GeV}`$. In order to select $`B\overline{B}`$ events, we require at least five charged tracks in the event with two well reconstructed and oppositely charged tracks identified as leptons. The lepton pair is required to have a mass within three standard deviations from the the mass of the $`B^0`$ meson. We use the beam constrained mass, $`M(B)=\sqrt{E_{beam}^2(\stackrel{}{p}_{l^+}+\stackrel{}{p}_l^{})^2}`$. The resolution on $`M(B)`$ is 2.8 MeV for the $`B^0e^+e^{}`$ decay and 2.6 MeV for the $`B^0\mu ^+\mu ^{}`$ and $`B^0e^\pm \mu ^{}`$ decays. For signal events, $`M(B)`$ is required to be within 3 standard deviations of the $`B^0`$ meson mass.

We calculate the energy difference between the candidate $`B^0`$ meson and the beam energy $`\mathrm{\Delta }E=E_{l^+}+E_l^{}E_{beam}`$ and require that the difference should be $`|\mathrm{\Delta }E|<75\text{ MeV}`$ for $`B^0\mu ^+\mu ^{}`$ decay, and for $`B^0e^+e^{}`$ and $`B^0e^\pm \mu ^{}`$ decays it should satisfy $`100\text{ MeV}<\mathrm{\Delta }E<75\text{ MeV}`$. The asymmetric $`\mathrm{\Delta }E`$ requirement in $`B^0e^+e^{}`$ and $`B^0e^\pm \mu ^{}`$ channels takes into account energy loss due to final state radiation from the electrons.

The main background comes from the processes $`e^+e^{}q\overline{q}(q=u,c,s,d)`$ and $`e^+e^{}\tau ^+\tau ^{}`$. Such events usually exhibit a two jet structure and produce high momentum, approximately back-to-back tracks which satisfy the requirements imposed on our candidate events. We calculate the angle $`\theta _T`$ between the thrust axis of the $`B^0`$ candidate tracks and the thrust axis of the remaining charged and neutral tracks in the event. The distribution of $`|\mathrm{cos}\theta _T|`$ is flat for signal events and strongly peaked at 1 for background events. We require $`|\mathrm{cos}\theta _T|<0.9`$.

## III Signal simulation

In order to take into account final state radiation (FSR) from the final state leptons, the matrix element for the decay $`J/\psi e^+e^{}(\gamma )`$ from Ref. , modified in order to simulate $`B^0l^+l^{}(\gamma )`$ with $`l^+l^{}=e^+e^{},\mu ^+\mu ^{}`$ is included in the event generator. The detector simulation program is based on GEANT , the simulated events are processed as the data. The overall efficiencies for detecting the decays $`B^0e^+e^{}`$ and $`B^0\mu ^+\mu ^{}`$ are 37.3% and 46.2%, respectively, without FSR. Including FSR reduces these efficiencies to 31.1% and 42.4%, respectively. We use the efficiencies obtained with FSR included for the derivation of limits on branching fractions.

The efficiency for the decay $`B^0e^\pm \mu ^{}`$ is 43.6% without FSR. There is no matrix element available for the decay $`B^0e^\pm \mu ^{}(\gamma )`$. In this mode we use Monte Carlo without FSR. Since final state radiation reduces the efficiencies for $`B^0e^+e^{}`$ and $`B^0\mu ^+\mu ^{}`$ by 20% and 9%, respectively, we estimate the effect of FSR in $`B^0e^\pm \mu ^{}`$ to be 15% and include this in the derivation of the upper limit on the branching fraction as a systematic uncertainty.

## IV Results

The background from the process $`e^+e^{}q\overline{q}`$ results from misidentification of hadrons as leptons and has been determined as follows. First, the analysis is redone without imposing the lepton identification requirements to determine the total number of 2-track events passing all other selection criteria in the signal and sideband regions ($`5.260\text{ GeV}<M(B)<5.270\text{ GeV}`$ and $`|\mathrm{\Delta }E|<200\text{ MeV}`$) for both the on-resonance and off-resonance data samples. These are used to determine the expected number of events in the on-resonance signal region. (This is done by multiplying the number of events seen in the signal region in off-resonance data by the ratio of numbers of events seen in the sideband region in on- and off-resonance data.) The probabilities for these tracks to be identified as leptons is determined from the relative abundances of particles ($`\pi ,K`$, proton, electron, muon) and the fake rates. The particle abundances for the $`B^0`$ candidate tracks in $`q\overline{q}`$ events have been determined using Monte Carlo and were found to be $`0.58`$ for $`\pi ^\pm `$, $`0.25`$ for $`K^\pm `$, $`0.15`$ for $`p,\overline{p}`$, and $`0.004`$ for electrons and muons. For $`\pi ^\pm `$ and $`K^\pm `$ the fake rates are determined experimentally from samples of $`D^+D^0\pi ^+K^{}\pi ^+\pi ^+`$ decays. Proton fake rates are taken from Monte Carlo. The expected number of backgrounds events from $`e^+e^{}q\overline{q}`$ are $`0.01\pm 0.01`$, $`0.14\pm 0.03`$ and $`0.09\pm 0.03`$ for the $`ee,\mu \mu `$ and $`e\mu `$ channels, respectively.

Backgrounds from $`e^+e^{}\tau ^+\tau ^{}`$ events are estimated using a Monte Carlo sample analyzed without lepton identification, using the procedure described above to determine the expected rate of events with the two candidate tracks identified as leptons. The particle abundances for the $`B^0`$ candidate tracks in simulated $`\tau ^+\tau ^{}`$ events were found to be $`0.74`$ for $`\pi ^\pm `$, $`0.04`$ for $`K^\pm `$, and $`0.11`$ for electrons and muons. The expected numbers of background events are $`0.10\pm 0.07`$ , $`0.08\pm 0.06`$ and $`0.40\pm 0.20`$ events for the $`ee`$, $`\mu \mu `$ and $`e\mu `$ channels, respectively. The background from $`B\overline{B}`$ events is estimated from Monte Carlo and is found to be less than 0.01 events in all modes.

We find no events in the signal regions for the $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ modes and two events in the signal region for the $`e^\pm \mu ^{}`$ mode. These events are consistent with the expected background. To conservatively account for the uncertainty in the efficiency determination, we reduce the efficiency by the estimated systematic error of 7.6% in the $`e^+e^{}`$ and $`\mu ^+\mu ^{}`$ modes and 16.3% in the $`e^\pm \mu ^{}`$ mode when calculating the upper limits. The latter error is dominated by the uncertainty from final state radiation. A summary of the efficiencies, event yields and limits is given in Table I.

## V Limits on Pati-Salam leptoquarks

One of the simplest models that assumes symmetry between quarks and leptons is the Pati-Salam model. This model predicts heavy spin-one gauge bosons, called leptoquarks (LQ), that carry both color and lepton quantum numbers. While most leptoquark models assume that the leptoquarks couple within one generation only, Pati-Salam model allows cross-generation couplings. Consequently, Pati-Salam leptoquarks can mediate the decay $`B^0e^\pm \mu ^{}`$ . The relationship between $`(B^0e^\pm \mu ^{})`$ and leptoquark mass is

$$\mathrm{\Gamma }(B^0e^\pm \mu ^{})=\pi \alpha _s^2(M_{LQ})\frac{1}{M_{LQ}^4}F_{B^0}^2m_{B^0}^3R^2$$

where

$$R=\frac{m_{B^0}}{m_b}\left(\frac{\alpha _s(M_{LQ})}{\alpha _s(m_t)}\right)^{\frac{4}{7}}\left(\frac{\alpha _s(m_t)}{\alpha _s(m_b)}\right)^{\frac{12}{23}}.$$

The first attempt to constrain the Pati-Salam leptoquark masses by measuring the branching fraction for $`B^0e^\pm \mu ^{}`$ was made by the CDF collaboration. The result obtained was $`M_{LQ}>21.7\text{ TeV}`$. We follow the CDF procedure by using $`F_B=175\pm 30`$ MeV for $`B^0`$ decay constant , $`m_{B^0}=5279.1\pm 0.7\text{(stat.)}\pm 0.3\text{(syst.)}`$ MeV, $`m_t=176.0\pm 6.5`$ GeV  and obtain a new limit on the Pati-Salam leptoquark mass of $`M_{LQ}>27\text{ TeV}`$ at 90% confidence level.

In conclusion we find no evidence for the decays $`B^0e^+e^{},\mu ^+\mu ^{}`$ and $`e^\pm \mu ^{}`$ and obtain upper limits on the corresponding branching fractions of $`(B^0e^+e^{})<8.3\times 10^7`$, $`(B^0\mu ^+\mu ^{})<6.1\times 10^7`$, $`(B^0e^\pm \mu ^{})<15\times 10^7`$ at 90% confidence level. The limit on $`(B^0e^\pm \mu ^{})`$ corresponds to a lower limit of 27 TeV on the Pati-Salam leptoquark mass at 90% CL.

We gratefully acknowledge the effort of the CESR staff in providing us with excellent luminosity and running conditions. I.P.J. Shipsey thanks the NYI program of the NSF, M. Selen thanks the PFF program of the NSF, M. Selen and H. Yamamoto thank the OJI program of DOE, M. Selen and V. Sharma thank the A.P. Sloan Foundation, M. Selen and V. Sharma thank the Research Corporation, F. Blanc thanks the Swiss National Science Foundation, and H. Schwarthoff and E. von Toerne thank the Alexander von Humboldt Stiftung for support. This work was supported by the National Science Foundation, the U.S. Department of Energy, and the Natural Sciences and Engineering Research Council of Canada.


format: text

text:
# Dynamics of evaporative cooling in magnetically trapped atomic hydrogen

## Abstract

We study the evaporative cooling of magnetically trapped atomic hydrogen on the basis of the kinetic theory of a Bose gas. The dynamics of trapped atoms is described by the coupled differential equations, considering both the evaporation and dipolar spin relaxation processes. The numerical time-evolution calculations quantitatively agree with the recent experiment of Bose-Einstein condensation with atomic hydrogen. It is demonstrated that the balance between evaporative cooling and heating due to dipolar relaxation limits the number of condensates to $`9\times 10^8`$ and the corresponding condensate fraction to a small value of 4% as observed experimentally.

Pacs. 03.75.Fi, 05.20.Dd, 32.80.Pj 11footnotetext: Present address: Department of Physics, Faculty of Science, Kumamoto University, Kurokami 2-39-1, Kumamoto 860-8555, Japan.

The recent development of atom-manipulation techniques  has realized Bose-Einstein condensation (BEC) in magnetically trapped alkali-metal atoms  and atomic hydrogen . In all these experiments, evaporative cooling was adopted at the final stage of the cooling procedures and was essential in obtaining the extremely low temperatures needed for the quantum degeneracy. The cooling mechanism of this powerful method is based on both the selective removal of energetic atoms through evaporation and collisional rethermalizations among the remaining atoms . Evaporative cooling itself is stimulating theoretical studies of the process of condensate formation in nonequilibrium atomic gases .

In BEC experiments with atomic hydrogen, evaporative cooling has been implemented by lowering the potential height of the saddle point at one end of a magnetic trap with cylindrical symmetry . This approach suffers from a reduction of the dimension of evaporation at low temperatures due to low elastic collision rate  since the $`s`$-wave scattering length, $`a`$, of a hydrogen atom is anomalously small, about two orders of magnitude smaller than that in alkali-metal atoms. Heating caused by the dipolar spin relaxation then easily retards the evaporative cooling and prevents further cooling before reaching the critical temperature of BEC. To overcome this problem, Fried et al. adopted the evaporative cooling induced by radio-frequency (rf) magnetic field just after the “saddle-point” evaporative cooling . The atom ejection technique utilizing the transition between trapped and untrapped spin states enables efficient three-dimensional evaporation even at low temperatures and in a highly anisotropic magnetic potential . BEC was finally achieved as a result of the heating-cooling balance in this efficient “rf-induced” evaporative cooling.

There are several characteristic features of BEC in atomic hydrogen as compared with that in alkali-metal atoms . The small mass, $`m`$, of hydrogen resulted in the transition temperature of $`50\mu `$K, the highest among BEC experiments. While the condensate fraction, i.e., the ratio of condensates to the whole trapped atoms, remained a small value of about $`5`$%, a huge number of condensates containing about $`10^9`$ atoms were observed and the peak density of the condensates was about 25 times higher than that of the noncondensed atoms. Furthermore, with respect to the loss of condensates, both the three-body recombination process and the background gas collisions in a cryogenic environment of trap are negligible in the hydrogen system. Only dipolar spin relaxation becomes the dominant loss mechanism, and causes serious heating .

In this paper, the kinetic theory for evaporative cooling of a magnetically trapped Bose gas  is extended to include dipolar relaxation loss. This theory enables the consistent and quantitative investigations on the whole evaporative cooling process of the experiments, from the classical regime at the early stage of cooling to the quantum degenerate regime after the BEC transition. We applied it to the rf-induced evaporative cooling in the recent BEC experiment with atomic hydrogen  and quantitatively investigated the characteristic features of hydrogen BEC mentioned above. We expect further that our calculation studying the dynamics of trapped atoms in cooling process will give a good estimation for the growing of condensates in the experiment, since there is room for discussions on the number of condensates measured in the spectroscopic way .

Here, we briefly mention the formulation of the theory . During evaporative cooling, the applied rf-magnetic field effectively truncates the trapping potential, $`U(𝐫)`$, at the energy, $`ϵ_t`$, determined by the rf-field frequency. The thermalized distribution of noncondensed atoms in such a truncated potential is well approximated by the truncated Bose-Einstein distribution function

$$\stackrel{~}{f}(𝐫,p)=\frac{1}{\stackrel{~}{\xi }^1(𝐫)\mathrm{exp}(ϵ_p/k_BT)1}\mathrm{\Theta }(A(𝐫)ϵ_p),$$
(1)
where $`T`$ is the temperature, $`ϵ_p=p^2/2m`$ is the kinetic energy of atoms, and $`\mathrm{\Theta }(x)`$ is the step function. $`\stackrel{~}{\xi }(𝐫)`$ denotes the local fugacity including the mean field potential energy; $`\stackrel{~}{\xi }(𝐫)=\mathrm{exp}\{[\mu U(𝐫)2v\stackrel{~}{n}(𝐫)]/k_BT\}`$, where $`\mu `$ is the chemical potential, $`v=4\pi a\mathrm{}^2/m`$ is interaction strength between trapped atoms, and $`\stackrel{~}{n}(𝐫)`$ is the density profile of atoms in a truncated potential. The step function in Eq. (1) eliminates the momentum states whose kinetic energies exceed the effective potential height $`A(𝐫)=ϵ_tU(𝐫)2v\stackrel{~}{n}(𝐫)`$. The momentum integration of this truncated distribution function gives the density profile in a self-consistent way such that $`\stackrel{~}{n}(𝐫)=4\pi h^3\stackrel{~}{f}(𝐫,p)p^2dp`$, and similarly the internal energy density such that $`\stackrel{~}{e}(𝐫)=4\pi h^3ϵ_p\stackrel{~}{f}(𝐫,p)p^2dp+v\stackrel{~}{n}^2(𝐫)+U(𝐫)\stackrel{~}{n}(𝐫)`$. The total number of atoms, $`\stackrel{~}{N}`$, and the total internal energy, $`\stackrel{~}{E}`$, in a truncated potential are evaluated respectively by the spatial integrations of these density functions, $`\stackrel{~}{N}=\stackrel{~}{n}(𝐫)d𝐫`$ and $`\stackrel{~}{E}=\stackrel{~}{e}(𝐫)d𝐫`$. At low temperatures, the density profile in the condensed region is described as the sum of condensates, $`n_0(𝐫)`$, and saturated noncondensed atoms, $`\stackrel{~}{n}_n`$, such that $`\stackrel{~}{n}(𝐫)=n_0(𝐫)+\stackrel{~}{n}_n`$. Condensates obey the Thomas-Fermi distribution $`n_0(𝐫)=n_pU(𝐫)/v`$, where $`n_p`$ represents the peak density of condensates at the center of the potential. The density $`\stackrel{~}{n}_n`$ is evaluated through the truncated Bose-Einstein distribution function under the condition $`\stackrel{~}{\xi }(𝐫)=1`$.

The dynamics of trapped atoms during evaporative cooling is investigated on the basis of kinetic theory of a Bose gas . Hydrogen atoms are removed from the trapping potential through both evaporation and dipolar spin relaxation processes . The change rates (i.e., loss rates) of density functions, $`\dot{n}_{\mathrm{loss}}`$ and $`\dot{e}_{\mathrm{loss}}`$, are evaluated respectively as the sum of the contributions of both processes, such that

$`\dot{n}_{\mathrm{loss}}(𝐫)`$ $`=`$ $`\dot{n}_{\mathrm{ev}}(𝐫)G_2K(𝐫)\stackrel{~}{n}^2(𝐫),`$ (2)
$`\dot{e}_{\mathrm{loss}}(𝐫)`$ $`=`$ $`\dot{e}_{\mathrm{ev}}(𝐫)G_2K(𝐫)\stackrel{~}{e}(𝐫)\stackrel{~}{n}(𝐫),`$ (3)
where $`\dot{n}_{\mathrm{ev}}`$ and $`\dot{e}_{\mathrm{ev}}`$ are the evaporation rates of density functions, $`G_2`$ is the dipolar decay rate constant, and $`K`$ is the correlation function which describes the second-order coherence of trapped atoms . Both $`\dot{n}_{\mathrm{ev}}`$ and $`\dot{e}_{\mathrm{ev}}`$ were derived from a general collision integral of a Bose gas system in Ref. , and we adopt here the opposite sign of the notations of these rates defined in this reference. The bosonic feature of evaporation process is therefore included as the strong enhancement of these evaporation rates due to the stimulated scattering of atoms. Regarding the dipolar spin relaxation process, we assume the correlation function of an ideal Bose gas such that $`K(𝐫)=(n_0^2(𝐫)+4n_0(𝐫)\stackrel{~}{n}_n+2\stackrel{~}{n}_n^2)/2\stackrel{~}{n}^2(𝐫)`$ . The spatial integration of Eq. (3) gives important thermodynamic quantities such as the loss rate of total number of trapped atoms, $`\dot{N}_{\mathrm{loss}}`$, and that of the total internal energy, $`\dot{E}_{\mathrm{loss}}`$:

$`\dot{N}_{\mathrm{loss}}`$ $`=`$ $`{\displaystyle \dot{n}_{\mathrm{loss}}(𝐫)d𝐫}=\dot{N}_{\mathrm{ev}}\dot{N}_{\mathrm{dip}},`$ (4)
$`\dot{E}_{\mathrm{loss}}`$ $`=`$ $`{\displaystyle \dot{e}_{\mathrm{loss}}(𝐫)d𝐫}=\dot{E}_{\mathrm{ev}}\dot{E}_{\mathrm{dip}},`$ (5)
where $`\dot{N}_{\mathrm{ev}}(\dot{E}_{\mathrm{ev}})`$ and $`\dot{N}_{\mathrm{dip}}(\dot{E}_{\mathrm{dip}})`$ represent the evaporation rate and dipolar relaxation loss rate, respectively. These loss rates are complicated functions of temperature $`T`$, chemical potential $`\mu `$, and truncation energy $`ϵ_t`$. Finally, the system obeys the coupled differential equations

$$\frac{\stackrel{~}{N}}{t}=\dot{N}_{\mathrm{loss}},\frac{\stackrel{~}{E}}{t}=\dot{E}_{\mathrm{loss}}.$$
(6)
Time-evolution calculations were performed numerically on the assumption of quick rethermalization, which is appropriate for the slow evaporative cooling used in the usual BEC experiments . Thus, the system always stays in a quasi-thermal equilibrium state and is described by the truncated Bose-Einstein distribution function in Eq. (1).

Precise knowledge of experimental conditions is also necessary for our quantitative calculations. The magnetic trap adopted in the BEC experiment with atomic hydrogen  is modeled by the Ioffe-Pritchard potential $`U(\rho ,z)=\sqrt{(\alpha \rho )^2+(\beta z^2+\theta )^2}\theta `$ with the radial potential gradient $`\alpha =2.2\times 10^{23}`$ J/m, the axial curvature $`2\beta =6.8\times 10^{24}\mathrm{J}/\mathrm{m}^2`$, and the bias energy $`\theta =4.8\times 10^{28}`$ J (i.e., $`\theta /k_\mathrm{B}=35\mu `$K). At low energies, this potential is well approximated by the harmonic one with the radial oscillation frequencies $`\omega _\rho =\alpha /\sqrt{m\theta }=2\pi \times 3.9`$ kHz and the axial frequency $`\omega _z=\sqrt{2\beta /m}=2\pi \times 10`$ Hz. It should be noted here that, even after the BEC transition, many noncondensed hydrogen atoms still distribute in the high-energy region, where such harmonic approximation fails. We treated the Ioffe-Pritchard potential exactly for the quantitative description of noncondensed atoms. Accordingly, all density functions were calculated in the cylindrical coordinate.

As is well known, the efficiency of evaporative cooling strongly depends on the way the rf-field frequency is swept (i.e., forced evaporative cooling) . The high efficiency is obtained by the slow evaporative cooling which continues for the time on the order of 10 seconds. In the BEC experiment with atomic hydrogen, the rf-field frequency was swept according to the functional form  $`\nu (t)=\nu _\mathrm{i}\left({\displaystyle \frac{\nu _\mathrm{f}}{\nu _\mathrm{i}}}\right)^{\tau ^g}`$, where $`\nu _\mathrm{i}=23`$ MHz is the initial rf-frequency, $`\nu _\mathrm{f}=2`$ MHz is the final one, $`\tau =t/t_0`$ denotes the normalized time with the specific time $`t_0=25`$ s, and $`g=1.5`$ is the fixed parameter. The truncation energy $`ϵ_t`$ then changes according to the relation $`ϵ_t(t)=h\nu (t)\theta `$.

Furthermore, we employed the following experimental parameters : s-wave scattering length $`a=6.48\times 10^2`$ nm , the dipolar decay rate constant $`G_2=1.1\times 10^{15}`$ cm3/s , initial temperature $`T=120\mu `$K, initial peak density $`\stackrel{~}{n}(\mathrm{𝟎})=5.0\times 10^{13}`$ cm-3, and the initial number of trapped atoms in a truncated potential $`\stackrel{~}{N}=1.16\times 10^{11}`$. The time-evolution of evaporative cooling was calculated for 25 s, i.e., the normalized time $`\tau `$ changes from 0 to 1.

In Figs. 1(a) and 1(b), we show the time evolution of the total number of trapped atoms $`\stackrel{~}{N}`$ and that of temperature $`T`$, respectively. $`\stackrel{~}{N}`$ decreases monotonically from $`1.16\times 10^{11}`$ to $`8.7\times 10^8`$ during a 25-s evaporative cooling. The BEC transition occurs at 15.4 s ($`\tau =0.62`$) as indicated by the arrow in the figure. The decrease of temperature (from 120 to $`22\mu `$K) is not so large in comparison with that has been demonstrated in alkali-metal atoms . We can see that the cooling speed strongly slows down after the BEC transition in Fig. 1(b), indicating considerable heating occurs. Here we show the calculated results at the BEC transition point: temperature $`T=56\mu `$K, peak density $`\stackrel{~}{n}(\mathrm{𝟎})=2.0\times 10^{14}`$ cm-3, number of trapped atoms $`\stackrel{~}{N}=3.3\times 10^{10}`$, and truncation energy $`ϵ_t/k_\mathrm{B}=302\mu `$K. These parameters are quite consistent with the experiment , which confirms the validity of our numerical calculations.

The condensate formation process shows the characteristic feature depicted in Fig. 2. Condensates first grow very rapidly after the transition. We can explain this by the fact that bosonic stimulation in evaporative cooling process  strongly accelerates the growing speed, since the density of condensates becomes very high in atomic hydrogen. The number of condensates $`N_0`$ reaches the maximum value at 17.2 s ($`\tau =0.69`$), and then decreases gradually due to dipolar decay. The calculated parameters at this maximum point are: temperature $`T=49\mu `$K, peak density of condensates $`n_p=4.3\times 10^{15}`$ cm-3, peak density of noncondensed atoms $`\stackrel{~}{n}_n=1.7\times 10^{14}`$ cm-3, total number of condensates $`N_0=9.0\times 10^8`$, total number of noncondensed atoms $`\stackrel{~}{N}_n=2.0\times 10^{10}`$, and truncation energy $`ϵ_t/k_\mathrm{B}=239\mu `$K. The corresponding condensate fraction at the maximum point of $`N_0`$ is calculated as $`f=N_0/(N_0+\stackrel{~}{N}_n)`$ to be a small value of 4.3%. These results quantitatively agree with the recent observation of BEC with atomic hydrogen .

The number of trapped atoms has been evaluated experimentally by measuring the density-dependent frequency shift of the two-photon 1S-2S transition in atomic hydrogen . The condensate number, $`N_010^9`$ in Ref. , was obtained on the assumption that the frequency shift for condensates is as large as that for noncondensed atoms with the same density. If we expect the disappearance of exchange effects in the excitation of condensates as discussed in Ref., the frequency shift for condensates becomes half of that for noncondensed atoms, and many more condensates, $`N_06\times 10^9`$, with an unreasonably higher condensate fraction, $`f=25\%`$, are given. The calculated results in Fig. 2 support the evaluation of condensate number which assumes the equal frequency shift for condensates and noncondensed atoms.

Next, we discuss the influence of dipolar loss on evaporative cooling. From Eq. (3), the dipolar relaxation process preferably removes atoms in the region with higher density. Since the removed atoms have lower internal energy in the trapping potential, the dipolar relaxation loss results in system heating. The time evolution of the dipolar loss rate $`\dot{N}_{\mathrm{dip}}`$ is shown in Fig. 3(a). The curve exhibits the jump structure around the BEC transition point indicated by the arrow. The serious heating after the transition in Fig 1(b) is caused by such drastic increase of dipolar loss known as “relaxation explosion” in hydrogen system . To clearly demonstrate the competition between evaporation and dipolar relaxation processes during evaporative cooling, the normalized rates $`\dot{N}_{\mathrm{ev}}/\stackrel{~}{N}`$ and $`\dot{N}_{\mathrm{dip}}/\stackrel{~}{N}`$ are plotted together in Fig. 3(b). The evaporation rate is larger than the dipolar loss rate until just after the BEC transition, which indicates that evaporative cooling is sufficiently efficient in this region. The explosively enhanced dipolar-loss rate after the transition finally exceeds the evaporation rate at around 17.5 s, and the resultant heating becomes serious. We note here that this time almost corresponds to the point where the number of condensates shows the maximum value in Fig. 2. The heating-cooling balance therefore limits the producible number of condensates and the corresponding condensate fraction. One can expect that the trapping potential with weaker confinement and the optimized evaporative cooling will move this balance in the direction of a larger condensate production .

In conclusion, we have investigated the rf-induced evaporative cooling of magnetically trapped atomic hydrogen on the basis of the kinetic theory of a Bose gas. This approach can be applied over the whole temperature region of evaporative cooling while the earlier investigations based on the classical kinetic theory  fail around the critical point of BEC transition. The calculated results quantitatively agree with the recent BEC experiment of atomic hydrogen. This proves our calculations very useful, and the important future work would be to optimize the cooling trajectory for a larger production of condensates . Our theory, on the other hand, assumes that the system always stays in the quasi-thermal equilibrium states during evaporative cooling. A precise study of the deviation from the quasi-thermal equilibrium  would give the limit of the approximation in the theory.

We thank L. Willmann of MIT for his valuable suggestions and discussions.


text:
# Three Dimensional Magneto-Hydrodynamical Modeling of Planetary Nebuale II: The Formation of Bipolar and Elliptical Nebulae with Point-Symmetric Structures and Collimated Outflows

## 1 Introduction

Numerical simulations that reproduce the shapes of planetary nebulae (PNe) have helped to understand or give insights into the basic physics acting behind the formation of the different morphologies that are commonly observed in PNe. The two-wind model proposed by Kwok, Purton & Fitzgerlad (1978) has been applied and refined in a number of numerical simulations based on hydrodynamic models (e.g. Khan & West 1985; Mellema, Eulderink & Icke 1991; Frank & Mellema 1994; Dwarkadas, Chevalier & Blondin 1996). During the last decade, though, it has become evident that the production of highly collimated, bipolar, supersonic outflows, play a major role in the evolution of a substantial number, if not all, of PNe (e.g. López 1997 for a review).

The main nebular morphologies are usually well reproduced in the hydrodynamic models, however, the production of jets and ansae have met with certain difficulties; in particular, fast winds do not tend to converge into stable structures to form ansae and jets (e.g. Dwarkadas & Balick 1998).

On a somewhat different approach, Chevalier & Luo (1994) have explored the effects of a rotating star with a magnetized wind on the formation of aspherical bubbles. Following this scheme, Różyczka & Franco (1996) and García-Segura et al. (1997a, 1999), performed 2-D MHD simulations of this magnetized wind in cylindrical and spherical calculations, respectively, showing that magnetic tension may indeed be responsible for the generation of jets in PNe. Recently, García-Segura (1997) (paper I) presented full 3-D MHD models of PNe where jets and ansae are convincingly reproduced as the result of magnetic collimation of the post-AGB wind.

A particularly intriguing case in PNe morphologies are those that display point-symmetric structures. Point-symmetry has been recognized in a wide variety of PNe (e.g. Guerrero, Vázquez & López 1998). In some instances, where reliable kinematical information is available, the presence of bipolar, rotating, episodic jets or BRETS have been inferred to explain similar structures (c.f. López 1997).

Cliffe et al. (1995) performed 3-D numerical simulations introducing in their computational mesh the presence of a BRET and reproducing the main morphological features observed, in particular for the case of Fleming 1 (López, Meaburn & Palmer 1993). They also noted that some bipolar nebulae, such as Hb 5 present distinct, clumpy, S-shaped intensity distributions in their opposing lobes, which led them to argue that these structures could also be originated by a precessing jet. The numerical simulations presented in this paper have been mainly motivated by this group of objects, though the results expand a wide range of configurations, including the conditions for jet detection in PNe, ansae formation in point-symmetric modes and ‘regular’ bipolar and elliptical PNe presenting some degree of asymmetry.

In § 2 we discuss the initial considerations and the motivation for the present paper. § 3 describes the methodology. § 4 presents results for bipolar nebulae and § 5 for ellipticals. § 6 Discusses the results and a summary is given in § 7.

## 2 Initial considerations

As an illustration of the main type of morphological structures that are the focus of the present simulations, a mosaic of images of these type of nebulae is presented in Figure 1. The objects included in this representative sample are, (1a) Hb 5, (1b) He 2-103, (1c) He 2-429 and (1d) NGC 2371. Hubble 5 is a bipolar nebula with distinct internal S-shaped intensity distribution. Recent HST images of this object obtained by Balick, Icke and Mellema show a great deal of detail in the secondary, tilted outflow. He 2-103 is an elliptical nebulae where the equatorial density enhancement is clearly not perpendicular to the elliptical shell; furthermore, faint nebular extensions are discerned displaced from the main elliptical axis on both sides of the nebula. He 2-429 (see Guerrero et al. 1998 and the IAC Catalog, Manchado et al. 1996) is an elliptical nebula with twisted, point-symmetric extensions protruding outside the main nebular shell. Finally, NGC 2371 has a dumb-bell appearance with intensity maxima in the equator defining a minor axis which is not orthogonal to the main nebular axis.

The main, standard, working hypothesis that sets the initial conditions of the simulations is that the equatorial density enhancement (EDE) produced in the late AGB phase defines the main symmetry axis of the bipolar/elliptical wind outflow, i.e. the EDE is orthogonal to the main bipolar/elliptical structure.

In paper I, in order to reproduce point-symmetric structures, a wide binary nucleus was considered and the magnetic collimation axis (i.e. the rotation axis of the primary star) was assumed to precess, with a period of the same order of the kinematic age of the nebulae ($`10^3`$ yr).

In the present paper, instead of invoking an active rotation of the magnetically collimated outflow, a misaligned configuration is produced by introducing a steady tilt of the magnetic collimation axis with respect to the bipolar/elliptical wind outflow. Hence, this is a particular case of paper I, in which the precession period is nule or much larger than the kinematic age of the nebula.

Thus, in order to achieve this steady misalignment, the rotating core of the AGB, which will provide the magnetized fast wind, must have a rotation axis tilted with respect to the axis defined by the EDE. Our current limited understanding of the detailed processes that produce the EDE at the AGB phase precludes a full physical justification for this assumption. Binary nuclei could represent the best way to attain this effect given the diversity of paths that they can take during their evolution through the AGB (e.g. Yungelson, Tutukov & Livio 1993). For example, wide binary systems may induce this misalignment since the ejected envelope’s rotation axis of the primary star could be tidally altered by the secondary and the EDE formed in the orbital plane, with a resulting tilt with respect to the rotation axis of the primary. Alternatively, systems which undergo a common envelope evolution, in which the rotation axis of the primary star is not orthogonal to the orbital plane could also produce the same effect. Close binary systems with substellar companions that may lead to the formation of accretion disks (Reyes-Ruiz & López 1999) and mergers could also play a significant role here.

However, beyond these arguments, which badly need additional investigation, the best initial justification that has motivated at this time further exploring of this assumption is the apparent asymmetry observed in the equatorial regions of a substantial number of PNe and the corresponding point-symmetric condensations and shapes (see Figure 1).

The numerical models do not consider an explicit physical form for the source of the wind, just an outflowing magnetized wind from a spherical radius of $`r=0.005`$ pc after t=0. Regions below this limit lie beyond the current resolution of the simulations.

The general idea behind these MHD models is sketched in Figure 2. Here, a rotating star produces the outflowing magnetized wind, in which the toroidal component of the magnetic field is dominant respect to the poloidal contribution at large distances.

## 3 Three-dimensional MHD models

The simulations have been performed using the magneto-hydrodynamic code ZEUS-3D (version 3.4), developed by M. L. Norman and the Laboratory for Computational Astrophysics. This is a finite-difference, fully explicit, Eulerian code descended from the code described in Stone & Norman (1992). A method of characteristics is used to compute magnetic fields as described in Clarke (1996). ZEUS-3D does not include radiation transfer, but we have implemented a simple approximation to derive the location of the ionization front for arbitrary density distributions (see Bodenheimer et al. 1979; Franco et al. 1989, 1990, García-Segura & Franco 1996). This is done by assuming that ionization equilibrium holds at all times, and that the gas is fully ionized inside the H ii region.

The location of the ionization front in any given direction $`(\theta ,\varphi )`$ from the photoionizing source is given by $`n^2(r,\theta ,\varphi )r^2𝑑rF_{}/4\pi \alpha _B`$, where symbols have their usual meaning. The models include the Raymond & Smith (1977) cooling curve above $`10^4`$ K. For temperatures below $`10^4`$ K, the shocked gas region is allowed to cool down with the radiative cooling curves given by Dalgarno & McCray (1972) and MacDonald & Bailey (1981). Finally, the photoionized gas is always kept at $`10^4`$ K, so no cooling curve is applied to the HII regions (unless, of course, there is a shock inside the photoionized region). The minimum temperature allowed in all models is set to $`10^2`$ K. The computations are done in Cartesian coordinates (x, y, z), with outflowing outer boundary conditions.

The models have grids of $`100\times 100\times 200`$ equidistant zones, with an extent of $`0.1\times 0.1\times 0.2`$ pc in x, y and z respectively. The wind, whose origin is centered at the center of the computational volume, is set based on the rotating wind solutions given by Bjorkman & Cassinelli (1993), Ignace et al. (1996, 1998) and the approach to those equation given by García-Segura et al.(1999). These equations permit the introduction of the wind self-consistently. For the reasons given in the preceding section, in the present calculations we do not address a particular functional shapes for the EDE, but simply state its presence independently of its origin.

The magnetic field of an outflowing wind from a rotating star can be described by two components, $`B_\varphi `$ and $`B_\mathrm{r}`$ (Chevalier & Luo 1994; Różyczka & Franco 1996; paper I; García-Segura et al. 1999). The radial field component can be neglected since $`B_\mathrm{r}r^2`$ while the toroidal component $`B_\varphi r^1`$, and the field configuration obeys $`B=0`$ . In all the simulations, the toroidal field is set always perpendicular to the z axis, therefore only $`B_{\varphi _\mathrm{x}}`$ and $`B_{\varphi _\mathrm{y}}`$ are introduced in the source wind sphere, with $`B_{\varphi _\mathrm{z}}=0`$.

As in paper I, we will use two useful dimensionless parameters containing the basic information of the models. The first one $`\mathrm{\Omega }`$ is the velocity ratio of the stellar rotation to the critical rotation, and the second parameter $`\sigma `$ is the ratio of the magnetic field energy density to the kinetic energy density in the fast wind. The wind compression field (“WCField”) found by Ignace et al.(1998) is also introduced. However, since $`\mathrm{\Omega }`$ is less than 10% for the fast wind in all the calculations, this correction is of little relevance in these cases.

The AGB wind is set with a $`v_{\mathrm{}}=10\text{ km s}\text{-1}`$ and a $`\dot{M}=10^6\text{M}\text{}\text{ yr}\text{-1}`$ as a standard value in all the models. The AGB wind is set unmagnetized in all models, because the estimated values for $`\sigma `$ are very low in this regime (see paper I). Since we are interested in the early phases of the PN evolution, the velocity of the fast wind is set to $`100\text{ km s}\text{-1}`$. This value should reasonably match the conditions at the departure from the AGB phase.

A more realistic approach to this problem would need to include the evolving wind conditions with gradual transitions of the wind from $`v_{\mathrm{}}10`$  km s-1 up to $`1000`$  km s-1 (or higher), and mass-loss rates varying from $`\dot{M}=10^6`$ M yr-1 (or higher) down to $`10^9`$ (or lower). In this case time dependent effects invoke additional behaviors on the dynamics (e.g. Dwarkadas & Balick 1998). They are out of the scope of the present paper and are deferred to a future work.

The parameters of the runs are summarized in Table 1. These include, an unmagnetized reference case (model A0); bipolars with different inclinations of the bipolar wind outflow symmetry axis with respect to the magnetic collimation axis and two different mass loss rates (models B and C, respectively) and ellipticals with a single constant tilt and three different mas loss rates (models E)

## 4 Bipolar Nebulae

Bipolar nebulae are mainly characterized by an equatorial waist with a substantial pole to equator density contrast. In order to achieve such a density contrast, we have used a value of $`\mathrm{\Omega }=0.98`$ in the set up, which is obtained with a rotational velocity of 6.929  km s-1 .

### 4.1 The unmagnetized bipolar case

As a starting point, the response of a spherically symmetric, unmagnetized wind starting from a spherical shape is shown in model A0 (Fig. 3). This unmagnetized case shows how the wind behaves in such a medium, i.e., the wind is impeded to expand in the equator due to the EDE, but it easily expands along the polar direction. This would be the classical hydrodynamic case that has been thoroughly studied in the past (c.f. see Frank 1999 for a review on bipolar outflows). Because of the radiative conditions of the reverse or terminal shock of such a wind in this set up ( $`t_{\mathrm{post}\mathrm{shock}}10^5`$ K for $`v_{\mathrm{shock}}=100\text{ km s}\text{-1}`$), the bipolar bubble is momentum driven. This particular dynamics generates ram-ram pressure instabilities (García-Segura et al. 1996; García-Segura et al. 1997b), however, given the low resolution of the grid, these instabilities are not reproduced in the present study, and only some indications of such a behavior is present as the small corrugation that appears in the lobes of model A0.

Model A0 is important as a reference model since it shows that fast unmagnetized winds, starting from spherical conditions and expanding within an EDE cannot converge to form jets or ansae. In more adiabatic conditions, the lobes become two tangential spheres (García-Segura et al. 1999). It is precisely this behavior that rules under the presence of an EDE, i.e. the bipolar lobes will grow always perpendicular to the EDE, and is the EDE the responsible for the equatorial waist observed in bipolar nebulae.

### 4.2 Magnetized Bipolar Nebulae

Keeping in mind the behavior of the reference model A0, a magnetized wind is now introduced to visualize the strong change that effects the dynamics. Model B0 (Fig. 3) only differs from model A0 in the magnetized wind, with a $`\sigma =0.01`$. But the result is remarkable, a pair of jets form at both sides in the lobes. The explanation of this behavior is fully described in Różyczka & Franco (1996), and it results mainly from the action of the hoop stress caused by the toroidal magnetic field. The rest of the bubble behaves now adiabatically because the magnetic pressure does not allow the gas to pile up towards the inner edges of the lobes. This is also the main reason why model B0 expands much faster than model A0, since energy conserving bubbles in power laws of the type $`r^2`$ have different dynamics than momentum conserving bubbles (García-Segura & Mac Low 1995 and references therein).

Although the hoop stress is always present in this type of flows, not always the conditions are in favor to form observable jets. Such a case is computed in model C0 (Fig. 3), in which the mass-loss rate of the fast wind has been decreased by an order of magnitude. In such conditions, the fast wind piles up at the poles, but the radiative cooling, proportional to $`n^2`$, is not sufficiently strong to allow the formation of visible jets, and ansae-like structures are formed instead. This difference is more evident in the synthetic pictures shown in figure 7 for models Bs and Cs, which are discussed further below.

As mentioned before, the symmetry axis of the bipolar wind outflow is considered orthogonal to the EDE. Models B0 and C0 are computed with the magnetic collimation axis also orthogonal to the EDE. The following models break this condition by introducing different tilts between these axes to explore the resultant effects, which are shown to lead to point-symmetric shapes.

The tilt is introduced at the EDE for simplicity. Thus, the magnetic collimation axis remains parallel with the z axis, and the same set up described above is used for the magnetized wind.

In order to have a general picture of the effect caused by the tilt, we have computed several models in which the tilt angle varies from low values, $`5^{}`$ (models B5 and C5), middle values, $`15^{}`$ (models B15 and C15) and large values, $`45^{}`$ (models B45 and C45). The evolution of the two last models is shown in figures 4 and 5, respectively.

The dynamics of the bubbles in such conditions are quite complex and the limited resolution of the models hinder a full, quantitative description. Nevertheless, the effect of the tilt can be summarize qualitatively as follows.

At the beginning of every simulation, the spherically symmetric shocked flow is driven by the magnetic pressure gradient after crossing the reverse or terminal shock and is launched in the z direction on both sides of the wind source region. Since the magnetic field is toroidal, the magnetic pressure is nule towards the z axis in the free expanding wind (but not in the shocked region), while it is maximum in the equatorial plane, consequently building a large pressure gradient. Gradually, the flow notices the effect of the EDE and is deflected perpendicular to it. The magnetically pressurized bubble is distorted, and the flow, which now follows the tendency of the bipolar cavity, pushes and deflects the column of gas in which the hoop stress is larger. When this column is dense enough, the formed jets are deflected an acquire the appearance of precessing jets, which eventually will collide with the bipolar cavity (see e,g, model B45). Thus, there are two competing effects, one is the tendency of the gas to flow perpendicular to the EDE, and the other is the tendency of the gas to flow parallel to the z axis as imposed by the magnetic pressure and the hoop stress. When the tilt is nule, i.e., when both effects do not compete with each other, the bubble acquires the largest expansion velocity at the polar regions.

## 5 Magnetized Elliptical Nebulae

Elliptical nebulae are characterized by a moderate pole to equator density contrast. We have adopted a value of $`\mathrm{\Omega }=0.80`$ for the AGB wind, which is obtained with a rotational velocity of 5.65  km s-1 . The AGB wind is set with a $`v_{\mathrm{}}=10\text{ km s}\text{-1}`$ and a $`\dot{M}=10^6\text{M}\text{}\text{ yr}\text{-1}`$ as a standard value for models E1 and E2, while $`\dot{M}=10^8\text{M}\text{}\text{ yr}\text{-1}`$ is used for E3 and E4. The mass loss rate of the fast wind is varied from values of $`\dot{M}=10^7\text{M}\text{}\text{ yr}\text{-1}`$ (E1) to $`\dot{M}=10^9\text{M}\text{}\text{ yr}\text{-1}`$ (E3 and E4). A tilt of $`15^{}`$ is used in this case for all models for simplicity.

An example of the model evolution is shown in Fig. 6 for model E1. The evolution is similar to the models described above for bipolar nebulae, with the difference that the bubble grows further at the equator since the EDE is not as restrictive as before. Figure 7 (right panels) shows the apparent differences between models E1, E2, E3 and E4, dictated by varying the mass loss rates of the fast wind from relatively large to lower values. The polar, piled-up gas in model E1 ($`\dot{M}=10^7\text{M}\text{}\text{ yr}\text{-1}`$) is able to cool down efficiently in a dynamical time scale comparable to the computed time and the axial jets become readily apparent. In model E2 ($`\dot{M}=10^8\text{M}\text{}\text{ yr}\text{-1}`$) the cooling is less efficient but still present, giving rise to the formation of ansae-type structures. Models E3 and E4 ($`\dot{M}=10^9\text{M}\text{}\text{ yr}\text{-1}`$) do not show any apparent signs of piled-up gas (see Figure 7). It should be noted, however, that given the limited numerical resolution of these models, the cooling regions are not resolved and the mass loss rates should be taken as estimative values. However, since these models have energy conserving dynamics (as imposed by the magnetic pressure) rather than momentum conserving, these estimations for the mass loss rates are quite independent of the terminal wind velocities. In the absence of a magnetic field, winds of the order of $`100\text{ km s}\text{-1}`$ have strongly radiative terminal or reverse shocks (see for example Frank et al. 1996), and the radiative conditions are strongly dependent on the wind velocity and not only of the mass loss rate.

## 6 Analysis of the models

In order to test the numerical results of the models, it would be convenient to run convergence tests at lower and higher resolution. The latter becomes prohibitively demanding for our computational resources. An experiment at half the resolution, which gives a qualitative indication on the reliability of the model is presented in Figure 6b. This figure shows a comparison for model E1 with two different resolutions. Althought, the individual details in both simulations in that figure differ by obvious reasons, the main features, such as the shape of the bubble and the formation and extension of the jets remain rather similar. The densities are larger in the high resolution run, indicating that in order to make quantitative estimates of physical parameters, one needs twice the resolution at least to be sure that the dynamics has been resolved, i.e., to find a convergence. For example, what is the post-shock cooling length and how does this compare with the resolution of the simulations ?

In order to resolve a cooling distance, one needs basically to get $`d_{\mathrm{cool}}>x_{\mathrm{max}}/n_\mathrm{x}`$ , where $`x_{\mathrm{max}}`$ is the physical size of the grid, and $`n_\mathrm{x}`$ the number of zones (i.e. $`d_{\mathrm{cool}}`$ bigger than one zone). However, several zones are needed to resolve detailed cooling processes (e.g., any sort of microphysics). A rough estimate of the cooling timescales, can be obtained from equilibrium line cooling. In this case, the cooling time of a parcel of shocked gas is given by $`t_{\mathrm{cool}}=3^{3/2}v_\mathrm{s}^3/(2^8A\rho _0)0.3v_8^3n_6^1`$ yr (Franco et al. 1994), where $`v_8=v_\mathrm{s}/10^8`$ cm s-1 is the shock velocity, $`\rho _0`$ is the pre-shock mass gas density and $`n_6`$ is this density in units of $`10^6`$ cm-3, $`A=3.9\times 10^{32}q_\mathrm{p}^{5/2}`$ cm6 s-4 g-1, and $`q_\mathrm{p}`$ is the mass per particle in units of $`m_\mathrm{H}`$ (i. e., $`q_\mathrm{p}=1.4/2.3`$ for a fully ionized gas with cosmic abundances). For typical numbers in our simulations ($`v_\mathrm{s}100\text{ km s}\text{-1}`$, $`n10`$ cm-3) we obtain $`t_{\mathrm{cool}}30`$ yr, which give us $`d_{\mathrm{cool}}2\times 10^{15}`$ cm, just slightly smaller than the zone zise of $`3\times 10^{15}`$ cm. One should also consider that in MHD models conditions are such that cooling distances are more favourable resolved perpendicular to the field lines given the tendency to adiabatic behaviors in this case.

Since the MHD simulations behave similar to the hydro solutions along the Z axis (parallel to the field lines), except for the extra pressure (hoop stress) that the magnetic field produces towards the rotation axis of the central star, it is possible to use the relations of conical converging flows (Cantó, Tenorio-Tagle & Różyczka 1988; or the specific solutions for PPNe by Frank, Balick & Livio 1996), as an approximation for the parameters in the jets. This is in fact a conical converging, magnetized flow problem. Thus, for example, the density in the jets can be estimated by following formula (23) in Cantó, Tenorio-Tagle & Różyczka (1988), where $`\rho _\mathrm{j}6\times 10^3sin^4\theta \left[v_0/100\text{ km s}\text{-1}\right]^4\rho _0`$ . For the case of our E1 model, one gets (for $`\theta 50^{}`$) $`\rho _\mathrm{j}1\times 10^3\rho _0`$. On the other hand, the computed models give $`\rho _\mathrm{j}1\times 10^1\rho _0`$ for the low resolution run and $`\rho _\mathrm{j}1\times 10^2\rho _0`$ for the higher resolution one. Since the resolution is not big enough (as discused above), and in both cases, $`\rho _01\times 10^{23}`$  g cm-3 , it is expected to find a jump of three orders of magnitude on the density ratio at twice the resolution, which would match the analytical value.

## 7 Discussion

The computational survey reveals that the formation of jets and ansae, or none, can be achieved by varying the mass-loss rate within a common framework of magnetized winds. In principle this would allow to have an estimate of the mass-loss rates involved in the formation of the different types of collimated outflows found in PNe. In particular, it is found that for the formation and detection of jets, the mass-loss rates must have values $`\dot{M}10^7\text{M}\text{}\text{ yr}\text{-1}`$. In the case of ansae a $`\dot{M}10^8\text{M}\text{}\text{ yr}\text{-1}`$ is required. Below these values, i.e. for typical cases of $`\dot{M}10^9\text{M}\text{}\text{ yr}\text{-1}`$ or lower, no collimated features are formed or revealed by the models. These numbers, however, must be taken with caution given the low resolution used in these 3-D calculations. Nevertheless, there are additional studies that indicate that jets in PNe must be formed during a brief lapse of high mass-loss rate in the pre-planetary nebula stage. For example, Reyes-Ruiz & López (1999) have concluded that when an accretion disk is involved in the formation of collimated outflows in PNe, the necessary high accretion rates, $`\dot{M}10^7`$, needed to launch an axis-symmetric wind (as in the case of young stellar objects, e.g. Calvet 1998) are attained only for a few thousand years during the late AGB or pre-planetary phase.

Figure 7 shows different models in which jets have formed. The projected shapes in the tilted cases give the impression of a varying direction in the outflow, closely resembling the point-symmetric structures often observed in a number of PNe. However here, no active precession or rotation in the source of the outflows is involved. These shapes result by the steady misalignment of the magnetized wind’s axis with respect to the axis of the bipolar wind outflow. This configuration can lead to a hydrodynamical deflection of the magnetized wind on the bipolar cavity. This is a novel feature revealed by the simulations and which can be tested in the future by detailed kinematical and spectropolarimetric studies. For example, in this case the deflected gas is expected to show a kinematic behavior that deviates from the rest of the nebular expansion and line emission in the magnetically collimated regions to show polarization effects. Hydrodynamical deflection can only account for S-shaped morphologies, helicoidal type structures would require some sort of rotation or wobbling from the wind source.

Although the PNe that were originally classified as point-symmetric comprise a very reduced group, the presence of point-symmetric structures in PNe has been shown now to span a much wider sample (Guerrero, Vázquez & López 1998). For example, signatures of point-symmetry can be appreciated in nearly 20% of the sample contained in the IAC morphological catalog (Manchado et al. 1996). Furthermore, practically every PN that has been recently observed at high spatial resolution with the Hubble Space Telescope shows some degree of point-symmetry, indicating that the mechanisms that produce it should be rather common. The representative models shown in Figures 2 - 7 reveal the variety of point-symmetric configurations that can be formed under the assumptions and conditions of the present study. The inclusion of binary nuclei in this analysis can only but strengthen the conditions for the development of point-symmetry.

## 8 Summary

We have made a computational survey of 3-D MHD simulations of young planetary nebulae in which are introduced different steady misalignments of the magnetic collimation axis with respect to the symmetry axis of the bipolar/elliptical wind outflow. The latter defined as perpendicular to the equatorial density enhancement. This steady tilt can also be interpreted as the result of a very large precession period. The simulations show that in these cases a hydrodynamical deflection of the magnetized collimated wind on the bipolar/elliptical cavity can produce morphologies which may resemble the presence of precessing or rotating sources. In this way, the resulting morphologies are able to reproduce a considerable number of PNe which present some degree of point-symmetry in their structures. The survey also suggests that the formation and detection of jets and ansae, or none, strongly depends on the mass-loss rate of the driving, magnetized wind. In particular, the presence of jets is detected in the models only for mass-loss rates $`\dot{M}10^7\text{M}\text{}\text{ yr}\text{-1}`$; below these rates the jets tend to vanish in the models, keeping for $`\dot{M}10^8\text{M}\text{}\text{ yr}\text{-1}`$ only ansae type structures.

Acknowledgements

GGS is grateful to J. Franco, N. Langer, M.-M. Mac Low, A. Manchado and M. Peimbert for useful discussions. We thank our anonymous referee for his comments. We also thank Michael L. Norman and the Laboratory for Computational Astrophysics for the use of ZEUS-3D. The computations were performed at the Instituto de Astronomía-UNAM. This work has been supported by grants from DGAPA-UNAM (IN 114199) and CONACyT (32214-E).

List of References

Bjorkman, J. E., & Cassinelli, J. P. 1993, ApJ, 409, 429

Bodenheimer, P., Tenorio-Tagle, G. & Yorke, H. W. 1979, ApJ, 233, 85

Calvet, N. 1998, in Accretion Processes in Astrophysical Systems: Some like it hot!, ed. S.S. Holt & T.R. Kallman, (AIP press: New York), p. 495

Cantó, J., Tenorio-Tagle, G., & Różyczka, M. 1988, A&A, 192, 287

Chevalier, R. A., & Luo, D. 1994, ApJ, 421, 225

Cliffe, J. A., Frank, A., & Jones, T.W. 1996, M.N.R.A.S., 282, 1114

Clarke, D. A. 1996, ApJ, 457, 291

Dalgarno, A. & McCray, R. A. 1972, ARAA, 10, 375

Dwarkadas, V. V., Chevalier, R. A. & Blondin, J.M. 1996, ApJ, 457, 773

Dwarkadas, V. V. & Balick, B. 1998, ApJ, 497, 267

Franco, J., Miller, W. W., Arthur, S. J., Tenorio-Tagle, G. & Terlevich, R. 1994, ApJ, 435, 805

Franco, J., Tenorio-Tagle, G. & Bodenheimer, P. 1989, RMxAA, 18, 65

Franco, J., Tenorio-Tagle, G. & Bodenheimer, P. 1990, ApJ, 349, 126

Frank, A. 1999, New AR, 43, 31

Frank, A. & Mellema, G. 1994, A&A, 289, 937

Frank, A., Balick, B., & Livio, M. 1996, ApJ, 471, L53

García-Segura, G. 1997, ApJL, 489, L189 (paper I)

García-Segura, G., & Franco, J. 1996, ApJ, 469, 171

García-Segura, G., Langer, N., and Mac Low, M.-M. 1997b, ”LBV Outbursts: The Effects of Rotation”, in A.S.P. Conference Series, Luminous Blue Variables: Massive Stars in Transition, ed. A. Nota and H. Lammers, Vol 120, 332

García-Segura, G., Langer, N., Różyczka, M., & Franco, J. 1999, ApJ, 517, 767

García-Segura, G., Langer, N., Różyczka, M., Mac Low, M.-M., & Franco, J. 1997a, in IAU Symp. No 180 “Planetary Nebulae”, eds. H.J. Habing & H.J.L.M. Lamers, (Kluwer Academic Publishers, Dordrecht), 226

García-Segura, G., Mac Low, M.-M., & Langer, N. 1996, A&A, 305, 229

García-Segura, G., & Mac Low, M.-M. 1995, ApJ, 455, 145

Guerrero, M.A., Vázquez, R., & López, J.A. 1998, ApJ, 117, 967

Ignace, R., Cassinelli, C. P., & Bjorkman, J. E. 1996, ApJ, 459, 671

Ignace, R., Cassinelli, C. P., & Bjorkman, J. E. 1998, ApJ, 505, 910

Khan, F.D. & West, K. 1985, MNRAS, 212, 837

Kwok, S., Purton, C. R. & Fitzgerald, P. M. 1978, ApJ 219, L125

López, J. A. 1997, in IAU Symp. No 180 “Planetary Nebulae”, eds. H.J. Habing & H.J.L.M. Lamers, (Kluwer Academic Publishers, Dordrecht), 197

López, J.A., Meaburn, J. & Palmer, J. 1993, ApJ, 415, L135

MacDonald, J. & Bailey, M. E. 1981, MNRAS, 197, 995

Manchado, A., Guerrero, M., Stanghellini, L., & Serra-Ricart, M. 1996, “The IAC Morphological Catalog of Northern Galactic Planetary Nebulae”, ed. Instituto de Astrofísica de Canarias

Mellema, G., Eulderink, F. & Icke, V. 1991, A&A, 252, 718

Raymond, J. C. & Smith, B. W. 1977, ApJS, 35, 419

Reyes-Ruiz & López 1998, ApJ, 524, 952

Różyczka, M. & Franco, J. 1996, ApJL, 469, L127

Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 753

Yungelson, L. R., Tutukov, A. V., & Livio, M. 1993, ApJ, 418, 794

Figure Captions

Figure 1. A representative sample of the type of nebulae that are the focus of this study. 1(a) Hb 5; 1(b) He 2-103; 1(c) He 2-429; 1(d) NGC 2371. See models B45, E3, E1 and E4 in Figure 7, respectively for comparison.

Figure 2. Schematic view of the models. The magnetic collimation axis (dashed line) is tilted with respect to the symmetry axis of the bipolar wind outflow, which is perpendicular to the equatorial density enhancement.

Figure 3. Evolution of the logarithm density distribution (g cm-3) of a two-dimensional slice in the x-z plane which crosses through the middle of the computational volume. Top panels correspond to the unmagnetized case, model A0 ($`\dot{M}=10^7\text{M}\text{}\text{ yr}\text{-1}`$) for t = 1250 and 2500 yr. Middle panels correspond to the magnetized case, model B0 ($`\sigma =0.01`$ and $`\dot{M}=10^7\text{M}\text{}\text{ yr}\text{-1}`$) for t = 900 and 1800 yr. Bottom panels, magnetized case, model C0 ($`\sigma =0.01`$ and $`\dot{M}=10^8\text{M}\text{}\text{ yr}\text{-1}`$) for t = 1250 and 2500 yr. In all these cases the symmetry axis of the bipolar wind outflow and the magnetic collimation axis are coincident.

Figure 4. Evolution of the logarithm density distribution (g cm-3) of a two-dimensional slice in the x-z plane which crosses through the middle of the computational volume for model B45 ($`\sigma =0.01`$ and $`\dot{M}=10^7\text{M}\text{}\text{ yr}\text{-1}`$). Here the symmetry axis of the bipolar wind outflow is tilted by $`45^{}`$ with respect to the magnetic collimation axis (t = 350, 700, 1000, 1350, 1650, 2000 yr).

Figure 5. Same as Figure 4 but for model C45 where $`\dot{M}=10^8\text{M}\text{}\text{ yr}\text{-1}`$ (t = 400, 850, 1250, 1650, 2100, 2500 yr).

Figure 6. Same as Figure 4 for model E1 ($`\sigma =0.01`$ and $`\dot{M}=10^7\text{M}\text{}\text{ yr}\text{-1}`$). Here the symmetry axis of the bipolar wind outflow is tilted by $`15^{}`$ with respect to the magnetic collimation axis (t = 250, 550, 800, 1050, 1350, 1600 yr).

Figure 6b. Two-dimensional slices of Model E1 for $`50\times 50\times 100`$ (top) and for $`100\times 100\times 200`$ (bottom) equidistant zones in x, y and z respectively.

Figure 7. Two-color, composite, synthetic images of the models where red represents the shock-dominated (non-photoionized) material and green the photoionized gas. The line of sight corresponds to the y axis, being the x-z plane the plane of the sky. Bipolar models in the left panels show the magnetized case with no tilt between the symmetry axis of the bipolar wind outflow and the magnetic collimation axis (model B0) and the magnetized cases for three different tilts (B5, B15 and B45). In these models $`\dot{M}=10^7\text{M}\text{}\text{ yr}\text{-1}`$ and $`\sigma =0.01`$. The central panels correspond to the bipolar models, as before, but with $`\dot{M}=10^8\text{M}\text{}\text{ yr}\text{-1}`$. Elliptical models are shown in the right panels. A tilt of $`15^{}`$ is present in these four models. E1 ($`\dot{M}=10^7\text{M}\text{}\text{ yr}\text{-1}`$, $`\sigma =0.01`$), E2 ($`\dot{M}=10^8\text{M}\text{}\text{ yr}\text{-1}`$, $`\sigma =0.01`$), E3 ($`\dot{M}=10^9\text{M}\text{}\text{ yr}\text{-1}`$, $`\sigma =0.01`$) and E4 ($`\dot{M}=10^9\text{M}\text{}\text{ yr}\text{-1}`$, $`\sigma =0.99`$). Each run in the figure corresponds to: B0 and B5, 1800 yr; B15 and B45, 2000 yr; C0, C5, C15 and C45, 2500 yr; E1, 1600 yr; E2, 1550 yr; E3, 2000 yr and E4, 350 yr.