text:
# Deuterium-burning in substellar objects

## 1 Introduction

The search for very low mass stars and substellar objects (SSOs) progresses nowadays at an astonishing pace. Many brown dwarfs (BDs) have now been discovered either in the field or in young clusters, confirming their predicted general evolution in various observable diagrams (see e.g. Basri 2000; Chabrier & Baraffe 2000 for recent reviews). This evolution relies on both the properties of dense degenerate matter in the interior, which determines the mechanical properties (mass, radius, central density and temperature), and on the complex thermochemistry and radiation processes taking place in the atmosphere, which characterize the thermal properties (effective temperature, spectral distribution) and thus the observable signature, magnitude and color. In the meantime several tens of exoplanets orbiting G-M stars have now been discovered by Doppler techniques, spaning a range from about 11 Jupiter mass ($`0.011M_{}`$) down to about a Saturn mass ($`3\times 10^4M_{}`$). By definition SSOs never reach thermal equilibrium since their core temperature is not hot enough to sustain hydrogen-burning and prevent further gravitational contraction. For this reason they do not define a main sequence, as stars do, and age adds another degree of freedom besides mass for a given observational signature.

To circumvent this difficulty, the abundance of lithium in the atmosphere of low-mass objects has been suggested as an independent diagnostic to determine the mass and/or age of an object below the bottom of the main sequence, the so-called ”lithium-test” (Rebolo et al. 1992; Basri et al. 1996). The temperature for lithium-burning through the $`{}_{}{}^{7}Li(p,\alpha )^4He`$ reaction is lower than the one required for proton fusion. Objects below 0.06 $`M_{}`$ (Chabrier & Baraffe 1997) are not massive enough - and thus not hot enough - to host lithium burning and thus retain their initial lithium abundances. Furthermore, these objects are fully convective and convection reaches high in the atmosphere, so that the observation of atmospheric lithium through the Li-I line reflects the core abundance, since the evolutionary timescale is orders of magnitude larger than the convection timescale. Therefore, the observation of lithium in an object older than $``$108 yrs - the timescale for the destruction of lithium in the lowest mass stars - signifies the lack of hydrogen burning. A complementary outcome of the lithium test is the determination of the age of open clusters by finding the lithium-depletion boundary which separates objects without lithium from slightly fainter objects with lithium (see e.g. Basri 1998, Stauffer et al. 1998).

As mentioned above, the timescale for the destruction of lithium in the lowest mass stars is $``$108 yr so that the lithium-test as a direct identification of the substellar nature of an object can not be applied to clusters younger than this age (although it can still be used to locate the substellar boundary in age or luminosity). Deuterium-burning, which occurs at a cooler temperature (see below) provides a new test at younger ages ($`<\mathrm{\hspace{0.17em}5}\times 10^7`$ yr, see §2) to determine the age and/or mass of an object (Béjar, Zapatero Osorio & Rebolo 1999). From a more fundamental viewpoint, deuterium-burning has been suggested to play a key role in the formation of star-like objects (Shu, Adams & Lizano 1987), following the collapse of a protostellar cloud, so that D-burning may define the minimum mass of brown dwarfs and thus of star-like object formation.

The aim of the present letter is to calculate the depletion of deuterium along evolution for different SSO masses, i.e. to derive \[2D\]/\[2D0\]-age-mass relations, where \[2D0\]=$`2\times 10^5`$ is the initial deuterium mass-fraction, and to relate this quantity to observable quantities, magnitude and color, in various optical and infrared passbands. We also characterize the spectroscopic signature of deuterium in the spectra of these objects.

## 2 Results

The basic theory entering the present calculations has been described in detail elsewhere (Chabrier & Baraffe 1997; Baraffe et al 1997, 1998). The atmosphere models are the so-called NextGen models (Hauschildt, Allard & Baron 1999). Evolutionary models based on these inputs successfully reproduce observations in various diagrams (color-magnitude, color-color, mass-magnitude) for main-sequence stars (see references above) and pre-main-sequence stars (PMS) (Baraffe et al., 2000). The surface gravity $`g`$ of all the objects of interest in the present study (SSOs at $`t<\mathrm{\hspace{0.17em}10}^8`$ yr) is always larger than $`g=10^3`$ cm s-2 so that the atmospheric structure is not affected by sphericity effects. The central density and temperature are in the range $`\rho _c1100`$ g cm-3 and $`T_C5\times 10^510^6`$ K. Under these conditions correlation effects between ions in the plasma largely dominate the kinetic contribution and lower the Coulomb barrier between two fusing particles, yielding an enhancement of the reaction rate, as identified initially by Schatzman (1948) and Salpeter (1954). Moreover, for conditions characteristic of SSO interiors, the electron screening length is of the order of the interionic spacing, so that the ion-ion interaction is screened by the polarized electron background (Chabrier & Baraffe 1997). For the mass-range of interest, the electronic and ionic screening factors are indeed of comparable importance. The present calculations include such screening corrections arising from both ion and electron polarization to thermonuclear reaction rates of D-burning through the $`{}_{}{}^{2}D(p,\gamma )^3He`$ reaction. The present calculations, based on the afore-outlined theory, represent significant improvements upon previous calculations of D-burning in low-mass stars, based in particular on approximate treatments of the atmosphere (e.g. D’Antona & Mazzitelli, 1994).

Calculations for PMS stars have been performed over the mass range $`0.010M_{}M/M_{}1.2M_{}`$, for solar metallicity, for ages $`t10^6`$ yr (Baraffe et al., 2000). In the present paper, we focus on objects near the D-burning limit, i.e. $`0.012M_{}M/M_{}0.1M_{}`$. At such young ages, central D-burning is followed by a rapid contraction of the object and thus a significant variation of the surface gravity within a few Myr (see Table 1). For the characteristic range of effective temperatures $`1800\mathrm{K}<T_{\mathrm{eff}}<\mathrm{\hspace{0.17em}3000}`$ K, convection reaches, or extends close to, the atmospheric layers where deuterated water forms. Therefore, the observed D-abundance should reflect the bulk abundance of the object. As shown recently by Chabrier et al. (2000) the core of massive and old brown dwarfs becomes dense and cold enough for electronic conduction to become the main energy transport mechanism. The onset of conduction, however, occurs only for BDs older than 1 Gyr, and thus does not occur under conditions characteristic of the deuterium-test.

As for hydrogen, the ignition temperature translates into a deuterium-burning minimum mass (DBMM). With the afore-mentioned screening factors, these are found to be $`T_D4.0\times 10^5`$ K and $`m_D0.012M_{}`$, in agreement with previous calculations (Saumon et al. 1996; Burrows et al. 1997). Objects above the DBMM destroy deuterium on a very short time scale in the age range $``$ 4 - 50 Myr, for $`m0.070.015M_{}`$, respectively, whereas objects below this mass preserve their original deuterium abundance.

Table 1 displays the age for which the initial D-abundance has been decreased by a factor 2 (50% depletion) and 100 (99% depletion) in the mass-range of interest, with the corresponding gravity, effective temperature, luminosity and absolute magnitudes. Figure 1 illustrates the D-depletion in the K-band with the corresponding effective temperature. Superposed are the evolutionary tracks of several masses for SSOs, above the DBMM. As illustrated on this figure and Table 1, knowledge of deuterium-depletion provides a powerful determination of the mass of an object, for a given age. Conversely, the position of the 50%-depletion boundary (undistinguishable from the 0%-depletion boundary), and the inferred magnitude, provide an independent method of age determination for clusters younger than $`50`$ Myr, the deuterium-burning timescale of a 0.015 $`M_{}`$ BD, just above the DBMM. It is worth noting that for an age $`t1`$ Myr, this boundary corresponds to the hydrogen-burning minimum mass, $`m_{HBMM}0.07M_{}`$ (Chabrier & Baraffe 1997), so that detection of deuterium in an object older than this age ensures its substellar nature. This provides a very useful diagnostic in very young clusters where the lithium-test can not be applied.

Figure 2 displays the borderline for 50% and 99% D-depletion in a $`L`$-$`T_{\mathrm{eff}}`$ diagram, with isochrones of 1 to 50 Myr. The crossing point between a depletion curve and an isochrone indicates the location of the (L,$`T_{\mathrm{eff}}`$) point which corresponds to the corresponding amount of D-depletion for this age. If we choose the 99%-depletion curve, objects cooler and fainter than the crossing point for a given age will all preserve deuterium.

The observable spectroscopic signature of D-depletion is a challenging task. The most distinct features of the presence of deuterium occur in the range $`\lambda =2.8`$ \- 10 $`\mu `$m. Figure 3 shows the relative differences between a spectrum with deuterium set to cosmic abundance and a spectrum with full D-depletion in this wavelength-range for $`T_{\mathrm{eff}}=2500`$ K, $`\mathrm{log}g=3.5`$ for a solar composition. As shown in the figure, spectral signatures of deuterium are in principle observable. Given the lack of accurate linelists, however, it is presently difficult to quantify exactly the amount of flux absorption due to the deuterium molecules transitions (HDO, OD). The present calculations were done using the most recent AMES linelists (Partridge & Schwenke 1997). The situation should improve significantly in the near future and will yield more precise predictions. Spectral resolution of 30000 and very high signal-to-noise ratio are required to observe this signature (assuming rotation of the object is small enough). Such conditions are within reach with NIRSPEC on KeckII and with instruments planed for the VLT. The observation of bands of deuterated methane CH3D at 1.6 $`\mu `$m could also be used in principle to identify the presence of deuterium. Methane, however, can not be used for the deuterium test as an age-indicator: at the temperature methane forms ($`T<\mathrm{\hspace{0.17em}1800}`$ K), SSOs above the D-burning minimum mass are old enough for all deuterium to have been burned.

## 3 Conclusion

In this Letter, we have derived deuterium-depletion abundances as a function of age for various masses in the substellar domain, and we have characterized the effective temperature, luminosity and observable quantities - absolute magnitudes and colors - corresponding to these depletion factors. We have also characterized the deuterium-burning minimum mass, temperature and luminosity as well as the spectral signatures of the presence of deuterium. These results are based on state-of-the-art, consistent evolutionary calculations between the interior structure, the atmosphere profile and the emergent spectrum. They provide the necessary ingredients to apply the deuterium test to identify the age or mass of young substellar objects and to provide an independent age determination of very young clusters. They also serve to identify objects at the deuterium-burning limit, as the lithium-test has been used to characterize objects at the hydrogen-burning limit. The unambiguous identification of free floating objects below this deuterium-burning minimum mass would indeed suggest that the formation of star-like objects extends below the deuterium-burning minimum mass. Conclusions when applying this test, however, must be drawn with caution, for exoplanets formed in multiple systems may be ejected from the system and end up as free floating bodies.

Note: For the transformation of observables into an $`L`$-$`T_{\mathrm{eff}}`$ diagram, a consistent solution is to use the NextGen grid of synthetic spectra. Spectral type vs $`T_{\mathrm{eff}}`$ relations for young objects are not reliable and using such methods will yield results of dubious validity. The NextGen grid is available by anonymous ftp. Please contact the authors for access details.

Acknowledgments: One of us (PH) acknowledges an invited visiting professor position supported by the Pôle Scientifique de Modélisation Numérique at ENS-Lyon. Part of this work was supported by NSF grant AST-9720704, NASA ATP grant NAG 5-8425 and LTSA grant NAG 5-3619, as well as NASA/JPL grant 961582 to the University of Georgia, NASA LTSA grant NAG5-3435 and NASA EPSCoR grant NCCS-168 to Wichita State University. Some of the calculations presented in this paper were performed on the CNUSC IBM SP2, the IBM SP2 and SGI Origin 2000 of the UGA UCNS, on the IBM SP of the San Diego Supercomputer Center (SDSC), with support from the National Science Foundation, and on the Cray T3E of the NERSC with support from the DoE. We thank all these institutions for a generous allocation of computer time.

## References

Baraffe, I., Chabrier, G., Allard, F., Hauschildt, P.H., 1997, A&A, 327, 1054

Baraffe, I., Chabrier, G., Allard, F., Hauschildt, P.H., 1998, A&A, 337, 403

Baraffe, I., Chabrier, G., Allard, F., Hauschildt, P.H., 2000, in preparation

Basri, G., 1998, in Cool stars in clusters and associations: magnetic activity and age indicators, Mem. Soc. Astr. Ital., 68, 917

Basri, G., Marcy, G.W.M. and Graham, J.R., 1996, ApJ, 458, 600

Basri, G., 2000, ARA&A, 38

Béjar, V.J.S., Zapatero Osorio, M.R., Rebolo, R., 1999, ApJ, 521, 671

Burrows, A., Marley, M., Hubbard, W.B., Lunine, J.I., Guillot, T., et al. 1997, ApJ, 491, 856

Chabrier, G., & Baraffe, I., 2000, ARA&A, 38, 339

Chabrier, G., & Baraffe, I., 1997, A&A, 327, 1039

Chabrier, G., Baraffe, I., Allard, F., Hauschildt, P.H., 2000, ApJ, in press

D’Antona, F., & Mazzitelli, I., 1994, ApJS, 90, 467

Hauschildt, P.H., Allard, F., Baron, E., 1999, ApJ, 512, 377

Partridge, H., & Schwenke, D.W., 1997, J. Chem. Phys., 106, 4618

Rebolo, R., Martín, E.L., Magazzù, A, 1992, ApJ, 389, 83

Salpeter, E.E., 1954, Austr. J. Phys., 7, 373

Saumon, D., Hubbard, W.B., Burrows, A., Guillot, T., Lunine, J.I., Chabrier G., 1996,ApJ, 460, 993

Schatzman, E., 1948, J. Phys. Rad. 9, 46

Shu, F.H., Adams, F.C., Lizano, S., 1987, ARA&A, 25, 23

Stauffer, J.R., Schultz, G., Kirkpatrick, J.D., 1998, ApJ, 499, L199

FIGURE CAPTIONS

Figure 1 : Deuterium-depletion curves as a function of age in the K-band. The dashed lines correspond to different SSO masses, while the solid lines correspond to the 50% and 99% D-depletion limit, respectively. The inset displays the corresponding $`T_{\mathrm{eff}}(t)`$.

Figure 2 : $`L`$-$`T_{\mathrm{eff}}`$ diagram for objects above the D-burning mass. Solid lines correspond to the 50% (right curve) and 99% (left curve) D-depletion limit, respectively, dotted lines to isochrones (in Myr) and dashed lines to cooling sequences for several masses.

Figure 3 : Relative difference in per-cent between infrared spectra calculated for $`T_{\mathrm{eff}}=2500`$K and $`\mathrm{log}(g)=3.5`$ models with complete deuterium depletion and no deuterium depletion $`=(F_0F_D)/(F_0+F_D)`$, where $`F_D`$ is the flux with solar D-abundance and $`F_0`$ is the flux with zero D-abundance. The spectral resolution is about 30,000 (reduced from a computed resolution of about 300,000). The vast majority of differences are due to HDO lines seen on top of a background of H2O lines. This figure was generated using the AMES linelists for water and HDO, which have the currently largest density of lines.


format: text

text:
# Fermilab-Conf-00/229Sept2000 Charged Kaons at the Main Injector (CKM)

## 1 The Physics

The branching ratio for the decay $`K^+\pi ^+\nu \overline{\nu }`$ was first calculated by Inami and Lin  as an isospin rotation from $`K^+\pi ^0e^+\nu `$. The measurement of the branching ratio of $`K^+\pi ^+\nu \overline{\nu }`$ alone will allow one to obtain the magnitude of $`V_{td}`$, the element in the Cabbibo Kobayashi Maskawa matrix, which controls Standard Model CP violation. In conjunction with the neutral channel decay $`K^0\pi ^0\nu \overline{\nu }`$ and measurements in the B sector, this will provide an overconstrained measurement of both the real and the imaginary part of $`V_{td}`$. This will address the question on CP violation: Can it be completely accounted for by the Standard Model. Figure 1 illustrates the proposed sensitivity of various experiments that will measure $`\rho `$ and $`\eta `$, the CP violating parameters of the standard model.

The Standard Model predicts a branching ratio for the decay of $`[0.8\pm 0.3]10^{10}`$ in Buras’s unitarity triangle . One event has been observed by the BNL experiment E787 in a stopped $`K^+`$ experiment. We propose  to measure the decay in flight of $`K^+\pi ^+\nu \overline{\nu }`$ and obtain $`>`$70 events in 2 years of data taking . This will meet the precision suggested by the theory. An intermediate stage, the observation of 5-10 events, is the goal of an approved BNL experiment E949, which is the continuation of E787. Some members of the CKM and E949 are now working in collaboration with each other as two stages of a physics program.

## 2 The Background

One of the challenges of the experiment is to control the background which is 10 orders of magnitude larger than the expected signal.

### 2.1  The $`K^+\pi ^+\pi ^0`$ background

This decay, with branching ratio $``$ 21%, can mimic the signal if the two $`\gamma `$ from the decay of the $`\pi ^0`$ are undetected. It will be controlled by a combination of kinematic rejection and a photon veto system. The CKM detector must veto photons with a very high efficiency.

### 2.2  The $`K^+\mu \nu `$ background

In this decay , with a branching ratio $``$63.5%, the $`\mu `$ can be misidentified as a $`\pi `$ .This is controlled by a muon veto and RICH systems.

### 2.3  Background from Interactions

The interactions of the undecayed beam or daughter $`\pi `$s contribute to the background. In order to minimize this contribution, as little material as possible should be encountered by the particles. This dictates the choice of the various detector elements.

## 3 The Apparatus

In Figure 2 are shown the detector elements which will be detailed below. To detect $`K^+`$ hadronic interactions before the decay volume entrance, a hadron calorimeter surrounding the beam (not described)is used .

### 3.1 The RICH and the Spectrometers

The charged particles in the reaction, namely the K and the $`\pi `$ are both measured by Ring Imaging CHerenkov counters (RICH) as velocity vectors and as momentum vectors in upstream and downstream spectrometers. The gas in the pion RICH is chosen in order to meet the low dispersion required for good resolution. In order to allow as little material in the beam as possible rather thin mirrors and windows are mandatory. This require the RICHs to operate at or near atmospheric pressure. These considerations have contributed to the choice of the K beam to be 22GeV energy and the radiator gas of the decay pion RICH to be neon. The $`\pi `$ momentum resolution achieved by the RICH is better than 1% for this configuration. The K and $`\pi `$ magnetic spectrometers use conventional trackers, silicon strips and straw tubes respectively. Straws have been chosen to minimize the material in the beam. They are 5mm in diameter with a 20 $`\mu `$m sense wire, 5 to 8 straw layers per station (x,y,u,v). $`410^{}4`$ $`X_0`$, with 150$`\mu `$m hit resolution per layer. The straws are located in the vacuum decay volume.

### 3.2 The $`\gamma `$ veto system

As already mentioned one of the main tasks of the experiment is to control the background from the $`K^+\pi ^+\pi ^0`$ decay. This points to the extreme importance of the $`\gamma `$ veto system. The veto system has a component surrounding the decay volume, where the decay takes place and a forward component, after the $`\pi `$ RICH to deal with high energy $`\gamma `$’s. We want to minimize interactions which are an important component of the non-Gaussian tails of the background. Therefore, the component surrounding the decay volume will be located inside the vacuum vessel with a vacuum of $`10^6`$ torr and we call this Vacuum Veto System (VVS). The VVS is made of a sandwich of plastic scintillator and lead 5mm/1mm each, stacked in modules 50 cm long and 40cm radially. There are over 80 such sandwich layers per module and 34 such modules in the VVS. There is a void of 50cm between each module which corresponds to an angle of incidence of $`>40`$ degrees for $`\gamma `$. This angle is not reached by the photons coming from the $`K^+\pi ^+\pi ^0\gamma \gamma `$ decay kinematics . This allows us to reduce the volume and the cost by a factor two.

The required VVS detection inefficiency is given in Figure 3. In the region of medium energies the requirements have been achieved by E787 with a comparable veto system.

For the Forward Veto System (FVS) a choice of various techniques are avalaible, and are still under study. We are assuming, for the moment, a setup of Sc/Pb 5mm/1mm sandwiches for this calorimeter as well.

### 3.3 The Muon Veto

Another important component of the background is the $`K^+\mu ^+\nu `$ which has a branching ratio of $``$63.5%. A muon veto follows the FVS. The non-interacting beam goes through a hole in the FVS and muon veto system and is dumped afterward by a magnetic field into the beam dump.

## 4 The Beam

A very important component of our experiment is the beam. In order to get $``$ 6MHz of $`K^+`$ decays in the fiducial volume we need a $``$30 MHz $`K^+`$ beam with a ratio of $`K/\pi `$ of$`>2`$, whereas nature tends to provide far more $`\pi `$’s than K’s, shown Figure 4 and Figure 5. We will construct a new RF separated beam based upon 3.9 GHz transverse mode superconducting RF cavities, being developed at Fermilab for this purpose. We use a set of 2 RF cavities properly spaced from each other in order for the $`\pi `$s to be deflected by the 1st cavity, and then deflected back by the same amount, by the 2nd cavity. The $`\pi `$s are absorbed by a beam stopper. With the cavities phased this way, protons which arrive at the second cavity just one RF period behind the $`\pi `$s are also undeflected. Ks, however, arrive $`180^{}`$ in RF phase behind the $`\pi `$ and are deflected again with the same sign, causing most of them to miss the beam stopper. A cartoon of the operation is shown in Figure 6.

In Figure 4 we show the behavior of the $`\pi `$ component of the beam which after having been deflected back by the second cavity is stopped. Shown in the same figure are the $`\mu `$s coming from $`\pi `$ decays in the beam line.

In Figure 5 the K component is shown with the decay $`\mu `$s. The Ks which avoid the stopper are shown reaching the detector.

## 5 Status

We are in the process of building prototypes for the photon veto system, the straw tubes and the super-conducting RF cavities. Simulations are being done in parallel for both the straws and the veto systems as well as for the beam line design.

## 6 Conclusion

Figure 7 shows the expected signal and $`K^+\pi ^+\pi ^0`$ background after the photon vetos have been applied from two years of running as a function of the missing mass squared of the neutral system recoiling against the $`\pi ^+`$. 72 events survive in a signal region with 3 unvetoed and mismeasured $`K^+\pi ^+\pi ^0`$ event. The total background estimate for this sample, including $`K^+\mu ^+\nu `$ , interactions and accidentals is $`<8`$ events. Work is now underway to refine and improve our proposal with the intention of seeking approval in the spring 2001.


text:
# Monte Carlo simulations of star clusters – II. Tidally limited, multi–mass systems with stellar evolution.

## 1 INTRODUCTION.

Understanding the dynamical evolution of self–graviting stellar systems is one of the grand challenges of astrophysics. Among its many applications, the one which motivates the work presented here is the study of real globular clusters. Unfortunately, dynamical modelling of globular clusters and other large collisional stellar systems (like galactic nuclei, rich open clusters and galaxy clusters) still suffers from severe drawbacks. These are, on the one hand, due to the very high (and presently unfulfilled) hardware requirements needed to model realistic, large stellar systems by use of the direct $`N`$–body method, and partly due to the poor understanding of the validity of assumptions used in statistical modeling based on the Fokker–Planck and other approximations. The Fokker–Planck method with finite differences, which has recently been greatly improved, can now be used to simulate more realistic stellar systems. It can tackle: anisotropy, rotation, a tidal boundary, tidal shocking by galactic disk and bulge, a mass spectrum, stellar evolution and dynamical and primordial binaries (Takahashi 1995, 1996, 1997, Takahashi & Portegies Zwart 1998, 1999 hereafter TPZ, Drukier et al. 1999, Einsel & Spurzem 1999, Takahashi & Lee 1999 ). Unfortunately, the Fokker–Planck approach suffers, among other things, from the uncertainty of differential cross–sections of many processes which are important during cluster evolution. It can not supply detailed information about the formation and movement of all objects present in clusters. Additionally, the introduction of the mass function, mass loss and tidal boundary into the code is very approximate. Direct $`N`$–body codes are the most natural methods that can be used in simulations of real star clusters (NBODY – Aarseth 1985, 1999a, 1999b, Spurzem & Aarseth 1996, Portegies Zwart et al. 1998 – KIRA). They practically do not suffer from any restrictions common in the Fokker–Planck method. But unfortunately, even special–purpose hardware (Makino et al al. 1997 and references therein) can be used effectively only in $`N`$–body simulations which are limited to a rather unrealistic number of stars (i.e. too small for globular clusters). Another possibility is to use a code which is very fast and provides a clear and unambiguous way of introducing all the physical processes which are important during globular cluster evolution. Monte Carlo codes, which use a statistical method of solving the Fokker–Planck equation, provide all the necessary flexibility. They were developed by Spitzer (1975, and references therein) and Hénon (1975, and references therein) in the early seventies, and substantially improved by Marchant & Shapiro (1980, and references therein) and Stodółkiewicz (1986, and references therein) and recently reintroduced by Giersz (1998, hereafter Paper I, Giersz 2000, see also Giersz in Heggie et al. 1999) and Joshi et al. (1999a), Joshi et al. (1999b, hereafter JNR). The Monte Carlo scheme takes full advantage of the established physical knowledge about the secular evolution of (spherical) star clusters as inferred from continuum model simulations. Additionally, it describes in a proper way the graininess of the gravitational field and the stochasticity of real $`N`$–body systems and provides, in a manner as detailed as in direct $`N`$–body simulations, information about the movement of any objects in the system. This does not include any additional physical approximations or assumptions which are common in Fokker–Planck and gas models (e.g. for conductivity). Because of this, the Monte Carlo scheme can be regarded as a method which lies between direct $`N`$–body and Fokker–Planck models and combines most of their advantages.

Very detailed observations of globular clusters have extended our knowledge about their stellar content, internal dynamics and the influence of the environment on cluster evolution (Janes 1991, Djorgovski & Meylan 1993, Smith & Brodie 1993, Hut & Makino 1996, Meylan & Heggie 1997). They suggest that the galactic environment has a very significant effect on cluster evolution. Gnedin et al (1999) showed that, in addition to the well known importance of gravitational shock heating of the cluster due to passages through the disk or close to the bulge ”shock–induced relaxation” (the second order perturbation) also has a crucial influence on the cluster destruction rate. A wealth of information on ”peculiar” objects in globular clusters (blue stragglers, X–ray sources (high– and low–luminosity), millisecond pulsars, etc.) suggests a very close interplay between stellar evolution, binary evolution and dynamical interactions. This interplay is far from being understood. Moreover, recent observations suggest that the primordial binary fraction in a globular cluster can be as high as $`15\%`$$`38\%`$ (Rubenstein & Bailyn 1997). Monte Carlo codes provide all the necessary flexibility to disentangle the mutual interactions between all physical processes which are important during globular cluster evolution.

The ultimate aim of the project described here is to build a Monte Carlo code which will be able to simulate the evolution of real globular clusters, as closely as possible. In this paper (the second in the series) the Monte Carlo code (which was discussed in detail in Paper I) is extended to include stellar evolution (as described by Chernoff & Weinberg (1990), hereafter CW or Tout et al. 1997), multi–component systems (simulated by continuous, power–law mass function), the tidal field of the Galaxy (simulated by a tidal radius and an appropriate escape criterion) and binary–binary energy generation (as introduced by Stodółkiewicz 1986). This is the first major step in the direction of performing numerical simulations of the dynamical evolution of real globular clusters submerged in the Galaxy potential. The results of simulations will be compared with those of CW, Aarseth & Heggie (1998, hereafter AH) TPZ and JNR.

The plan of the paper is as follows. In Section 2, a short description of the new features introduced into the Monte Carlo code will be presented. In Section 3, the initial conditions will be discussed and results of the simulation will be shown. And finally, in Section 4 the conclusions will be presented.

## 2 MONTE CARLO METHOD.

The Monte Carlo method can be regarded as a statistical way of solving the Fokker–Planck equation. Its implementation presented in Paper I is based on the orbit-averaged Monte Carlo method developed in the early seventies by Hénon (1971) and then substantially improved by Stodółkiewicz (1986, and references therein), Giersz (1998) and recently by JNR. The code is described in detail in Paper I, which deals with simulations of isolated, single–mass systems. In this section, additional physical processes (not included in Paper I) will be discussed: multi–mass systems and stellar evolution (§2.1), mass loss through the tidal boundary (§2.2), formation of three–body binaries and their subsequent interactions with field stars in multi–component systems (§2.3), and interactions between binaries (§2.4). The implementation of an arbitrary mass spectrum in the Monte Carlo method is very straightforward and will not be discussed here (see for example Stodółkiewicz 1982).

### 2.1 Mass Spectrum and Stellar Evolution

Observations give more and more evidence that the mass function in globular/open clusters and for field stars as well is not a simple power–law, but is rather approximated by a composite power–law (e.g. Kroupa, Tout & Gilmore 1993). However, in the present study the simple power–law is used, for simplicity and to facilitate comparison with other numerical simulations (CW, Fukushige & Heggie 1995, Giersz & Heggie 1997, AH, TPZ, JNR). The following form for the initial mass function was assumed:

$$N(m)dm=Cm^\alpha dm,m_{min}mm_{max},$$
(1)
where C and $`\alpha `$ are constants. The detailed discussion of $`\alpha `$, $`m_{min}`$ and $`m_{max}`$ will be presented in §3.1. To describe the mass loss due to stellar evolution the same simplified stellar evolution model as adopted by CW was used. More sophisticated models (Portegies Zwart & Verbunt 1996, Tout et al. 1997) give very similar results for the time which a star spends on the main–sequence, and also slightly smaller masses of remnants of stellar evolution (see Table 1).

This should not significantly change the results of simulations. In CW it is assumed that a star at the end of its main–sequence lifetime instantaneously ejects its envelope and becomes a compact remnant (white dwarf, neutron star or black hole). This is a good approximation, since (from the point of view of cluster evolution) the dominant mass loss phase occurs on a time scale shorter than a few Myrs, and this is much less than the cluster evolution time scale (proportional to the relaxation time), which is of the order 1 Gyr. Mass loss due to stellar winds is neglected for main–sequence stars. According to the prescription given in CW; main–sequence stars of mass $`m>8M_{}`$ finish their evolution as neutron stars of mass $`1.4M_{}`$, while stars of mass $`m<4M_{}`$ end as white dwarfs of mass $`0.58+0.22(m1)`$. Stars of intermediate masses are completely destroyed. For stars with masses smaller than $`0.83M_{}`$ the main–sequence lifetime is linearly (log – log) extrapolated, which is in very good agreement with the scaling $`m^{3.5}`$ used by TPZ and JNR. The initial masses of stars are generated from the continuous distribution given in equation (1). This ensures a natural spread in their lifetimes. To ensures that the cluster remains close to virial equilibrium during rapid mass loss due to stellar evolution, special care is taken that no more than $`3\%`$ of the total cluster mass is lost during one overall time–step.

### 2.2 Tidal Stripping

For a tidally truncated cluster the mass loss from the system is dominated by tidal stripping – diffusion across the tidal boundary. This leads to a much higher rate of mass loss than in an isolated system, where mass loss is attributed mainly to rare strong interactions in the dense, inner part of the system. In the present Monte Carlo code, a mixed criterion is used to identify escapers: a combination of apocenter and energy–based criteria. In the apocenter criterion a star is removed from the system if

$$r_a(E,J)>r_t,$$
(2)
and in the energy criterion a star is removed if

$$E>E_tGM/r_t,$$
(3)
where $`r_a(E,J)`$ is the apocenter distance of a star with energy $`E`$ and angular momentum $`J`$, $`r_t`$ is the tidal radius of the cluster of mass $`M`$, and $`G`$ is the gravitational constant. Recently, TPZ demonstrated that the energy–based criterion can lead to an overestimation of the escape rate from a cluster. Sometimes stars on nearly circular orbits well inside the tidal radius can be removed from the system . The use of the mixed criterion can even further increase the escape rate and shorten the time to cluster disruption. No potential escapers are kept in the system. This mixed criterion was mainly used to facilitate direct comparison with Monte Carlo simulations presented in the collaborative experiment (Heggie et al. 1999). Stars regarded as escapers are lost instantaneously from the system. This is in contrast to $`N`$–body simulations where stars need time proportional to the dynamical time to be removed from the system. Recently, Baumgardt (2000) showed that stars with energy greater than $`E_t`$ in the course of escape can again become bound to the system, because of distant interactions with field stars. This process can substantially influence the escape rate. The mass loss across the tidal boundary can become unstable, when too many stars are removed from the system at the same time. This is characteristic of the final stages of cluster evolution, and for clusters with initially low central concentration. To properly follow these stages of evolution the time–step has to be decreased to force smaller mass loss. JNR used an iteration procedure to determine the mass loss and the tidal radius. A similar procedure will be introduced into a future version of our Monte Carlo code.

### 2.3 Three–Body Binaries

In the present Monte Carlo code (as it was described in Paper I) all stellar objects, including binaries, are treated as single superstars. This allows one to introduce into the code, in a simple and accurate way, processes of stochastic formation of binaries and their subsequent stochastic interactions with field stars and other binaries. The procedure for single–component stellar systems was described in detail in Paper I. Now the procedure for multi–component systems will be discussed.

Suppose that the rate of formation, per unit volume, of three-body binaries with components of mass $`m_1`$ and $`m_2`$, by interaction with stars of mass $`m_3`$ is

$$\frac{dn_b}{dt}=Bn_1n_2n_3,$$
(4)
where $`n_1`$, $`n_2`$ and $`n_3`$ are the number densities of stars of mass $`m_1`$, $`m_2`$ and $`m_3`$, respectively. The coefficient B is a function of the masses of the interacting stars and the local mean kinetic energy (see for example Heggie 1975 and Stodółkiewicz 1986). Suppose also, that the formation of binaries within a radial zone (superzone, see Paper I) of volume $`\mathrm{\Delta }V`$ in a time interval $`\mathrm{\Delta }t`$ is considered. The number of binaries with components of mass $`m_1`$, $`m_2`$, formed by interaction with stars of mass $`m_3`$, in this zone, in this time interval is

$$\mathrm{\Delta }n_b=Bn_1n_2n_3\mathrm{\Delta }V\mathrm{\Delta }t.$$
(5)
Now suppose that all stars in this zone are divided into groups of three, and for each group a binary is created from the first two stars in the group with probability $`P`$. Suppose there are $`\mathrm{\Delta }N`$ shells in the zone. Then, there are $`\mathrm{\Delta }N/3`$ groups of three stars. Also, the probability that, in one given group, the masses of the three stars are $`m_1`$, $`m_2`$ and $`m_3`$ (as required), is $`n_1n_2n_3/n^3`$ (where $`n`$ is the total number density). Therefore, the average number of binaries formed in this zone in this time interval, with components of mass $`m_1`$ and $`m_2`$ (by interaction with a star of mass $`m_3`$) is

$$\mathrm{\Delta }n_b=\frac{\mathrm{\Delta }N}{3}\frac{n_1n_2n_3}{n^3}P.$$
(6)
This is equal to the required number (equation 5) if

$$P=\frac{3Bn^3\mathrm{\Delta }V\mathrm{\Delta }t}{\mathrm{\Delta }N}.$$
(7)
So, to determine the probability of the formation of a three–body binary with components of mass $`m_1`$, $`m_2`$ and $`m_3`$ one has to know only the local total number density instead of the local number density for stars of each mass. This procedure substantially reduces fluctuations in the binary formation rate. The determination of the local number density in the Monte Carlo code is a very delicate matter (Paper I). It is practically impossible if one has to compute the local number densities for each species of a multi–component system.

A binary living in the cluster is influenced by close and wide interactions with field stars. Close interactions are the most important from the point of view of cluster evolution. They generate energy, which supports post–collapse cluster evolution. For a single–component system the probability of a binary field star interaction can be computed using the well known simple semi–empirical formula of Spitzer (Spitzer 1987). For multi–component systems there is no simple semi–empirical formula which can fit all numerical data (Heggie et al. 1996). In the present Monte Carlo model, the following strategy was used to compute the total probability of interaction between a binary of binding energy $`E_b`$ consisting of mass $`m_1`$ and $`m_2`$ and a field star of mass $`m_3`$. The mass dependence of the probability was deduced from equations (6-23), (6-11) and (6-12) presented in Spitzer (1987) and results of Heggie (1975). Then, the coefficient was suitably adjusted, so that the total probability was correct for the equal–mass case. The formula obtained in such a way is as follows

$$P_{3b}=\frac{5\sqrt{2}\pi A_sG^2m_1^2m_2^2\sqrt{m_{123}}n}{8\sqrt{3}\sqrt{m_{12}}\sqrt{m_3}\sqrt{m_a}\sigma E_b}𝑑t,$$
(8)
where $`m_{12}=m_1+m_2`$, $`m_{123}=m_{12}+m_3`$, $`m_a`$ is an average stellar mass in a zone, and $`\sigma `$ is the one dimensional velocity dispersion. This procedure is, of course, oversimplified and in some situations can not give correct results, for example, when a field star is very light compared to the mass of the binary components. The changes of the binding energy of binaries and their velocities due to interactions with field stars are computed in the same way as in Paper I.

### 2.4 Binary–Binary interactions

It is well known that interactions between binaries can play an important role in globular cluster evolution, particularly when primordial binaries are present (Gao et al. 1991, Hut et al. 1992 and reference therein, Meylan & Heggie 1998 and reference therein, Giersz & Spurzem 2000). Binary–binary interactions, besides their role in energy generation, can also be involved in the creation of many different peculiar objects observed in globular clusters. The problem of energy generation in binary–binary interactions is very difficult to solve (see the pioneering work of Mikkola (1983, 1984)). The implementation of interactions between binaries in the Monte Carlo code is based on the method described by Stodółkiewicz (1986). Only strong interactions are considered, and only two types of outcomes are permitted: one binary (composed of the heaviest components of the interacting binaries) and two single stars, or two binaries in a hyperbolic relative orbit. The eccentricities of the orbits of both binaries, and their orientations in space, are neglected. Also stable three–body configurations are not allowed.

In the case when binaries are formed in dynamical processes, and only a few binaries are present at any time in the system, it is very difficult to use the binary density (Giersz & Spurzem 2000) to calculate the probability of a binary–binary interaction. Another approach has to be employed (Stodółkiewicz 1986). For a given binary, the pericenter, $`r_{}`$, and apocenter, $`r_+`$, distances of its orbit in the cluster are known. This binary can hit only binaries (regarded as targets) whose actual distance from the cluster center lie between $`r_{}`$ and $`r_+`$. To compute the probability of this interaction, the following procedure is introduced (a modification of the procedure proposed by Stodółkiewicz 1986). The rate of encounters of binaries within impact parameter $`p`$ is

$$f(\underset{¯}{\text{v}}_1)n_1(\underset{¯}{\text{r}})f(\underset{¯}{\text{v}}_2)n_2(\underset{¯}{\text{r}})|\underset{¯}{\text{v}}_1\underset{¯}{\text{v}}_2|\pi p^2d^3\underset{¯}{\text{v}}_1d^3\underset{¯}{\text{v}}_2d^3\underset{¯}{\text{r}},$$
(9)
where $`f`$ is a distribution function normalized to unity, $`n_1`$ and $`n_2`$ are the number densities of binaries, and $`\underset{¯}{\text{v}}_1`$ and $`\underset{¯}{\text{v}}_2`$ are their barycentric velocities. Suppose there is a binary at a given position $`\underset{¯}{\text{r}}`$ with velocity $`\underset{¯}{\text{v}}^{}`$. The rate of encounters with other binaries at an impact parameter less than $`p`$ is

$$\dot{N}=f(\underset{¯}{\text{v}})n(\underset{¯}{\text{r}})|\underset{¯}{\text{v}}\underset{¯}{\text{v}}^{}|\pi p^2d^3\underset{¯}{\text{v}},$$
(10)
where $`f`$ is normalized to unity and $`n`$ is normalized to the total number of binaries. The number of encounters in time $`\mathrm{\Delta }t`$ is

$$\mathrm{\Delta }N=\mathrm{\Delta }tf(\underset{¯}{\text{v}})n(\underset{¯}{\text{r}})|\underset{¯}{\text{v}}\underset{¯}{\text{v}}^{}|\pi p^2d^3\underset{¯}{\text{v}}.$$
(11)
The number of encounters with one specified binary is

$$\mathrm{\Delta }N_1=\mathrm{\Delta }tf(\underset{¯}{\text{v}})n_s(\underset{¯}{\text{r}})|\underset{¯}{\text{v}}\underset{¯}{\text{v}}^{}|\pi p^2d^3\underset{¯}{\text{v}},$$
(12)
where $`n_s`$ is normalized to unity. Equation 12 can be written in the form

$$<\mathrm{\Delta }N_1>=n_s(\underset{¯}{\text{r}})<|\underset{¯}{\text{v}}\underset{¯}{\text{v}}^{}|>\pi p^2\mathrm{\Delta }t.$$
(13)
Now, $`n_s(\underset{¯}{\text{r}})4\pi r^2`$ is the probability density of $`r`$, i.e.

$$4\pi r^2n_s(\underset{¯}{\text{r}})dr=\{\begin{array}{ccc}\frac{dr/|\dot{r}|}{\underset{r_{}}{\overset{r_+}{}}\frac{dr}{\dot{r}}}& for& r_{}<r<r_+\hfill \\ 0& & otherwise\hfill \end{array}$$
(14)
So, $`4\pi r^2n_s(\underset{¯}{\text{r}})=2/|\dot{r}|/T`$, where $`T`$ is the orbital period of the binary in the cluster. Hence

$$n_s(\underset{¯}{\text{r}})=\frac{1}{2\pi r^2T|\dot{r}|}.$$
(15)
Substituting equation (15) into equation (13) we may obtain a formula for the probability $`P_{3b3b}`$ of the encounter between the two binaries, i.e.

$$P_{3b3b}=\frac{w}{|v_r|}\frac{p^2}{2r^2T}\mathrm{\Delta }t,$$
(16)
where $`w=<|\underset{¯}{\text{v}}\underset{¯}{\text{v}}^{}|>`$, and $`|v_r|=|\dot{r}|`$. It is interesting to note that equation (16) is very similar to equation (36) in Stodółkiewicz (1986), which was obtained in a very approximate way. It differs only by the factor $`w/|v_r|`$, which takes into account in a proper way the geometry of the encounter. To compute the maximum impact parameter, a value equal to $`2.5`$ times the semi–major axis of the softer binary (according to Stodółkiewicz’s (1986) prescription) for the minimum distance during the encounter was adopted. The probability $`P_{3b3b}`$ is evaluated every time step for each pair of binaries (when binaries are sorted against increasing distance from the cluster center). If this probability is smaller than a random number (drawn from a uniform distribution) the binaries are due to interact. The outcome of the interaction is as follows (see details in Stodółkiewicz 1986):

* Two binaries in a hyperbolic relative orbit - $`12\%`$ of all interactions. In this case it is assumed that the recoil energy received by both binaries is equal to $`0.4`$ times the binding energy of the softer binary, and this binding energy increases by the same value.
* One binary and two escapers - $`88\%`$ of all interactions. The total recoil energy released is equal to $`0.516(E_{b1}+E_{b2})`$ and distributed according to conservation of momentum. The softer binary is destroyed and the harder binary increases its binding energy by the amount equal to the recoil energy. $`E_{b1}`$ and $`E_{b2}`$ are the binding energies of interacting binaries.

The new values of the velocity components of the binaries/singles are computed as discussed in Stodółkiewicz (1986) and Giersz (1998).

The procedure described above is very uncertain in regard to the amount of energy generated by binaries in their interactions with field stars and other binaries (particularly for multi–component systems). To solve this problem, it is planned in future to introduce into the Monte Carlo code numerical procedures (based on Aarseth’s NBODY6 code) which can numerically integrate the motion of three– and four–body subsystems. This will at least ensure that the energy generated in these interactions will be calculated properly. A similar procedure was introduced, with success, into the Hybrid code (Giersz & Spurzem 2000). But still there will be remaining uncertainties in the determination of the overall probabilities for binary creation and binary–binary interactions.

## 3 RESULTS

In Paper I the first results of Monte Carlo simulations of the evolution of single–component systems were presented. Here, the Monte Carlo code is extended to include a power–law mass function, stellar evolution, tidal stripping and binary–binary interactions.

### 3.1 Initial Models

The initial conditions were chosen in a similar way as in the “collaborative experiment” (Heggie et al. 1999). The positions and velocities of all stars were drawn from a King model with a power–law mass spectrum. All standard models have the same total mass $`M=60000M_{}`$ and the same tidal radius $`R_c=30`$ pc. Masses are drawn from the power–law mass function according to equation (1). The minimum mass was chosen to be $`0.1M_{}`$ and the maximum mass be $`15M_{}`$. Three different values of the power–law index were chosen: $`\alpha =1.5`$, $`2.35`$ and $`3.5`$. The sets of initial King models were characterized by $`W_0=3`$, $`5`$ and $`7`$. To facilitate comparison with results of CW, AH, TPZ and JNR additional models of CW’s Family 1 were computed (minimum mass equal to $`0.4M_{}`$, $`\alpha =1.5`$, $`2.5`$, $`3.5`$ and total mass $`M=90685M_{}`$, $`99100M_{}`$, $`103040M_{}`$, respectively). In this paper the results of the standard models will mainly be presented. Models of Family 1 will not be shown (with one exception - Figure 12) and will only be discussed in cases, where their results are different from these of the standard models. All models are listed in Table 2.

The initial model is not in virial equilibrium, because of statistical noise and because masses are assigned independently from the positions and velocities. Therefore it has to be initially rescaled to virial equilibrium. During all the simulations the virial ratio is kept within $`<2\%`$ of the equilibrium value (for the worst case within $`<5\%`$). Standard $`N`$–body units (Heggie & Mathieu 1986), in which the total mass $`M=1`$, $`G=1`$ and the initial total energy of the cluster is equal to $`1/4`$, have been adopted for all runs. Monte Carlo time is equal to $`N`$–body time divided by $`N/ln(\gamma N)`$, where gamma was adopted to be $`0.11`$ as for single–component systems (Giersz & Heggie 1994). However, there are some results suggesting even smaller values of $`\gamma `$ ($`0.015`$) for multi–components systems (Giersz & Heggie 1996, 1997). In order to scale time to the physical units the following formula is used

$$\frac{t_{scale}}{10^6yrs}=14.91\left(\frac{R_c}{r_{t_M}}\right)^{1.5}\frac{N_T}{\sqrt{M}ln(\gamma N_T)},$$
(17)
where $`r_{t_M}`$ is the tidal radius in Monte Carlo units, $`M`$ and $`R_c`$ are in cgs units. $`r_{t_M}`$ and $`t_{scale}`$ are listed in Table 2. In the course of evolution, when the cluster loses mass, the tidal radius changes according to $`r_tM^{1/3}`$.

Finally, a few words about the efficiency of the Monte Carlo code are presented here. The simulations of the largest models (consist of $`360195`$ stars) took about one month on a Pentium II 300 MHz PC. This is still much shorter than the biggest direct $`N`$-body simulations performed on GRAPE–4 (Teraflop special–purpose hardware). Nevertheless, the speed of the code is not high enough to perform, in a reasonable time, a survey of models, as was done, for example, by CW. It is clear, that to simulate evolution of real globular clusters a substantial speed up of the code is needed. This can be done, either by parallelizing the code (in a similar way to JNR), introducing a more efficient way of determining the new positions of superstars, or by using a hybrid code (as was done by Giersz & Spurzem 2000). However, the Monte Carlo code has presently a great relative advantage over the $`N`$-body code for simulations with a large number of primordial binaries. Primordial binaries substantially downgrade the performance of $`N`$-body codes on supercomputers or on special-purpose hardware. They can be introduced into the Monte Carlo code in a natural way, practically without a substantial loss of performance (Giersz & Spurzem 2000).

### 3.2 Global evolution.

To facilitate comparison of standard models with recently obtained results of $`N`$–body (AH), Fokker–Planck (CW, TPZ) and Monte Carlo (JNR) simulations of globular cluster evolution, the following parameter (introduced by CW), which describes the initial relaxation time, is calculated. It is defined by

$$F\frac{M}{M_{}}\frac{R_g}{kpc}\frac{220kms^1}{v_g}\frac{1}{lnN},$$
(18)
where $`R_g`$ is the distance of the globular cluster to the Galactic center, and $`v_g`$ is the circular speed of the cluster around the Galaxy. For the initial parameters of the standard models ($`M=60000M_{}`$, $`r_t=30`$ $`pc`$, $`v_g=220`$ $`kms^1`$), $`R_g=3.89`$ $`kpc`$. Therefore according to equation 18, $`F`$ is equal to $`2.14`$x$`10^4`$, $`1.92`$x$`10^4`$ and $`1.82`$x$`10^4`$ for $`\alpha `$ equal to $`1.5`$, $`2.35`$ and $`3.5`$, respectively. This is about $`2.5`$ times smaller than for Family 1 of CW. Generally, the greater the value of $`F`$, the longer the relaxation time and the slower the evolution. Therefore one should expect, that for standard models the collapse or disruption times are shorter than for CW, AH, TPZ or JNR. Additionally, because the minimum stellar mass is $`0.1M_{}`$ instead of $`0.4M_{}`$ (the value adopted by CW, AH, TPZ and JNR), one should expect that mass segregation will proceed faster in standard models, and the collapse time should be further shortened. Indeed, as can be seen in Table 3, this is true with the exception of models with a flat mass function (W3235 and probably W715). For these models the dominant physical process during most of the cluster life is stellar evolution. The standard models contain a smaller number of massive main–sequence stars than Family 1 models. Therefore they should evolve more slowly.

The standard models show very good agreement with $`N`$body results (Heggie 2000). See columns labeled by G and H-0.1 in Table 3. Only models with a flat mass function show substantial disagreement. These models are difficult for both methods. Violent stellar evolution and induced strong tidal stripping lead to troubles with time–scaling for the $`N`$body model and with the proper determination of the tidal radius for the Monte Carlo model. Generally, the same is true for Monte Carlo models of Family 1. Results of these models show good agreement with the results of CW, AH and TPZ. JNR’s results, particularly for strongly concentrated systems, disagree with all other models. This may be connected with the fact that JNR’s Monte Carlo scheme is not particularly suitable for high central density and strong density contrast. The deflection angles adopted by JNR are too large and consequently the time-steps are too large, which can lead to too fast evolution for these models.

Figures 1 to 3 present the evolution of the total mass for standard models.

Three phases of evolution are clearly visible. The first phase is connected with the violent mass loss due to evolution of the most massive stars. Then there is the long phase connected with gradual mass loss due to tidal stripping. And finally there is the phase connected with the tidal disruption of a cluster. This phase is not well represented by most models presented here. The chosen overall time–step is too long to follow in detail the evolution, which proceeds basically on a dynamical time scale. Only for models W515 and W715 was a sufficiently small time–step adopted to properly follow this phase of evolution. Models with a steep mass–function ($`\alpha 2.35`$) and with different $`W_o`$ show very similar evolution during the phase of tidal stripping (except W3235). It seems that the initial mass loss due to stellar evolution is not sufficiently strong to substantially change the initial cluster structure. Models of Family 1 do not show this feature. It seems that the initial mass loss across the tidal boundary is sufficiently strong to change the structure of the system. These models contain initially a much larger number of massive main–sequence stars than standard models. Figures 1 to 3 present qualitatively very similar features as figures shown by TPZ and JNR.

Figures 4 to 6 show the amount of mass loss due to stellar evolution, $`evo`$, and tidal stripping, $`esc`$, for standard models.

The final ratio of the amount of mass loss due to tidal stripping to the amount of mass loss due to stellar evolution is smallest for the shallow mass function ($`\alpha =1.5`$). The steeper the mass function, the higher the ratio. Mass loss connected with stellar evolution dominates the initial phase of cluster evolution, as could be expected. Then the rate of stellar evolution substantially slows down and escape due to tidal stripping takes over. During this phase of evolution the rate of mass loss is nearly constant, and higher for shallower mass functions. Energy carried away by stellar evolution events dominates the energy loss due to tidal stripping, even though the tidal mass loss is higher (see Figure 7). The higher the initial cluster concentration, the smaller the absolute value of energy loss due to tidal stripping. For a more concentrated cluster the overall mass loss due to tidal stripping is smaller, and so the energy loss is smaller. It is interesting to note, that for model W735 the energy loss due to tidal stripping is positive (see Figure 7) during the second phase of evolution. This means that stars with positive energy are preferentially lost from the cluster, as in isolated systems. A tidally limited system behaves like an isolated one, where strong binary–single and binary–binary interactions are responsible for escapers with a large positive energy. These relatively small numbers of escapers energetically dominate the numerous escapers (due to tidal stripping) with small negative energy.

The typical evolution of the central potential and Lagrangian radii for the standard models is shown in Figures 8 and 9, respectively. For clusters which are not disrupted due to strong mass loss connected with stellar evolution of the most massive stars, the three different phases of evolution can be clearly distinguished. The first phase of violent mass loss due to stellar evolution leads to overall cluster expansion and an increase of the central potential. This phase is less pronounced in Figure 9, because it is superposed on the shrinking of the tidal radius (Lagrangian radii being calculated for smaller total mass). Because of mass segregation the $`1\%`$ Lagrangian radius is decreasing. The second phase is characterized by a decrease of the central potential (except for models with low concentration and shallow mass function, for which strong tidal stripping continuously increases the central potential – W325). The cluster undergoes core collapse and behaves as an ordinary isolated system. The higher the cluster concentration, the shorter this phase is. Then, in the third phase, post–collapse evolution is superposed on effects of tidal stripping. The central potential (on average) is continuously increasing. The cluster contracts nearly homogeneously. The central parts of the system show signs of gravothermal oscillations. Family 1 models show the same features as discussed above.

Generally, the agreement throughout the evolution between the results presented here and these by AH for $`N`$body models, TPZ for 2-D Fokker–Planck models and by JNR for Monte Carlo models is rather good. In all cases, the qualitative behaviour is identical, even though the models of AH, TPZ and JNR have longer lifetimes than the standard models.

### 3.3 Anisotropy evolution.

The degree of velocity anisotropy is measured by a quantity

$$\beta 2\frac{2\sigma _t^2}{\sigma _r^2},$$
(19)
where $`\sigma _r`$ and $`\sigma _t`$ are the radial and tangential one–dimensional velocity dispersions, respectively. For isotropic systems $`\beta `$ is equal to zero. Systems preferentially populated by radial orbits have $`\beta `$ positive and systems preferentially populated by tangential orbits have $`\beta `$ negative. The evolution of anisotropy is presented in Figures 10 to 12 for models W515, W725 and W735-4, respectively. As can be seen in Figure 10, for systems which can not survive the violent initial mass loss, the anisotropy stays close to zero. The system does not live long enough to develop strong positive or negative anisotropy. These clusters which undergo core collapse, are strongly concentrated and have a steep mass function (W535, W7235), develop a small positive anisotropy in the outer and middle parts of the system (see Figure 11). The amount of anisotropy in the outer parts of the system is reduced substantially by tidal stripping, as stars on radial orbits escape preferentially. As tidal stripping exposes deeper and deeper parts of the system, the anisotropy (for large Lagrangian radii) gradually decreases and eventually becomes slightly negative. Most stars on radial orbits were removed from the systems and most stars on tangential orbits remained. At the same time the anisotropy in the middle and inner parts of the system stays close to zero. Finally, just before cluster disruption, the anisotropy in the whole system becomes slightly positive again. For model W735 the evolution of anisotropy proceeds in a very similar way (at least long before the cluster disruption) as for isolated clusters. Positive anisotropy develops throughout most of the system. In the middle parts of the system it becomes nearly constant just after core bounce (post–collapse evolution). Only in the outer parts of the system does the anisotropy increase, showing the development of the cluster halo. When tidal stripping becomes important anisotropy will be reduced and it will behave as in the model W725 discussed above. For models less strongly concentrated and with a flatter mass function (W3235, W5235) the anisotropy in the outer parts of the system very quickly becomes negative. Strong stellar evolution and induced tidal stripping preferentially force these stars to stay in the system, which are on tangential orbits. Just before cluster disruption the anisotropy in the whole system again becomes slightly positive.

For all models of Family 1 (see, for example, Figure 12) the evolution of anisotropy for the outer parts of the system proceeds in a similar way as for standard models W3235 and W5235. Family 1 models contain more massive stars than the standard models. Therefore the stronger mass loss very quickly forces the anisotropy in the outer parts of the system to become negative. It stays negative until cluster disruption, when it becomes slightly positive. The anisotropy in the central parts of the system stays close to zero.

The anisotropy of the main–sequence stars shows very similar behaviour to that discussed above. Mass segregation forces white dwarfs and neutron stars to occupy the inner parts of the system preferentially. Therefore the anisotropy for them is much more modest than for the main–sequence stars, and stays close to zero. The evolution of the anisotropy agrees qualitatively with the results obtained by Takahashi (1997).

### 3.4 Mass segregation.

Three main effects may be expected to govern the evolution of the mass function in models of the kind studied in this paper. First, there is a period of violent mass loss due to evolution of the most massive stars. It takes place mainly during the first few hundred million years and its amplitude strongly depends on the slope of the mass function. Secondly, there is the process of rapid mass segregation, caused by two–body distant encounters (relaxation). It takes place mainly during core collapse and, as can be expected, is relatively unaffected by the presence of a tidal field. Thirdly, there is the effect of the tidal field itself, which becomes important after core collapse or even earlier (for models with a shallow mass function and with low concentration). Because the tides preferentially remove stars from the outer parts of the system, which (because of mass segregation) are mainly populated by low–mass stars, one can expect that the mean mass should increase as evolution proceeds, making allowance for stellar evolution.

The basic results are illustrated in Figures 13 to 15 for the total average mass inside the $`10\%`$ and $`50\%`$ Lagrangian radii and the tidal radius, and in Figures 16 to 18 for the overall average mass of main–sequence stars and white dwarfs. The violent and strong mass loss due to stellar evolution is characterized by an initial decrease of the average mass. This is best seen in Figure 13 for model W715. The initial average mass in this model is much larger than for other models. It contains a relatively small number of low mass stars and a large number of massive stars. Therefore the evolution of the most massive stars will remove a substantial amount of mass from the system and, consequently, greatly lower the average mass in the whole system. This behaviour is visible (with smaller amplitude) in other models (with steeper mass function) as well, but with one exception. The average mass inside the $`10\%`$ Lagrangian radius is increasing instead of decreasing (see Figures 14 and 15). For these models, there is enough time for mass segregation to force the most massive stars (in this situation less massive main–sequence stars, neutron stars and massive white dwarfs) to sink into the center and increase the average mass there.

Family 1 models show a much stronger initial decrease of the average mass than the standard models. This is connected with the fact that they contain a much larger number of massive stars than the standard models. The initial mass segregation is only visible for model W735-4 for the 10% Lagrangian radius.

For the standard models with $`\alpha =2.35`$ and $`3.5`$, mass segregation substantially slows down in the inner parts of the system after the core collapse. This is in good agreement with results obtained by Giersz & Heggie (1997) for small $`N`$-body simulations, but the reason for that behaviour is still unclear. During the post–collapse phase, for systems with a steeper mass function, mass segregation proceeds further, but at a smaller rate. The effect of tides manifests itself by a gradual increase of the average mass. This increase, as expected, is fastest in the halo and slowest in the core. There is one exception to this rule models which are disrupted before core collapse show a nearly constant rate of increase of the average mass inside the whole system (see Figures 13 and 14). This is probably connected with the fact, that for these disrupting systems stars are removed from the whole body of the cluster.

Family 1 models which enter the long post–collapse evolution phase show a steady, slow decrease of the average mass in the whole system. The total initial average mass for these models is much larger than for the standard models, and therefore stellar evolution is important for a longer time and is competitive with tidal stripping.

The evolution of the overall average mass for main–sequence stars and white dwarfs in the standard models shows basically the same features as discussed above (see Figures 16 to 18). When less and less massive main–sequence stars finish their evolution as less and less massive white dwarfs, the average mass of white dwarfs decreases with time. It seems, that in nearly the whole first Gyr the rate of decrease of white dwarf average mass is practically independent of the central concentration and the initial mass function. The differences start to build up later and become most evident close to the collapse time. This can be explained by the fact that models with a steeper mass function contain a smaller number of very massive main–sequence stars. Therefore the total mass of white dwarfs is smaller for these models than for models with a shallower mass function. When less and less massive main–sequence stars finish their evolution, the newly created less massive white dwarfs affect the average mass of white dwarfs more strongly for the steeper mass function than for the shallower one. At the time around core bounce the average mass of white dwarfs increases, particularly for low concentration models with a shallow mass function. Binaries start to be created mainly from the most massive stars (neutron stars and white dwarfs) deep in the core. In interacting with field stars these binaries remove them from the system – preferentially less massive white dwarfs. At the time of core bounce less concentrated models with a shallower mass function are on the verge of disruption. They contain only a small number of stars, and therefore removal of some low mass white dwarfs can lead to substantial changes of the average mass.

For models which enter post–collapse evolution the average mass of main–sequence stars behaves as expected. There is a small initial decrease of the average mass, connected with stellar evolution of the most massive stars. Then there is a period of the gradual increase of the average mass due to tidal stripping and the preferential removal of the least massive stars. At late phases of cluster evolution the increase of average mass speeds up. A different behaviour of the average mass of main–sequence stars is observed for clusters which are disrupted just before core collapse (W515 and W715). For these models removal of low mass main–sequence stars is so strong (because of violent mass loss due to stellar evolution of the most massive stars and induced tidal stripping) that the average mass increases slightly (see Figure 17 and 18). Then the average mass quickly decreases because of stellar evolution and nearly homologous mass removal from the system – most stars are on radial orbits (see Figures 12 and 13).

Generally, Family 1 models show the same features as discussed above. Only for models which enter the long post–collapse evolution phase does the average main–sequence mass decrease slightly instead of increasing. This is in agreement with the evolution of the total average mass for Family 1 models discussed above.

## 4 CONCLUSIONS.

This paper is continuation of Paper I, in which it was shown that the Monte Carlo method is a robust scheme to study, in an effective way, the evolution of a large $`N`$–body systems. The Monte Carlo method describes in a proper way the graininess of the gravitational field and the stochasticity of the real $`N`$–body systems. It provides, in almost as much detail as $`N`$–body simulations, information about the movement of any object in the system. In that respect the Monte Carlo scheme can be regarded as a method which lies between direct $`N`$–body and Fokker–Planck models and combines most advantages of both these methods. The code has been extended to include stellar evolution, multi–component systems described by a power–law mass function, the tidal field of a parent galaxy and generation of energy by binary–binary interactions. This is the first major step in the direction of simulating the evolution of a real globular cluster.

A small survey (similar to these presented by CW, AH, TPZ and JNR) on the evolution of globular clusters in the Galactic tidal field was carried out. It was shown that the results obtained are in qualitative agreement with these presented by CW, AH, TPZ and JNR. Particularly good agreement is obtained with AH’s $`N`$-body simulations. JNR’s Monte Carlo results, mainly for strongly concentrated models, disagree with the other models. The Monte Carlo scheme discussed in this paper has some problems with low concentration models and these with a flat mass function. It seems that the collapse times for these models are too short, in comparison with the results of the other methods. However, these models are problematic for all methods; for example, the $`N`$–body method has problems with proper time scaling. All standard models, for which mass loss due to violent stellar evolution of the most massive stars does not induce quick cluster disruption, evolve in a very similar way. The rate of mass loss, the evolution of the central potential, and the evolution of the average mass, do not depend much on the central concentration of the system. They depend strongly on the index of the mass function. Models of Family 1, on the contrary, show dependence on the initial concentration, as well. The very high initial mass loss across the tidal boundary, connected with the evolution of the most massive stars (for models of Family 1, there are more massive stars than for standard models) forces to substantial changes in the structure of the system and in consequence different evolution of the total mass, anisotropy, etc. Models which are quickly disrupted show only small signs of mass segregation. Models with larger central concentration survive the phase of rapid mass loss and then undergo core collapse and subsequent post–collapse expansion in a manner similar to isolated models. The expansion phase is eventually reversed when tidal limitation becomes important. As in isolated models, mass segregation substantially slows down by the end of core collapse. After core bounce there is a substantial increase in the mean mass in the middle and outer parts of the system, caused by the preferential escape of stars of low mass and by tidal effects. Standard models, which are not quickly disrupted, show a modest initial build up of anisotropy in the outer parts of the system. As tidal stripping exposes the inner parts of the system, the anisotropy gradually decreases and eventually becomes slightly negative. The central part of the system stays nearly isotropic. From the very beginning models of Family 1 develop a modest negative anisotropy in the outer parts of the system. It stays negative until the time of cluster disruption, when it becomes slightly positive (during cluster disruption most stars are on radial orbits).

In order to perform simulations of real globular clusters several additional physical effects have to be included into the code. The tidal shock heating of the cluster due to passages through the Galactic disk, interaction with the bulge, shock–induced relaxation, primordial binaries, physical collisions between single stars and binaries are some of them. Inclusion of all these processes does not pose a fundamental theoretical or technical challenge. It will allow us to perform detailed comparison between simulations and the observed properties of globular clusters, and will also help to understand the conditions of globular cluster formation and explain how peculiar objects observed in clusters can be formed. These kinds of simulations will also help us to introduce, in a proper way, into future $`N`$–body simulations all the necessary processes required to simulate evolution of real globular clusters on a star–by–star basis from their birth to their death.

Acknowledgments I would like to thank Douglas C. Heggie and Rainer Spurzem for stimulating discussions, comments and suggestions. I also thank DCH for making the $`N`$–body results of standard model simulations available. This work was partly supported by the Polish National Committee for Scientific Research under grant 2–P03D–022–12.


text:
# References

An Explanation on Negative Mass-Square of Neutrinos

Tsao Chang

Center for Space Plasma and Aeronomy Research

University of Alabama in Huntsville

Huntsville, AL 35899

Email: changt@cspar.uah.edu

Guangjiong Ni

Department of Physics, Fudan University

Shanghai, 200433, China

Abstract: It has been known for many years that the measured mass-square of neutrino is probably negative. For solving this puzzle, we have further investigated the hypothesis that neutrinos are superluminal fermions. A new Dirac-type equation is proposed and a tachyonic quantum theory is briefly discussed. This equation is equivalent to two Weyl equations coupled together via nonzero mass while respecting the maximum parity violation, and it reduces to one Weyl equation when the neutrino mass becomes zero.

PACS number: 14.60.Lm, 14.60.Pq, 14.60.St

1. Introduction

The square of the neutrino mass is measured in tritium beta decay experiments by fitting the shape of the beta spectrum near endpoint. In many experiments, it has been found to be negative. Most recent data are listed in ”Review of Particle Physics, 2000”  and references therein. The weighted average from two experiments reported in 1999 \[2-3\] is

$$m^2(\nu _e)=2.5\pm 3.3eV^2$$
$`(1)`$
However, other nine measurements from different experiments in 1991-1995 are not used for averages. For instance, a value of $`m^2`$ ($`\nu _e`$) = -130 $`\pm `$ 20 $`eV^2`$ with 95$`\%`$ confidence level was measured in LLNL in 1995.

Furthermore, the pion decay experiment also obtained a negative value for $`\mu `$-neutrinos .

$$m^2(\nu _\mu )=0.016\pm 0.023MeV^2$$
$`(2)`$
The negative value of the neutrino mass-square simply means:

$$E^2/c^2p^2=m^2(\nu _e)c^2<0$$
$`(3)`$
The right-hand side in Eq. (3) can be rewritten as (- $`m_s^2c^2`$), then $`m_s`$ has a positive value. Eq. (1) and (2) suggests that neutrinos might be particles faster than light, no matter how small the $`m_s`$ is. This possibility is further investigated in this paper.

Based on special relativity and known as re-interpretation rule, superluminal particles were proposed by Bilaniuk et al in the Sixties \[6-8\]. The sign of 4-D world line element, $`ds^2`$, is associated with three classes of particles. For simplicity, let $`dy=dz=0`$, then

$$>0\mathrm{C}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{I}(\mathrm{subluminal}\mathrm{particles})$$

$$ds^2=c^2dt^2dx^2=0\mathrm{C}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{I}\mathrm{I}(\mathrm{photon})$$
$`(4)`$
$$<0\mathrm{C}\mathrm{l}\mathrm{a}\mathrm{s}\mathrm{s}\mathrm{I}\mathrm{I}\mathrm{I}(\mathrm{superluminal}\mathrm{particles})$$

For Class III particles, i.e. superluminal particle, the relation of momentum and energy is shown in Eq.(3). The negative value on the right-hand side of Eq. (3) for superluminal particles means that $`p^2`$ is greater than $`(E/c)^2`$. The velocity of a superluminal particle, $`u_s`$, is greater than the speed of light. The momentum and energy in terms of $`u_s`$ are as follows:

$$p=\frac{m_su_s}{\sqrt{u_s^2/c^21}},E=\frac{m_sc^2}{\sqrt{u_s^2/c^21}}$$
$`(5)`$
where the subscript $`s`$ means superluminal particle, i.e. tachyon. From Eq. (5), it is easily seen that when $`u_s`$ is increased, both of $`p`$ and $`E`$ would be decreased. This property is opposite to Class I particle.

Any physical reference system is built by Class I particles (atoms, molecules etc.), which requires that any reference frame must move slower than light. On the other hand, once a superluminal particle is created in an interaction, its speed is always greater than the speed of light. Neutrino is the most possible candidate for a superluminal particle because it has left-handed spin in any reference frame. On the other hand, anti-neutrino always has right-handed spin.

The first step in this direction is usually to introduce an imaginary mass, but these efforts could not reach a point for constructing a consistent quantum theory. Some early investigations of a Dirac-type equation for tachyonic fermions can be found in Ref.. An alternative approach was investigated by Chodos et al. . They examined the possibility that neutrinos might be tachyonic fermions. A form of the lagrangian density for tachyonic neutrinos was proposed. Although they did not obtain a satisfatory quantum theory for tachyonic fermions, they suggested that more theoretical work would be needed to determine physically acceptable tachyonic theory.

2. A new Dirac-type equation

In this paper, we will start with a different approach to derive a new Dirac-type equation for tachyonic neutrinos. In order to avoid introducing imaginary mass, Eq. (3) can be rewritten as

$$E=(c^2p^2m_s^2c^4)^{1/2}$$
$`(6)`$
where $`m_s`$ is called proper mass, for instance, $`m_s(\nu _e)`$ = 1.6 eV from Eq. (1). We follow Dirac’s search , Hamiltonian must be first order in momentum operator $`\widehat{p}`$ :

$$\widehat{E}=c(\stackrel{}{\alpha }\widehat{p})+\beta _sm_sc^2$$
$`(7)`$
with ($`\widehat{E}=i\mathrm{}/t,\widehat{p}=i\mathrm{}`$), where $`\stackrel{}{\alpha }=(\alpha _1,\alpha _2,\alpha _3)`$ and $`\beta _s`$ are 4$`\times `$ 4 matrix, which are defined as

$$\alpha _i=\left(\begin{array}{cc}0& \sigma _i\\ \sigma _i& 0\end{array}\right),\beta _s=\left(\begin{array}{cc}0& I\\ I& 0\end{array}\right)$$
$`(8)`$
where $`\sigma _i`$ is 2$`\times `$2 Pauli matrix, $`I`$ is 2$`\times `$2 unit matrix. Notice that $`\beta _s`$ is a new matrix, which is different from the one in the traditional Dirac equation. the relation between the matrix $`\beta _s`$ and the traditional matrix $`\beta `$ is as follows:

$$\beta _s=\beta \gamma _5$$
$`(9)`$
where

$$\beta =\left(\begin{array}{cc}I& 0\\ 0& I\end{array}\right),\gamma _5=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$$
$`(9a)`$
When we take square for both sides in Eq.(7), and consider the following relations:

$$\alpha _i\alpha _j+\alpha _j\alpha _i=2\delta _{ij}$$

$$\alpha _i\beta _s+\beta _s\alpha _i=0$$

$$\beta _s^2=1$$
$`(10)`$
the relation in Eq. (3) or Eq. (6) is reproduced. Since Eq. (6) is related to Eq. (5), this means $`\beta _s`$ is a right choice to describe neutrinos as superluminal particles.

Denote the wave function as

$$\mathrm{\Psi }=\left(\begin{array}{c}\phi (\stackrel{}{x},t)\\ \chi (\stackrel{}{x},t)\end{array}\right)with\phi =\left(\begin{array}{c}\phi _1\\ \phi _2\end{array}\right),\chi =\left(\begin{array}{c}\chi _1\\ \chi _2\end{array}\right)$$
$`(11)`$
From Eq. (7), the complete form of the new Dirac-type equation becomes

$$\widehat{E}\mathrm{\Psi }=c(\stackrel{}{\alpha }\widehat{p})\mathrm{\Psi }+\beta _sm_sc^2\mathrm{\Psi }$$
$`(7a)`$
Since $`\beta _s^2=1`$ in the Eq.(7a), the new Dirac-type equation is different from the traditional Dirac equation in any covariant representation in terms of the $`\gamma `$ matrices.

We now study the spin-1/2 property of neutrino as a tachyonic Fermion. Eq. (7a), can be rewritten as a pair of two-component equations:

$$i\mathrm{}\frac{\phi }{t}=ic\mathrm{}\stackrel{}{\sigma }\chi +m_sc^2\chi $$

$$i\mathrm{}\frac{\chi }{t}=ic\mathrm{}\stackrel{}{\sigma }\phi m_sc^2\phi $$
$`(12)`$
Eq. (12) is an invariant under the space-time inversion transformation with ($`\stackrel{}{x}\stackrel{}{x},tt`$) ,

$$\phi (\stackrel{}{x},t)\chi (\stackrel{}{x},t),\chi (\stackrel{}{x},t)\phi (\stackrel{}{x},t)$$
$`(13)`$
It means that Eq.(12) obeys CPT theorem. From the equation (12), the continuity equation can be derived:

$$\frac{\rho }{t}+\stackrel{}{j}=0$$
$`(14)`$
and we have

$$\rho =\phi ^{}\chi +\chi ^{}\phi ,\stackrel{}{j}=c(\phi ^{}\stackrel{}{\sigma }\phi +\chi ^{}\stackrel{}{\sigma }\chi )$$
$`(15)`$
where $`\rho `$ and $`\stackrel{}{j}`$ are the probability density and current; $`\phi ^{}`$ and $`\chi ^{}`$ are the Hermitian adjoint of $`\phi `$ and $`\chi `$ respectively.

Eq.(15) can be rewritten as

$$\rho =\mathrm{\Psi }^{}\gamma _5\mathrm{\Psi },\stackrel{}{j}=c(\mathrm{\Psi }^{}\gamma _5\stackrel{}{\alpha }\mathrm{\Psi })$$
$`(15a)`$
It is easy to see that the probability density $`\rho `$ is positive definite when the components in $`\phi `$ and $`\chi `$ are positive.

Considering a plane wave along the $`z`$ axis for a left-handed particle $`(\stackrel{}{\sigma }\stackrel{}{p})/p=1`$, the equations (12) yields the following solution:

$$\chi =\frac{cpm_sc^2}{E}\phi $$
$`(16)`$
We now consider a linear combination of $`\phi `$ and $`\chi `$,

$$\xi =\frac{1}{\sqrt{2}}(\phi +\chi ),\eta =\frac{1}{\sqrt{2}}(\phi \chi )$$
$`(17)`$
where $`\xi (\stackrel{}{x},t)`$ and $`\eta (\stackrel{}{x},t)`$ are two-component spinor functions. In terms of $`\xi `$ and $`\eta `$, Eq. (15) becomes

$$\rho =\xi ^{}\xi \eta ^{}\eta ,\stackrel{}{j}=c(\xi ^{}\stackrel{}{\sigma }\xi +\eta ^{}\stackrel{}{\sigma }\eta )$$
$`(18)`$
In terms of Eq. (17), the equation (12) can be rewritten in the Weyl representation:

$$i\mathrm{}\frac{\xi }{t}=ic\mathrm{}\stackrel{}{\sigma }\xi m_sc^2\eta $$

$$i\mathrm{}\frac{\eta }{t}=ic\mathrm{}\stackrel{}{\sigma }\eta +m_sc^2\xi $$
$`(19)`$
In the above equations, both $`\xi `$ and $`\eta `$ are coupled via nonzero $`m_s`$.

In order to compare Eq. (19) with the well known two-component Weyl equation, we take a limit $`m_s=0`$, then the first equation in Eq.(19) reduces to

$$\frac{\xi _\nu }{t}=c\stackrel{}{\sigma }\xi _\nu $$
$`(20)`$
and the second equation in Eq. (19) vanishes because $`\phi =\chi `$ for the massless limit.

Eq. (20) is the two-component Weyl equation for describing neutrinos, which is related to the maximum parity violation discovered in 1956 by Lee and Yang . They pointed out that no experiment had shown parity to be good symmetry for weak interaction. Now we see that, in terms of Eq.(19), once if neutrino has some mass, no matter how small it is, two equations are coupled together via the mass term while still respecting the maximum parity violation.

3. Remarks

In this paper, we have further investigated the hypothesis that neutrinos are tachyonic fermions. A new Dirac-type equation is proposed and a tachyonic quantum theory is briefly discussed, which can be used to solve the puzzle of negative mass-square of neutrinos.

Notice that the matrix $`\beta _s`$ in the new Dirac-type equation is not a $`4\times 4`$ hermitian matrix. However, based on the above study, we now realize that the violation of hermitian property is related to the violation of parity. Though a non- hermitian Hamiltoian is not allowed for a subluminal particle, it does work for superluminal neutrinos.

In the light of the above theoretical framework, the consequences derived from the new Dirac-type equation agree with all known properties of neutrinos. Therefore, we have reached a conclusion that neutrinos are tachyonic fermions with permanent helicity if the neutrino mass-square is negative. Therefore, more accurate tritium beta decay experiments are needed to further determine the neutrino mass-square.

According to special relativity , if there is a superluminal particle, it might travel backward in time. However, a re-interpretation rule has been introduced since the Sixties \[6-8\]. Another approach is to introduce a kinematic time under a non-standard form of the Lorentz transformation \[17-20\]. Therefore, special relativity can be extended to space-like region, and superluminal particles are allowed without causality violation.

We wish to thank S.Y. Zhu and Y. Takahashi for helpful discussions