text:
# The GRID and the Linux Farm at the RCF

## I BACKGROUND

The RHIC Computing Facility (RCF) is a large scale data processing facility at Brookhaven National Laboratory (BNL) for the Relativistic Heavy Ion Collider (RHIC), a collider dedicated to high-energy nuclear physics experiments.

RHIC’s first physics collisions occurred in the Summer of 2000, when all four experiments began recording data from the collisions. Year 3 of RHIC operations is currently underway.

The RCF provides for the computational needs of the five RHIC experiments (BRAHMS, PHENIX, PHOBOS, PP2PP and STAR), including batch, mail, printing and data storage. In addition, BNL is the U.S. Tier 1 Center for ATLAS computing, and the RCF also provides for the computational needs of the U.S. collaborators in ATLAS.

The Linux Farm at the RCF provides the majority of the CPU power in the RCF. It is currently listed as the $`3^{rd}`$ largest cluster, according to ”Clusters Top500” (http://clusters.top500.org). Figure 1 shows the rapid growth of the Linux Farm in the last few years.

All aspects of its development (hardware and software), operations and maintenance are overseen by the Linux Farm group, currently a staff of 5 FTE within the RCF.

## II HARDWARE

The Linux Farm is built with commercially available thin rack-mounted, Intel-based servers (1-U and 2-U form factors). Currently, there are 1097 dual-CPU production servers with approximately 917,728 SpecInt2000. Table 1 summarizes the hardware currently in service in the Linux Farm. Hardware reliability has not been an issue at the RCF. The average failure rate is 0.0052 $`failures/(machinemonth)`$, which translates to 5.7 hardware failures per month at its present size. Hardware failures are dominated by disk and power supply failures. A detailed breakdown of the hardware failures by category is shown in Figure 2.

## III SOFTWARE

The Linux Farm at the RCF uses a custom image of RedHat 7.2, modified to conform to the requirements of the RHIC experiments and to the security protocols of BNL. The customized image is installed via KickStart , the RedHat Linux automated installation tool.

The Linux Farm servers are equipped with a variety of compilers (gcc, PGI, Intel) and debuggers (gdb, Totalview, Intel) to provide a large degree of flexibility to its end users. In addition, the servers also support network file systems (AFS, NFS) and batch services, LSF and a RCF-designed software compatible with our MDS system. Figure 3 shows the GUI for the RCF-designed batch software.

Monitoring and control of the cluster hardware, software and infrastructure (power and cooling) is provided via a mix of open-source software, RCF-designed software and vendor-provided software. Figures 4 and 5 display some of our cluster monitoring tools, while figure 6 shows the historical temperature data of the cluster equipment, one of the tools we use to monitor the infrastructure supporting the Linux Farm.

## IV SECURITY

Because of general cyber-security standards at BNL and as part of the general security policy of the RCF, several measures were taken to insure that only authorized users can access the Linux Farm.

A firewall has been placed around the RCF (including the Linux Farm) to minimize security breaches. In addition, users are only allowed interactive access to the Linux Farm via dedicated gatekeeper systems whose software is kept up-to-date to minimize the exploitation of security weaknesses in the software. No other method is provided for interactive access to the Linux Farm. The operating system in the Linux Farm servers has been modified to accomodate these security measures and to enhance our ability to detect unauthorized access.

The RCF has also begun the deployment of a Kerberos 5 single sign-on system that will eventually replace our current authentication system.

## V GRID-LIKE CAPABILITIES

GRID-like technology has evolved from conceptual designs to promising prototypes with real capabilities in the last few years. It is a natural fit to increasingly powerful Linux clusters coupled with geographically diverse end-users and increasingly large data samples typical of large scale high energy & nuclear physics experiments. GRID-like technology is also making significant inroads into industrial applications and has attracted the interest and support of well-known software manufacturers, such as Platform Computing .

The Linux Farm has begun to investigate, install and support (where appropriate) prototypes of GRID-like software that possess capabilities that are of interest to our users, such as Ganglia, Condor and GLOBUS.

### V.1 GANGLIA

Ganglia  is a full-feature, open source distributed monitoring software for high-performance computer clusters. It is based on a hierarchical design targeted at federation of clusters, and it supports clusters up to 2000 servers in size. Early prototypes are already equipped with a end-user Web interface, historical data information and clustering of remote systems. The monitoring data collected by ganglia can be used as the basis of a batch system job scheduler mechanism, although this has not yet been tested.

The Linux Farm has deployed a ganglia prototype for the STAR experiment. In the prototype, the ganglia collector has been configured to gather information from each of the nodes where a ganglia client daemon is running. This master collector then makes this information available to qualified external collectors. Figure 9 illustrates the ganglia deployment within the RCF.

Security issues with this ganglia prototype are being investigated. Of particular concern, currently there is no user-friendly method to restrict the type and amount of information transmitted to external collectors. Wrap-around scripts written by RCF staff were used to restrict the information (see figures 7 and 8). In addition, as more servers are added to the ganglia master collector, scalability issues will become a major concern as well. The Linux Farm group plans to continue to test and expand the scope of the ganglia prototype where appropriate.

### V.2 CONDOR

Condor is a open-source batch software created and supported by the University of Wisconsin . Condor is a full-feature batch software that include features such as job queuing mechanism, configurable scheduling policy, priority scheme, checkpoint capability and resource monitoring & management. Condor can be used to connect remote clusters at geographically diverse locations, so it is a natural fit to the GRID computational philosophy. Condor has an interface to the GRID via Condor-G.

The Linux Farm group is in the midst of upgrading its MDS-compatible batch system to improve reliability & scalability and add functionality. As part of the upgrade, Condor is being evaluated as a job scheduler for the new MDS-compatible batch system. Since media-based MDS systems is expected to play a considerable role at the RCF for the foreseable future, an effort is being made to integrate Condor with the MDS-interface API software. The current batch system does not have an interface to GRID-like architectures, and Condor can add this missing functionality via Condor-G. The basic design of the new batch system is shown in figure 10.

Once the prototype of the upgraded MDS-compatible batch system is installed, Condor scalability studies will be done to understand how performance is affected under the expected heavy usage.

### V.3 GLOBUS & LSF

The ability of users to submit jobs to remote clusters has been one of the principal motivations for the Linux Farm group to explore interfacing our batch system with GRID-like software.

The Linux Farm has a prototype GLOBUS  gatekeeper server that interfaces with LSF for the ATLAS experiment. Authorized users at remote sites can submit jobs to the gatekeeper. The gatekeeper interprets the GLOBUS commands and submits jobs to the proper LSF queues running on the ATLAS Linux cluster at the RCF. A diagram of the prototype is shown in figure 11. Figure 12 shows actual LSF jobs submitted by remote users via the GLOBUS gatekeeper.

Currently, the system is being expanded to include both PHENIX and STAR experiments in RHIC, and additional GLOBUS gatekeepers are being brought on-line.

## VI Near-Term Plans

In the near-term, the current prototypes are expected to expand and slowly mature into production tools, as the outstanding issues are resolved.

Ganglia has already gone through two upgrades within the RCF, and the Linux Farm group is expecting that it will become part of the standard software packages on all its production servers in the near future.

Condor has been evaluated continuously on a small number of servers as the future job scheduler of the upgraded MDS-compatible batch system. Many of the oustanding issues have been resolved or are being studied, and we expect the batch system upgrade to be a year-long project.

The Linux Farm is currently using LSF v.4.2, and we plan to upgrade it to LSF v.5.x together with an OS upgrade in the next few months. New LSF features such as advance resource reservation and GRID membership protocols match well with the GRID computational architecture and can further integrate the GLOBUS gatekeepers with the LSF batch system.

###### Acknowledgements.

The authors wish to thank the Information Technology Division and the Physics Department at BNL for their support to the RCF mission.


format: text

text:
# Using globular clusters to test gravity in the weak acceleration regime

## 1 Introduction

Newton’s law of gravity is routinely used to describe the physical properties of galaxies, even though its validity has been fully and directly verified only within the boundaries of the solar system. At the distance of Pluto, the sun produces an acceleration of $`3\times 10^4`$ cm s-2. By contrast, the acceleration experienced by the sun in the gravitational field of the Milky Way is $`2\times 10^8`$ cm s-2. Thus, any time Newton’s law is applied to galaxies (e.g., to infer the existence of dark matter, hereafter DM), its validity is extrapolated by several orders of magnitude. Although there are in principle no reasons to distrust Newton’s law in the weak acceleration regime, testing it would be interesting on its own. That such a test may be important is also suggested by the strong observational evidence that all spacecraft in the outer solar system and beyond are experiencing an anomalous, unexplained additional acceleration toward the sun (Anderson et al. (1998)).

In addition, a particular modification of Newtonian dynamics known as MOND (Milgrom (1983)) claims to be able to describe many properties of galaxies without invoking DM (McGaugh & de Block (1998); Mortlock & Turner (2001); see Sanders & McGaugh (2002) for a recent review). MOND postulates a breakdown of Newton’s law of gravity (or inertia) below few times $`a_0=1.2\times 10^8`$ cm s-2 (Begeman et al. (1991)). If this is correct any physical system, not just galaxies, should deviate from Newton’s law below this acceleration. Confirming such a breakdown would be of major importance regardless of the validity of MOND.

All these considerations prompted us to look for an experiment involving gravitational acceleration of the order of $`a_0`$. To reach such low accelerations in a laboratory is clearly difficult, so we focused our attention on a class of astrophysical objects useful for this purpose. Globular clusters (GC) are among the smallest virialized structures in the universe and it is currently believed that DM is not affecting their internal dynamics. This belief is based on solid observational results:

1) the existence of tidal tails (Leon et al. (2000)) provides the strongest evidence against large amounts of DM, either baryonic or non-baryonic, in GCs. If present, the potential well would be so deep that no stars could escape the cluster. Numerical simulations (Moore (1996)) show indeed that in order to form tidal tails the maximum mass-to-light ratio consistent with observations is $`M/L<2.5`$ in solar units.

2) In GC studies good agreement is found between dynamical and luminous masses for an old stellar population (Mandushev et al. (1991)).

3) Up to now micro-lensing studies have failed to detect the dense dark objects of GC dimension postulated by cold dark matter hierarchical model and predicted to exist in the Milky Way halo (Navarro et al. (1997); Ibata et al. (2002)).

It is a fortunate coincidence that stars located in the outskirts of GCs are subjected to a gravitational acceleration comparable to the one experienced by stars in the outskirts of galaxies. Therefore, these stars are ideal test particles for probing Newton’s law down to very weak acceleration regimes.

If as suggested by MOND, Newton’s law really breaks down in the low acceleration limit and a different law controls the motion of stars, then the deviation from expectations must follow a systematic pattern. That is, GC dynamics must reproduce exactly the dynamics of galaxies (provided the regime of acceleration is the same). In particular, compact GCs should behave like high-surface-brightness galaxies, i.e., the velocity dispersion profile should be flat at large radii (Carollo et al (1995); Simie & Prugniel (2002)). Loose GCs, on the other hand, should behave like low-surface-brightness galaxies and show large mass discrepancies at any radius (Mateo (1998)).

In this letter, we present the result of a pilot experiment reaching accelerations as low as $`10^8`$ cm s-2.

## 2 A Pilot Experiment on $`\omega `$ Centauri

We studied the dynamical properties of stars located in the outskirt of $`\omega `$ Cen, a nearby, well studied GC, chosen only for being one of the very few clusters for which internal proper motions are known, and because it is visible from the south.

Located at a distance of 6.4 kpc from the Galaxy center, $`\omega `$ Cen is the largest GC known, with mass estimates $`>10^6`$ M, and tidal radius of $`70`$ pc (Meylan & Mayor (1986)). In spite of the large dimension and relative proximity to the Galaxy center, the density of the Galactic halo (Bahcall & Soneira (1980)) is so low that within its tidal radius $`2\times 10^4`$ M of DM are expected, which is dynamically negligible. The cluster is therefore suitable for testing Newton’s law down to the smallest acceleration probed by the motion of its stars.

Available data allow the determination of the dynamical properties of the cluster central regions. In Fig. 1 the radial velocity dispersion profile, as derived using data for $`400`$ stars up to $`30`$ pc from the center (Meylan & Mayor (1986)), is compared to the proper motion dispersion profile derived using data for many thousands stars up to 25 pc from the cluster center (van Leeuwen et al. (2000)). The three components of the velocity dispersion agree well with each other, with no indication of anisotropy. This is very important for two reasons: i) it implies the cluster has reached dynamical equilibrium, making it irrelevant whether the cluster is the result of merging; ii) radial velocities, which are the easiest to determine and the simplest to interpret, are representative of the cluster’s kinematics.

To extend the dispersion profile and reach the desired acceleration regime, we have selected from the list of van Leeuwen et al. 91 stars at distance $`>30`$ pc from the center and membership probability $`>90`$%. Selected stars are located only on the west side of the cluster, covering $`150`$ degrees in azimuth. This limitation is due to the technical details of val Leeuwen et al. observations and is in no way related to the physical properties of $`\omega `$ Cen. Thus this should not bias our result toward large velocity dispersions.

Stars were subsequently observed during August 2001 with the ESO Very Large Telescope (VLT), in order to measure their radial velocity. Spectra were obtained with UVES (UV-visual Echelle Spectrograph) with resolution R=40000. The derived radial velocities have an average error of 1 km s-1, as verified using telluric lines. With these new data, one can in principle study the velocity dispersion of all 3 components of the velocity. However, due to the low velocity of stars in the outskirt of the cluster, proper motion data have average errors comparable to the dispersion we are trying to measure ($`8`$ km s-1). Therefore, though useful for selecting probable cluster members, these data are not suitable for our purpose and we considered radial velocity data only. Having verified the isotropy of the velocity ellipsoid, this limitation does not affect the generality of the result.

## 3 Results for $`\omega `$ Centauri

The average radial velocity of $`\omega `$ Cen is 232 km s-1 (Meylan & Mayor (1986)), allowing the discrimination of cluster members from foreground stars (Fig. 2). Of the initial candidates, 75 were found to be true cluster members (consistent with the assigned membership probability). The remaining have radial velocity of few km s-1 and are most likely part of the Galactic disk.

The radial velocity dispersion as a function of distance is given in Table 1. Columns 1 and 2 give the bin width in pc and the number of stars per bin. Column 3 gives the observed velocity dispersion in km s-1. In order to trace the cluster potential, these values were combined in quadrature with the values for the cluster rotation given in Column 4 (Meylan & Mayor (1986)). This correction is marginal and is relevant only for comparing our data with published data. The final corrected values of the velocity dispersion $`\sigma `$ are in Column 5 These new data extend the study of the kinematics from 30 to 45 pc from the cluster center (projected distances, Fig. 1). We found the dispersion does not decrease monotonically with distance, as expected in Newtonian dynamics, remaining constant at large radii. Invoking a rapid increase of the velocity anisotropy outward of 25 pc from the center to make the profile falling as $`1/\sqrt{r}`$, implies that the other two components of the velocity dispersion have only $`70`$% of the energy of the radial component. Such a large and rapid increase of the velocity anisotropy seems very unlikely. Therefore, though limited to the west side of the cluster, the apparent constancy of the velocity dispersion is intriguing.

## 4 Discussion

This project was started as a pilot experiment to see if at least in one globular cluster interesting results, either in favor or against Newton’s law of gravity, could be found. After the study of $`\omega `$ Cen was completed, we became aware of the existence of a similar set of data showing that also in M15 the velocity dispersion profile remains flat (within errors) at large radii (see Fig. 8 in Drukier et al. (1998)).

The flattening of the velocity dispersion profile at large radii observed in both $`\omega `$ Cen and M15 can be due i) to tidal heating, as explicitly suggested for M15 (Drukier et al. (1998)); ii) to a DM halo surrounding the clusters; iii) to a breakdown of Newton’s law in the weak regime limit.

The first two scenarios are certainly possible. Despite the fact that they require ad hoc assumptions to maintain the profile flat, they can not be ruled out by the present data. Therefore, the purpose of this work is not to disprove these two scenarios, but to fully investigate the third one.

In an attempt to study Newton’s law in a regime were its validity is not proved, the use of masses (required to compute accelerations) derived from dynamical models based on this law would lead to a circular argument. In particular masses derived from the virial theorem or King models may be overestimated (one should keep in mind that both man-made spacecrafts and MOND suggest gravity is pulling harder than expected). To estimate the mass of the clusters we therefore are better off converting the cluster luminosities into masses adopting a suitable mass-to-light ratio $`\tau `$.

The total absolute V magnitude of $`\omega `$ Cen and M15 are $`10.29`$ and $`9.17`$ (Harris (1996)), corresponding to a mass $`M=1.08\tau `$ and $`3.87\tau `$ millions of solar masses, respectively. The mass so estimated can be safely assumed to be all inside the radius at which the velocity dispersion profile becomes flat. In the case of $`\omega `$ Cen the profile flattens at $`27\pm 3`$ pc from the center, where the acceleration is $`2.1_{0.4}^{+0.5}\times 10^8\times \tau `$ cm s-2. In M15 the profile flattens at $`18\pm 3`$ pc (that is $`6\pm 1`$ arcmin from the center, see Fig. 8 of Drukier et al. (1998)), corresponding to $`1.7_{0.5}^{+0.7}\times 10^8\times \tau `$ cm s-2. Very interestingly, provided $`\tau `$ in the same in the two clusters, both profiles start to be flat for the same value (within errors) of the acceleration.

These two clusters have very little in common: they have different masses, different positions and orbits in the galactic halo. Their dynamical evolution was also different. Moreover, it has been claimed that $`\omega `$ Cen could contain two stellar populations being the result of a merging of two cluster (Lee et al. (1999)), or that it is the remnant of a dwarf galaxy (Hilker & Richtler (2000); Carraro & Lia (2000)). Thus we find the fact that the two profiles are so similar a significant one.

To further investigate if this result is in agreement with the claim that deviations from Newtonian dynamics should appear at few times $`a_0`$ (Milgrom (1983)), we need to explicitly adopt a value for $`\tau `$. For a Salpeter’s initial mass function with slope $`x=1.3`$, updated theoretical evolutionary models (Cassisi et al. (1998)) give $`M/L=1.4`$ for simple stellar population. This value would apply to a GC that has retained all its initial mass. Real clusters do experience mass losses because of tidal interaction and evaporation (e.g. Aguilar, Hut, & Ostriker (1988), Smith & Burkert (2002)). The dynamical evolution preferentially removes low mass stars that are the major contributors to the cluster mass, while contributing little to the luminosity. This effect has been verified in those GCs where the present day mass function has been determined, finding a slope $`x=0.7`$ at low mass end (De Marchi (1999); Piotto et al. (1997); Andreuzzi et al. (2000)). With this mass function the same models give $`M/L=0.6`$. Considering that both values depend on the adopted low mass cutoff, we assume $`M/L=1`$ as a fair estimate of the true value. For this $`\tau `$ the velocity dispersion profile of both clusters starts to be flat at $`2a_0`$, consistent with MOND prediction.

GCs are hundreds of times smaller and thousands of time lighter than galaxies, nonetheless the velocity dispersion profile of at least these two clusters precisely mimics, both in shape and absolute acceleration, the one of elliptical galaxies (explained invoking DM; Carollo et al (1995)). There is no reason for the flattening in GCs and galaxies to occur for the same value of the acceleration. If tidal heating is responsible for the observed profile, then we should conclude that the tidal action of the Milky Way not only conspires to make $`\omega `$ Cen and M15 very similar to each other, but also produces exactly the same effect as observed in galaxies. The same argument applies if a DM halo surrounding these two clusters is invoked to explain the effect. As we said, both scenarios are possible but require ad hoc assumptions to explain these coincidences, making them somewhat weak.

The parallel between GCs and galaxies can be pushed even further if we reanalyze a recent result for the GC Palomar 13. It has been reported that this GC has a velocity dispersion of 2 km/s, corresponding to $`M/L40`$ (Côté et al. (2002)). This high dispersion can simply indicate that the cluster is out of dynamical equilibrium, in which case the virial theorem is not applicable and the derived M/L meaningless. Palomar 13, however, moves in an hardly disturbed orbit (Siegel et al. (2001)) and its crossing time ($`2`$ My for an internal dispersion of 2 km s-1) is short enough to ensure stars reach dynamical equilibrium quickly. We therefore argue there are good reasons to believe Palomar 13 is in equilibrium. If this is the case, the data are either reflecting the effects of a massive DM halo or a break down of Newton’s law. Palomar 13 is very loose, having central density of only $`4`$ M pc-3 (Harris (1996)), corresponding to an internal acceleration of gravity $`<10^9`$ cm s-2 all the way to the cluster center. It is therefore similar to low-surface-brightness galaxies which are characterized by very large $`M/L`$ (Mateo (1998)). If our interpretation of the result for $`\omega `$ Cen and M15 is correct, it is therefore not surprising that Palomar 13 shows a large mass discrepancy because this is expected in the case of a break down of Newton’s law below $`a_0`$. Indeed, it as been shown that both Palomar 13 Scarpa (2003) and low-surface-brightness galaxies McGaugh & de Block (1998) follows the predictions of MOND.

## 5 Conclusions

We have presented new measurements of the velocity dispersion of the GC $`\omega `$ Cen, which allow us to trace the gravitational potential down to an acceleration of $`8\times 10^9`$ cm s-2. We have found that the dispersion profile remains flat well inside the tidal radius as soon as the acceleration of gravity approaches $`a_0`$. A similar behavior is also observed in M15, though it was previously ascribed to tidal heating (Drukier et al. (1998)).

This result is surprising and may suggest a failure of Newton’s law at low accelerations. A conclusion also supported by data for Palomar 13.

Although present data can not rule out the standard dark-matter scenario, nor tidal heating effects, the coincidence we pointed out is remarkable. Combined with the anomalous acceleration experienced by spacecrafts in the periphery of the solar system (Anderson et al. (1998)), this result prompts for accurate determination of the internal dynamics of as many GCs as possible. Beside the interest of testing a law of physics in a regime where its validity has never been verified, these observations have the potential of proving or disproving in a direct way the existence of non-baryonic DM.

We wish to thanks L. Pulone and L. Ciotti and S. Cassisi for helpful discussion and suggestion.


text:
# Rotational Evolution of Protoneutron Stars

## 1 Introduction

A neutron star (NS) is born after the iron core collapse of a massive progenitor star ($`>8M_{\mathrm{}}`$) which exhausts its nuclear fuel. The birth of the neutron star often results in a Type-II supernova explosion (Burrows, 2000). Immediately after its birth, the neutron star is hot and lepton-rich because the core is opaque to neutrinos. This young object is called a protoneutron star (PNS). As compared to the neutrino-free case (i.e. the cold neutron star), the trapped neutrinos significantly change the chemical equilibria between baryons and leptons, the fractions of all the compositions, and the equation of state (EOS). Therefore, the Kelvin-Helmholtz epoch of the evolution of the PNS, during which the PNS changes from a hot and lepton-rich compact star to a cold and neutrino-free one, is the most important evolutionary stage of the PNS. The time scale of the Kelvin-Helmholtz epoch completely depends on the neutrino interactions in the hot and dense matter (Reddy et al., 1998). Namely the diffusion time scale of the trapped neutrinos from the core to the outside is roughly several tens of seconds (Prakash et al., 1997). As pointed out in the literature, in principle, observations of supernova neutrinos may constrain the EOSs and the interior composition of neutron stars. (Prakash et al., 1995; Pons et al., 1999, 2001b, 2001a).

The properties of PNSs which contain only ordinary nuclear matter have been investigated in detail by many authors (Takatsuka et al., 1994; Bombaci et al., 1995; Prakash et al., 1997; Hashimoto et al., 1994; Goussard et al., 1997, 1998; Strobel et al., 1999; Sumiyoshi et al., 1999). Because of the uncertainties of the composition of the matter in the interior of neutron stars, the effects of some exotic states, such as hyperonic matter (Prakash et al., 1997; Pons et al., 1999), quark matter (Pons et al., 2001b), kaon condensation (Pons et al., 2001a) and the quark-hadron phase transition (Prakash et al., 1995), on the evolution of PNSs have also been explored. Because the trapped neutrinos increase the chemical potential of electrons, the common consequence of the inclusion of the possible exotic states which contain negatively charged components in the PNS is the existence of a metastable star, the mass of which is larger than the maximum mass for a cold neutrino-free neutron star. Inevitably, the metastable PNS will collapse to a black hole at some time during the Kelvin-Helmholtz epoch. Meanwhile, the neutrino signal from the PNS will cease suddenly, because the time scale of the collapse is on the order the free-fall time which is much shorter than the evolutionary time scale of the PNS.

As the trapped neutrinos escape, the global structure of the PNS changes, especially the momentum of inertia which affects the stellar rotational frequency dramatically. Our key assumption is that Kelvin-Helmholtz timescale over which the neutrinos escape ($`10`$ s) is much longer than the dynamical timescale of the PNS ($`1`$ ms); therefore, the protoneutron star evolves through a series of hydrostatic equilibria with successively smaller neutrino fractions. In this work, we investigate the stellar global properties as a function of the number density of the trapped electron neutrinos. Especially, we contrast the rotational evolution of hyperonic PNSs with ordinary PNSs. The evolutionary behavior of a PNS in whose core the other possible exotic states exist is qualitatively similar to the results for the hyperonic star. In order to describe the hyperonic matter with trapped neutrinos, the relativistic mean field theory (RMFT) has been generally applied (Walecka, 1974; Chin, 1977; Serot, 1979; Serot & Walecka, 1986; Müller & Serot, 1996; Prakash et al., 1997; Yuan & Zhang, 1999; Glendenning, 2000). The potential model is also used sometimes (Prakash et al., 1997). Recently, the potential model with the Bruecker-Bethe-Goldstone many-body theory at zero temperature was employed to model the hyperonic matter with trapped neutrinos, and to investigate the global properties of the PNSs. As expected, the qualitative results which have been obtained in the both models are very similar (Vidaña et al., 2003). The previous works have indicated that the effect of the trapped neutrinos dominates over that of temperature on the internal composition and EOS (Prakash et al., 1997). Thus, for simplicity, we use the RMFT to describe the dense nuclear matter with hyperons at zero temperature, and explore the response of the rotational frequency to the escape of the neutrinos.

After its birth, the neutron star generally spins fast initially, even near to the Keplerian angular velocity, $`\mathrm{\Omega }_\text{K}`$ (Fryer & Heger, 2000). Stellar rotation is essential for the generation of the magnetic field (e.g. Thompson & Duncan, 1993) and affects the global properties of the PNS significantly (Akiyama et al., 2003). In order to to deal with the rotation of a rapidly rotating neutron star, a fully general relativistic method should be used. Several independent techniques exist in the literature (Stergioulas, 1998). Our calculation is based on the Stergioulas & Friedman’s KEH codes (Komatsu et al., 1989a, b), which is available as a public domain code (Stergioulas & Friedman, 1995).

This paper is organized as follows. A brief description of the RMFT for dense matter with hyperons, the microscopic properties of dense matter in $`\beta `$-equilibrium, and the global properties of the non-rotating and rapidly rotating PNSs are given in § 2. In § 3, the numerical results of the evolution of the rotational frequency is presented. The results and discussions will be given in § 4.

## 2 Equations of state for protoneutron stars

### 2.1 Description of the dense matter

In the framework of relativistic field theory, the interaction between nucleons and hyperons is mediated by the exchange of three meson fields, $`\sigma `$, $`\omega `$, and $`\rho `$ mesons. The coupling constants and meson masses can be algebraically related to the bulk properties of the symmetric nuclear matter (Glendenning, 2000). The Lagrangian of the hadrons is given by,

$``$ $`=`$ $`{\displaystyle \underset{i}{}}\overline{\psi }_i\left[i\gamma _\mu \left(^\mu +ig_{\omega i}\omega ^\mu +i{\displaystyle \frac{g_{\rho i}}{2}}𝝉𝝆^\mu \right)(m_ig_{\sigma i}\sigma )\right]\psi _i`$ (1)
$`+`$ $`{\displaystyle \frac{1}{2}}\left(_\mu \sigma ^\mu \sigma m_\sigma ^2\sigma ^2\right)U(\sigma )`$
$``$ $`{\displaystyle \frac{1}{4}}\omega _{\mu \nu }\omega ^{\mu \nu }+{\displaystyle \frac{1}{2}}m_\omega ^2\omega _\mu \omega ^\mu `$
$``$ $`{\displaystyle \frac{1}{4}}𝝆_{\mu \nu }𝝆^{\mu \nu }+{\displaystyle \frac{1}{2}}m_\rho ^2𝝆_\mu 𝝆^\mu ,`$
where,

$`\omega _{\mu \nu }`$ $`=`$ $`_\mu \omega _\nu _\nu \omega _\mu ,`$ (2)
$`𝝆_{\mu \nu }`$ $`=`$ $`_\mu 𝝆_\nu _\nu 𝝆_\mu +g_{\rho n}𝝆_\mu \times 𝝆_\nu .`$ (3)
Throughout this paper we will use the set of natural units where $`\mathrm{}=c=1`$ unless noted otherwise. Here $`\psi _i`$, $`\sigma `$, $`\omega `$, and $`𝝆`$ denote the fields of the baryon of species $`i`$ ($`i=n,p,\mathrm{\Lambda },\mathrm{\Sigma }^{},\mathrm{\Sigma }^0,\mathrm{\Sigma }^+,\mathrm{\Xi }^{},\mathrm{\Xi }^0`$), the mesons $`\sigma `$, $`\omega `$, and $`\rho `$ with the masses of $`m_i`$, $`m_\sigma `$, $`m_\omega `$, $`m_\rho `$, respectively. The constants $`g_{\sigma i}`$, $`g_{\omega i}`$, $`g_{\rho i}`$ are coupling constants for interactions between mesons and baryons. $`𝝉=(\tau _1,\tau _2,\tau _3)`$ denotes the $`2\times 2`$ Pauli isospin matrices and the dot and cross products are calculated over isospin space. The potential of the self-interaction of the scalar field which reads,

$$U(\sigma )=\frac{1}{3}bm_n(g_{\sigma n}\sigma )^3+\frac{1}{4}c(g_{\sigma n}\sigma )^4,$$
(4)
is introduced to produce reasonable incompressibility of nuclear matter (Boguta & Bodmer, 1977). Here the coefficients $`b`$ and $`c`$ denote self-coupling constants for the $`\sigma `$ meson field.

The dynamical equations of nucleons and mesons can be derived from the above Lagrangian. Generally, these equations are coupled. Applying the RMFT approximation (Walecka, 1974), the dynamical equations are decoupled, they are,

$`\mu _i`$ $`=`$ $`g_{\omega i}\omega _0+g_{\rho i}\rho _{03}I_{3i}+\sqrt{k^2+m_i^2},`$ (5)
$`m_\sigma ^2\sigma `$ $`=`$ $`g_{\sigma i}n_{\mathrm{s}i}{\displaystyle \frac{U(\sigma )}{\sigma }},`$ (6)
$`\omega _0`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{g_{\omega i}}{m_\omega ^2}}n_i,`$ (7)
$`\rho _{03}`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle \frac{g_{\rho i}}{m_\rho ^2}}I_{3i}n_i.`$ (8)
Here, $`\mu _i`$ is the chemical potentials of baryons, $`I_{3i}`$ is the third component of isospin for the baryons, and $`m_i^{}=m_ig_{\sigma i}\sigma `$ are the effective masses of baryons. Finally, the scalar density $`n_{\mathrm{s}i}\overline{\psi }\psi `$ is defined as

$$n_{\mathrm{s}i}=\frac{2}{\pi ^2}_0^{k_i^\mathrm{F}}\frac{m_i^{}}{\sqrt{k^2+(m_i^{})^2}}k^2𝑑k,$$
(9)
where $`k_i^\mathrm{F}`$ are the Fermi momenta of baryons. The condition of chemical equilibrium between baryons and leptons is given by,

$$\mu _i=b_i\mu _nq_i\stackrel{~}{\mu _l},$$
(10)
where $`b_i`$ is the baryon number of particle $`i`$ and $`q_i`$ is its charge. For the neutrino-free case, $`\stackrel{~}{\mu _l}=\mu _l`$, for the trapped-neutrino case, $`\stackrel{~}{\mu _l}=\mu _l\mu _{\nu _l}`$. Here $`\mu _l`$ $`(l=e^{},\mu ^{})`$ are the chemical potentials of the leptons. Because neutrinos are trapped, the lepton number per baryon $`Y_{\mathrm{L}_l}`$ of each flavor must be conserved on dynamical time scales,

$$Y_{\mathrm{L}_l}=Y_l+Y_{\nu _l}=\text{constant}.$$
(11)
Studies of the gravitational collapse calculations of the core of a massive star indicate that $`Y_{\mathrm{L}_e}0.4`$ (Prakash et al., 1997). In addition, because no muons appear when neutrinos become trapped, $`Y_{\mathrm{L}_\mu }=0`$. During the Kelvin-Helmholtz epoch, $`Y_{\nu _e}`$ changes from its initial value to zero. Therefore, if the chemical potential of neutrons and electrons are known, the fields of the mesons $`\sigma `$, $`\omega `$, and $`\rho `$ can be solved numerically, and then the momenta of all the baryons can be determined at the same time. Given the baryon number density $`n_\mathrm{B}`$, the additional condition, charge neutrality, is needed to determined $`\mu _n`$ and $`\mu _e`$ self-consistently.

According to the standard procedure of RMFT theory (Walecka, 1974), the energy density and the pressure of the baryons can be obtained,

$`ϵ_\mathrm{B}`$ $`=`$ $`{\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2+{\displaystyle \frac{1}{2}}m_\omega ^2\omega _0^2+{\displaystyle \frac{1}{2}}m_\rho ^2\rho _{03}^2+{\displaystyle \underset{i}{}}{\displaystyle \frac{1}{\pi ^2}}{\displaystyle _0^{p_i^\mathrm{F}}}\sqrt{p^2+m_i^{}}p^2𝑑p,`$ (12)
$`p_\mathrm{B}`$ $`=`$ $`{\displaystyle \frac{1}{2}}m_\sigma ^2\sigma ^2+{\displaystyle \frac{1}{2}}m_\omega ^2\omega _0^2+{\displaystyle \frac{1}{2}}m_\rho ^2\rho _{03}^2+{\displaystyle \underset{i}{}}{\displaystyle \frac{1}{\pi ^2}}{\displaystyle _0^{p_i^\mathrm{F}}}{\displaystyle \frac{p^4dp}{\sqrt{p^2+m_i^{}}}}.`$ (13)
The coupling constants $`g_{\sigma N}`$, $`g_{\omega N}`$, $`g_{\rho N}`$, $`b`$, $`c`$ ($`N=n,p`$) in RMFT are algebraically related to the bulk properties of symmetric nuclear matter at saturation density. In our calculation, we choose the so-called GM1 set (Glendenning & Moszkowski, 1991),

$$\left(\frac{g_{\sigma N}}{m_\sigma }\right)^2=11.79,\left(\frac{g_{\omega N}}{m_\omega }\right)^2=7.149,\left(\frac{g_{\rho N}}{m_\rho }\right)^2=4.411,b=2.947\times 10^3,c=1.07\times 10^3,$$
(14)
from which the nuclear properties arise as follows: the saturation density, effective mass of nuclear, incompressibility, and binding energy per nucleon are $`n_0=0.153`$ fm-3, $`m^{}=0.7`$, $`K=300`$ MeV, and $`B/A=16.3`$ MeV.

At the saturation density, the bulk properties are independent of the hyperon couplings which could be determined by reproducing the binding of the $`\mathrm{\Lambda }`$ hyperon in nuclear matter (Glendenning & Moszkowski, 1991). It is assumed that all hyperons in the octet have the same coupling which are expressed in terms of nucleon couplings,

$$x_{\sigma H}=\frac{g_{\sigma H}}{g_{\sigma N}},x_{\omega H}=\frac{g_{\omega H}}{g_{\omega N}},x_{\rho H}=\frac{g_{\rho H}}{g_{\rho N}}.$$
(15)
In the GM1 set (Glendenning et al., 1997; Glendenning & Moszkowski, 1991; Glendenning, 2000),

$$x_{\sigma H}=0.6,x_{\omega H}=x_{\rho H}=0.653.$$
(16)
### 2.2 The interior composition and the EOS

The evolution of the PNS with time should be investigated by solving the neutrino transport equations (Pons et al., 1999; Janka et al., 2003; Rampp & Janka, 2002). Due to the escape of the neutrinos, $`Y_{\nu _e}`$ decreases with time. Because our purpose in this work is to explore the rotational behavior of the PNS as the neutrinos escape, we study the rotational frequency and other interesting quantities as a function of $`Y_{\nu _e}`$. The ratio of the number fraction of electron neutrinos to its initial value is assumed to be a constant with baryon number density, for a given time. Because the density in the stellar core is almost constant, this is a good approximation. As the electron neutrinos escape, the muon neutrino appears due to the following reaction,

$$e^{}\mu ^{}+\overline{\nu }_\mu +\nu _e.$$
(17)
To describe the escape of the muon neutrinos, we assume that the ratio of the number of the escaped muon neutrinos to that of the muon neutrinos which should escape freely if the matter is transparent to neutrinos equals to $`Y_{\nu _e}/Y_{\nu _e}^i`$. Our calculations show the number density of muons is generally several orders of magnitude lower than that of electrons, so its effects on the EOS of PNS could be neglected, even though we have taken it into consideration.

Previous studies have shown that neutrinos are trapped above a density $`n_{\mathrm{env}}6\times 10^4`$ fm-3 (Burrows et al., 1995; Strobel et al., 1999). Therefore, before the neutrinos escape, we set $`Y_{\mathrm{L}_\mathrm{e}}=0.4`$ for densities above $`n_{\mathrm{env}}`$. For the description of neutrino-free subnuclear matter, the FPS EOS (Lorenz et al., 1993) is chosen, which is smoothly matched to the EOS of nuclear matter described in RMFT. At low baryon number density ($`n_\text{B}<0.1`$ fm-3), the nucleons are non-relativistic, thus the pressure contributed by the trapped neutrinos and the relativistic eletrons dominates. The EOS of subnuclear dense matter becomes much stiffer due to the trapped neutrinos (see Fig. 1).

The main effect of the neutrino trapping in the dense matter above nuclear density is that the fraction of electrons increases dramatically, regardless of the interior composition. Therefore, for a protoneutron star which contains only nucleons in its core, the fraction of protons is forced to increase in the same manner because of charge neutrality. Comparing to the neutrino-free case, the symmetry energy which is contributed by the $`\rho `$ meson and linearly depends on the difference between the baryon number density of neutrons and protons, decreases significantly (see Eq. 8 and Eq. 12-13). Consequently, the EOS in the normal stellar core is softer with trapped neutrinos than without them. For hyperonic matter, the higher chemical potential of electrons that results from trapped neutrinos increases the critical baryon number density of the appearance of $`\mathrm{\Sigma }^{}`$ particles. Because there are fewer species of baryons, the Fermi momenta of each species is larger, so the pressure of hyperonic matter with trapped neutrinos is larger than the pressure without neutrinos (see Fig. 2).

Figure 3 shows the composition of matter in the stellar interior which contains hyperons. The subfigures in Fig. 3 are the results for $`Y_{\nu _e}/Y_{\nu _e}^i=100\%,60\%,20\%,0`$, respectively. As noticed in previous work, the neutrino trapping causes the following results: (1) muons do not appear; (2) the number density of electrons and protons increases significantly; (3) the onset of the negatively charged hyperons takes place at higher baryon number density, chargeless hyperons appear at nearly the same baryon number density, while the onset of positively charged hyperons occurs at lower density (see Fig. 3). The whole effect is that the EOS of the hyperonic matter becomes stiffer at high density. With the escape of the trapped neutrinos, muons begin to appear and $`\mathrm{\Sigma }^{}`$ particles appear at lower density.

For the purpose of comparison, Fig. 4 shows the composition of the matter which includes only nucleons. As in the above case, the neutrino trapping suppresses the appearance of muons and increases the number density of electrons and protons.

The responses of the EOSs based on the hyperon model (NPH) and nucleon model (NP) as the neutrinos escape are shown in Fig. 5. In the NP model, as expected, the EOSs become more and more stiff as the neutrinos escape. The same occurs in the NPH model, below the critical density of the emergence of the hyperons. However, in the NPH model, at densities above the critical density, the EOSs become softer due to the appearance of the more and more additional components.

### 2.3 Properties of the non-rotating and rapidly rotating PNSs

In order to solve the Einstein field equations with the source terms, neutron-star matter is generally assumed to behave as a perfect fluid. For a static neutron star, the resulting equations are well known, the Tolman-Oppenheimer-Volkoff equations. It is complicated to deal with the rotation of a neutron star. If the stellar angular velocity $`\mathrm{\Omega }`$ is small compared to the critical value $`\mathrm{\Omega }_0=\sqrt{4\pi Gϵ_c}`$, where $`ϵ_c`$ is the mass density at the center of the star, the perturbation theory developed by Hartle and Thorne is accurate enough (Hartle, 1967; Hartle & Thorne, 1968; Chubarian et al., 2000). Hartle’s original perturbation theory was improved by considering the effects of centrifugal stretching and frame dragging in an attempt to calculate the properties of a rapidly rotating neutron star with $`\mathrm{\Omega }\mathrm{\Omega }_0`$ (Glendenning, 1992; Weber & Glendenning, 1992). Nevertheless, the exact treatment of the rotation of a compact star should be done in full general relativity. Several reasonable assumptions are generally made, they are, the spacetime and the matter are stationary and axisymmetric, and the star uniformly rotates. Under these assumptions, the metric of the spacetime can be expressed by

$$ds^2=e^{\gamma +\rho }dt^2+e^{2\alpha }(dr^2+r^2d\theta ^2)+e^{\gamma \rho }r^2\mathrm{sin}^2\theta (d\varphi \omega dt)^2,$$
(18)
where the metric potentials $`\gamma ,\rho ,\omega `$ and $`\alpha `$ are functions of $`r`$ and $`\theta `$ only. Several independent methods to numerically calculate the metric potentials exist. After solving the metric potentials, any coordinate-invariant physical quantity can be calculated.

Figures 6-8 show the global properties of the neutron stars. In Fig. 6, we plot the gravitational mass ($`M`$) as a function of the central energy-density ($`ϵ_c`$) for PNSs with different fractions of trapped neutrinos $`Y_{\nu _e}/Y_{\nu _e}^i`$ and with different rotational frequencies $`\mathrm{\Omega }=0`$, $`\mathrm{\Omega }=4\times 10^3`$, and $`\mathrm{\Omega }=\mathrm{\Omega }_\mathrm{K}`$ in both the NP and NPH model. As predicted in the previous works, the maximum mass of the static PNS decreases as the neutrinos escape in the NPH model, which results in the possible existence of a metastable star. In the NP model the maximum mass increases as the neutrinos escape. These results hold even when the PNS spins very rapidly.

Figures 7-8 show the radii and the moments of inertia as a function of the rest mass respectively. In the NP model, if the stellar rest mass $`M_0`$ is near the maximum mass, the radii $`R`$ and the moments of inertia $`I`$ increase as the neutrinos escape, while decrease in the NPH model. Evidently, this behavior is dominated by the change of the EOSs in the stellar core. The trapped neutrinos always make the EOS in the crust stiffer (see Fig. 1) because of the dominance of the lepton pressure. When the initial mass is small, firstly, the PNS shrinks as the neutrino escape the crust, then begins to expand when the pressure of the $`\rho `$ meson which determines the symmetry energy is comparable to that of the leptons.

Figure 8 show the moments of inertia as a function of the rest mass. In the NP model, if the stellar rest mass $`M_0`$ is near the maximum mass, the moments of inertia $`I`$ increase as the neutrinos escape, as the radii $`R`$ do. In the NPH model, the moments of inertia decrease. When the initial mass is small, firstly, $`I`$ decreases then increases while the neutrinos escape. The change in behavior is not as significant as what happens to the radii, because gravitational mass increases as the neutrinos escape the core. In some sense, this compensates the decrease of the radius.

## 3 Rotational evolution of protoneutron stars

As investigated in the above section, the global structure of the PNS changes in different ways as the neutrinos escape depending on whether hyperons exist in the core. This may affect the rotational evolution of the PNS. During the evolution of the PNS, the baryonic mass is fixed, if there is no accretion. Neglecting the loss of angular momentum due to the radiation of gravitational and electromagnetic waves, the loss of the total angular momentum of the PNS only results from the escaping of the neutrinos. To determine how the loss of angular momentum carried away by the neutrinos affects the rotational evolution, we also assume that the neutrinos do not carry any angular momentum as they escape.

Figure 9 shows the stellar radii and moment of inertia as a function of neutrino fraction. For the NP model, the radii decrease first mainly because the EOSs of the subdense matter in the stellar crust becomes softer as the neutrinos escape. The radii then increase mainly because the EOS of the dense matter in the stellar core becomes stiffer, as contribution from symmetry energy increases. For the NPH model, the radii keep decreasing because both the decrease of the neutrino pressure and the appearance of more and more hyperons in the stellar core make the EOS softer. If $`M_0=1.5M_{\mathrm{}}`$, hyperonic matter does not appear in the stellar core at all. If $`M_0=2.2M_{\mathrm{}}`$, the PNS is metastable, and it collapses to a black hole at the end of the Kelvin-Helmholtz epoch. The evolution of the total moment of inertia is very similar to that of the radii.

Figure 10 shows how the rotational frequency changes as the neutrinos escape in the NPH and NP model. Early in the Kelvin-Helmholtz epoch, the PNS always spins up independent of the presence of hyperons, because as neutrinos escape the envelope, the PNS shrinks. The crustal EOS affects the rotational evolution of a protoneutron star significantly (Cheng et al., 2002). Later in the evolution, roughly speaking, for the NPH model, the PNS spins up with the EOS in the stellar core becoming softer; while for the NP model, the PNS spins down with the EOS in the stellar core becoming stiffer. For $`M=2.0M_{\mathrm{}}`$, in the middle of the evolution, the NPH spins down for a short period. Even though the increased number of hyperons make the EOS softer, at the same time the fraction of the proton decreases, which makes the EOS stiffer. These two effects compete with each other, the spin down results from the increase of the symmetry energy. If the PNS is metastable, the star keeps spinning up significantly, and then collapses to a black hole at the end of its life. It is clearly shown in Fig. 10 that even though the escaped neutrinos carry away some of the angular momentum which decreases the rotational frequency of the PNS, the effect of the change of the global structure due to the escape of the neutrinos dominates the stellar rotational evolution.

The change of the rotational frequencies of the PNS after the neutrinos escape would change the initial distribution of the periods of neutron stars. Assuming the distribution at the beginnning of the Kelvin-Helmholtz epoch is uniform between zero and the Keplerian frequency ($`5344s^1`$ for $`M_0=2.0M_{\mathrm{}}`$), the new distribution is shown in Fig. 11. For PNSs without hyperons, the range of the periods is somewhat narrower than before, but roughly speaking, the distribution is still almost uniform. For the PNS with hyperons, a significant number of slowly rotating PNSs shift to higher spins; consequently, the presence of hyperonic matter at the center of neutron stars will skew the distribution of initial spin periods toward shorter periods.

Figure 12 shows the evolution of the integrated flux of the escaped electron neutrinos. The total energy of the trapped electron neutrinos $`E_\nu ^{\mathrm{tot}}`$ may be estimated as follows: the average baryon number density $`n_\mathrm{B}M/(\frac{4}{3}\pi R^3m_n)`$, the corresponding neutrino number density $`n_{\nu _e}=n_BY_{\nu _e}`$, and the Fermi energy of the trapped neutrinos $`E_{\nu _e}^\mathrm{F}=p_{\nu _e}^\mathrm{F}=(6\pi ^2n_{\nu _e})^{1/3}0.2\text{GeV}`$, therefore, $`E_\nu ^{\mathrm{tot}}`$ is given by,

$$E_{\nu _e}^{\mathrm{tot}}\frac{3}{4}n_{\nu _e}E_{\nu _e}^\mathrm{F}\times \frac{4}{3}\pi R^32.12\times 10^2\left(\frac{M}{2.0M_{\mathrm{}}}\right)^{4/3}\left(\frac{R}{1.5\times 10^6\mathrm{cm}}\right)^1\left(\frac{Y_{\nu _e}}{0.1}\right)^{4/3}M_{\mathrm{}}.$$
(19)
It is evident that the total neutrino flux is of the same order for the different initial masses. The flux of the electron neutrinos in both the NPH model and NP model is also similar, because the hyperons only appear at a higher density when neutrinos are trapped, which makes the hyperonic PNS and the normal PNS very similar.

## 4 Conclusions and discussion

In this paper, we investigate the rotational evolution of protoneutron stars which contain only nucleons with that of PNSs which contain hyperonic matter at zero temperature in full general relativity. It is found that PNSs contract and spin up at the early stage of its evolution as the trapped neutrinos escape. During this stage the contribution of leptons to the pressure in the stellar crust dominates. This conclusion is independent of the interior composition. As the neutrinos continue to escape, the evolution of the PNS mainly depends on the changing behavior of the EOS of the dense matter in the stellar core; therefore, it differs in the two different models. A PNS which contains only nucleons stops contracting and begins to expand because the EOS becomes stiffer as neutrinos escape from the core. PNSs that contain hyperonic matter keep shrinking, and their spin frequency increases because the EOS becomes softer with the escape of the neutrinos. At the end of the evolution, for PNSs without hyperons, the range of the spin periods becomes a little bit narrower than the initial one. However, the shape of the distribution of the spin periods is very similar to that of the initial distribution. For PNSs with hyperons, the distribution of the initial spin periods is skewed significantly toward shorter periods. If a PNS is metastable, it keeps spinning up significantly on the neutrino diffusion timescale before it collapses to a black hole.

In addition to the hyperonic matter, other exotic states, such as quark matter, kaon condensation and others might exist in the stellar core. Because the neutrino trapping makes the onset of these exotic states which contain negatively charged particles take place at the higher baryon number density, a PNS with such an exotic composition spins up as the neutrinos escape. This characteristic phenomenon is common among them.

In this work, the effect of the temperature is ignored. Certainly regardless of the interior composition, heat introduces another energy source which would make the EOS stiffer than in the zero temperature case. The effect of the cooling of the PNS on the evolution of the PNS could be predicted. However, since the change of the chemical equilibrium among the interior particles dominates the change of the EOS, the cooling of the PNS could not change the results we obtain qualitatively.

In principle, the rotation period of a PNS could be observed by a ground-based gravitational radiation detector, such as LIGO. Fryer et al. (2002) examine several avenues for gravitational wave emission from core collapse. If the core forms a bar or breaks in two, the analysis that we have presented here is inapplicable because we have assumed that the core remains axisymmetric (Eq. 18). Furthermore, the third mode that they considered is the ringdown of the black hole; in this case, the protoneutron-star stage is no longer visible. However, they also have estimated the role that $`r`$modes may play in gravitational wave generation (Ho & Lai, 2000). The frequency of the observed gravitational radiation reflects the spin frequency of the protoneutron star, so the waveform would reveal the internal structure of the protoneutron star as long as fallback did not strongly affect the spin of the core during the Kelvin-Helmholtz epoch. It is unclear whether the $`r`$modes will have a chance to grow during this early epoch (Fryer et al., 2002), but if they do LIGO-II could detect gravitational radiation from supernovae within 10 Mpc.

As estimated by Eq. 19, a huge energy ($`10^{52}`$ ergs) is released in the form of high energy neutrinos during the Kelvin-Helmholtz epoch. If the neutrino flux from the PNS is modulated with the stellar rotation, a neutrino telescope might observe temporal structure at the stellar period, which could provide strict constraints on models of dense matter. The typical energy of the emitted neutrinos is a few MeV (Arnett et al., 1989). Twenty neutrinos were detected from SN1987A (Hirata et al., 1987; Bionta et al., 1987) over about a ten-second interval. Currently, the total sensitivity of detectors sensitive to supernova neutrinos is thirty times larger than in 1987, so one would expect to detect about 15,000 neutrinos from a galactic supernova. Ransom et al. (2002) give straightforward formulae to estimate the minimum pulsed fraction detectable with a given certainty from a sample of arrival times. If we assume that the spin frequency of the PNS lies between zero and 1500 Hz and that the spin may vary by 200 Hz during the ten-second duration of the neutrino burst, a blind search for the spin frequency and its evolution would require about $`3\times 10^7`$ trials. If one wants to detect a signal with a false alarm probability of $`P_0`$, the minimum sinusoidal pulsed fraction that one could detect is

$$f_p=2\left[\frac{\mathrm{ln}(N_{\text{trial}}/P_0)1}{N}\right]^{1/2}=0.08$$
(20)
for a false alarm probability of $`10^3`$ and 15,000 detected neutrinos. Because the expected change in frequency is small compared to the expected spin frequency, an unaccelerated search only does marginally better with a minimum detectable pulsed fraction of 0.06. Although these pulsed fractions are several times larger than required for pulsar kicks (e.g. Spruit & Phinney, 1998; Arras & Lai, 1999), Earthbound neutrino detectors are most sensitive to neutrinos in the Wien tail of the thermal distribution where the pulsed fraction will be larger than the momemtum anisotropy.

A third observational probe of these models is the distribution of initial periods of neutron stars. This may be obtained from the distribution of current spin periods of radio pulsars along with an independent estimate of the ages of the neutron stars. One must assume that only magnetic dipole radiation has contributed to the spin down, and furthermore that the standard model for magnetic dipole radiation is correct (Kaspi & Helfand, 2002). If gravitational radiation is important, then such estimates of the initial spin distribution probe the gravitational radiation epoch (see Owen et al., 1998; Ho & Lai, 2000; Watts & Andersson, 2002; Arras et al., 2002, for a few viewpoints) rather than the protoneutron-star epoch discussed here. Only a few neutron stars have a good estimates of their initial periods, which are typically an order of magnitude longer than the Kepler period. If the emission of gravititional radiation is not important to the spin evolution of neutron stars, this would argue against a skewed distribution of neutron-star initial spins; however, so few neutron stars have estimates of their initial periods, and the contribution of gravitational radiations is uncertain, it is premature to make any definite conclusions.

The most spectacular manifestation of the rotational evolution of a protoneutron star would be the formation and subsequent collapse of a metastable stellar core. How the core collapses into a black hole and the observable consequences of that collapse are beyond the scope of this paper. We will examine the process in a subsequent paper. For a protoneutron star that rotates rapidly, an important issue is whether the entire protoneutron star can collapse on a dynamical time or portions of the star cannot fall directly into the black hole and form an accretion disk which evolves on a viscous timescale. This could provide a detailed model of the central engines of collapsars, a successful explanation of gamma-ray bursts (MacFadyen et al., 2001).

We have studied the structure of rapidly rotating axisymmetric protoneutron stars using a fully general relativistic treatment of the stellar structure and relativistic mean-field theory to model the equation of state of dense matter. We have contrasted the spin evolution of normal protoneutron stars with hyperonic protoneutron stars. In general the evolution protoneutron stars with other exotic species will be similar to that of hyperonic neutron stars. We have found that hyperonic protoneutron stars can be metastable; that is, the maximum mass of a stable hyperonic protoneutron star decreases as the neutrinos escape. Both normal and hyperonic protoneutron stars initially spin up as neutrinos escape the outer layers of the protoneutron star and the low-density region shrinks. As neutrinos begin to escape the core of a hyperonic protoneutron star, the core continues to spin up. The metastable protoneutron star continues to spin up as it begins to collapse to a black hole. This contrasts with behavior of a normal protoneutron star. As neutrinos stream out the core of a normal protoneutron star, the equation of state stiffens, the core expands and spins down. This spin evolution leaves a signature on both the neutrino emission and gravitational radiation from the stellar collapse that might be observable from nearby supernovae. Finally, if the processes leading to the formation of hyperonic and normal protoneutron stars are similar, we would expect that the hyperonic neutron stars initially to spin faster than normal neutron stars.

We would like to thank the discussions with Prof. Ramesh Narayan, Prof. K.S. Cheng, Prof. J.L. Zhang, Dr. Scott Ransom and Dr. John Bahcall. Y.F.Y acknowledges the hospitality of Harvard-Smithsonian Center for Astrophysics. Y.F.Y is partially supported by the Special Funds for Major State Research Projects, and the National Natural Science Foundation (10233030). Support for J.S.H. was provided by the National Aeronautics and Space Administration through Chandra Postdoctoral Fellowship Award Number PF0-10015 issued by the Chandra X-ray Observatory Center, which is operated by the Smithsonian Astrophysical Observatory for and on behalf of NASA under contract NAS8-39073.


text:
# Solitonic transmission of Bose-Einstein matter waves

## Abstract

We consider a continuous atom laser propagating through a wave guide with a constriction. Two different types of transmitted stationary flow are possible. The first one coincides, at low incident current, with the non-interacting flow. As the incident flux increases, the repulsive interactions decrease the corresponding transmission coefficient. The second type of flow only occurs for sufficiently large incident currents and has a solitonic structure. Remarkably, for any chemical potential there always exists a value of the incident flux at which the solitonic flow is perfectly transmitted.

The transport properties of matter confined to small structures display distinct quantum effects qualitatively different from those observed at macroscopic scales. These are grounded on global phase coherence throughout the sample and can be, in many cases, understood within single particle pictures, without referring to any specific details of the system. As a result, they arise in many different fields (electronic systems, atomic physics, electromagnetism, acoustics) . Some examples are (weak and strong) localization, Bloch oscillations and conductance quantization.

Recent experimental developments in the physics of Bose-Einstein condensation (BEC) of dilute vapor (in particular the microchip guiding technique) open up the prospect of studying coherent transport phenomena using guided atom lasers . Besides, because of the extraordinary control over these systems, they offer a unique opportunity to go beyond the single particle behavior, and to study specific effects induced by interaction. In the present article we focus on a simple situation, where a BEC matter wave propagates through a guide with a constriction . By an adiabatic approximation, the three–dimensional flow is reduced to one dimension, where the atoms now feel, due to the constriction, a longitudinal step–like potential of height $`V_0`$. In the absence of interaction, the transmission $`T`$ does not depend on the incident current but only on the beam’s energy; $`T`$ is always lower than unity and tends to this limit when the energy of the beam is large compared to $`V_0`$. In the following we consider atoms with a repulsive effective interaction characterized by a scattering length $`a_{sc}>0`$. The most salient features of the flow are all at variance with respect to the non-interacting case: (i) the transmission coefficient depends on the current, (ii) at given chemical potential, there exists a maximum transmitted current above which no stationary flow exists, (iii) at a given current, several distinct stationary solutions with different $`T`$ are possible and (iv) for any chemical potential larger than $`V_0`$, there is a particular value of the incident current which induces total transmission.

Consider a continuous atom laser incident on a constriction of a waveguide. Within the adiabatic approximation  the transverse motion is restricted to the lowest transverse eigen-state. The constriction affects the longitudinal motion via an effective step-like potential whose magnitude is fixed by the ground state energy of the transverse Hamiltonian. For a BEC system, the adiabatic approximation implies that the condensate wave function can be cast in the form (see Ref. ),

$$\mathrm{\Psi }(\stackrel{}{r},t)=\psi (x,t)\varphi (\stackrel{}{r}_{};n;x),$$
(1)
where $`\psi (x,t)`$ describes the motion along the axis of the laser (the beam is flowing along the positive $`x`$ direction). $`\varphi `$ is the equilibrium wave function (normalized to unity) in the transverse $`(\stackrel{}{r}_{})`$ direction. It depends parametrically on the longitudinal density $`n(x,t)=d^2r_{}|\mathrm{\Psi }|^2=|\psi (x,t)|^2`$. The beam is confined in the transverse direction by a trapping potential $`V_{}(\stackrel{}{r}_{};x)`$, which is $`x`$-dependent in the region of the constriction. Then, the longitudinal wave equation reads  (in units where $`\mathrm{}=m=1`$)

$$\frac{1}{2}_{xx}\psi +\left\{V_{}(x)+ϵ[n(x,t);x]\right\}\psi =i_t\psi .$$
(2)
In (2), $`V_{}(x)`$ represents an effective longitudinal potential due to the constriction. If, to be specific, we consider a transverse harmonic confinement with pulsation $`\omega _{}(x)`$, then $`V_{}(x)=\omega _{}(x)\omega _{}(\mathrm{})`$ (energy is measured with respect of the ground state energy of the non interacting transverse Hamiltonian far before the constriction). $`ϵ(n;x)`$ is a nonlinear term describing the mean field interaction averaged over a transverse slice of the beam. One has $`ϵ(n;x)=2\omega _{}(x)na_{sc}`$ in the low density regime ($`na_{sc}1`$), and $`ϵ(n;x)=2\omega _{}(x)\sqrt{na_{sc}}`$ in the high density regime ($`na_{sc}1`$) .

Our purpose is to determine the transmission of steady state solutions of (2) where $`\psi (x,t)=\mathrm{exp}\{i\mu t\}A(x)\mathrm{exp}\{iS(x)\}`$, with $`A`$ and $`S`$ real functions. The density is $`n=A^2`$ and the local velocity is $`v=dS/dx`$. From Eq. (2) one obtains (i) flux conservation: $`n(x)v(x)`$ is a constant that we denote $`J_{\mathrm{}}`$, and (ii) a Schrödinger-like equation for the amplitude:

$$\frac{1}{2}\frac{d^2A}{dx^2}+\left\{V_{}(x)+ϵ[n(x);x]+\frac{J_{\mathrm{}}^2}{2n^2(x)}\right\}A=\mu A.$$
(3)
To define the scattering problem one needs to study the asymptotic behavior of the flow far from the constriction. Far upstream $`V_{}(x\mathrm{})=0`$ and the non-linear term in (2) looses its explicit $`x`$ dependence, taking the simpler form $`ϵ[n(x,t)]`$. Thus, in this region, (3) admits a first integral of the form  :

$`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{dA}{dx}}\right)^2+W[n(x)]`$ $`=`$ $`E_{cl}`$ (4)
$`\text{with}W(n)`$ $`=`$ $`\epsilon (n)+\mu n+{\displaystyle \frac{J_{\mathrm{}}^2}{2n}},`$ (5)
where $`\epsilon (n)=_0^nϵ(\rho )𝑑\rho `$ and $`E_{cl}`$ is an integration constant. Eq. (4) has a simple interpretation in terms of classical dynamics. It expresses the energy conservation of a fictitious classical particle with “position” $`A`$ and “time” $`x`$, moving in a potential $`W`$; $`E_{cl}`$ being the total energy of this particle. Eq. (4) is thus integrable by quadrature, and the density profile can be deduced from the plot of Fig. 1. Small values of $`E_{cl}W(n_1)`$ correspond to small density oscillations, whereas the highest acceptable value is $`E_{cl}=W(n_2)`$, corresponding to a gray soliton.

In the far down-stream region, $`V_{}(x+\mathrm{})`$ also takes a constant value $`V_0=\omega _{}(+\mathrm{})\omega _{}(\mathrm{})>0`$. Hence, (3) admits, in this region, a first integral analogous of (4) where, due to the change in $`\omega _{}`$, $`ϵ(n)`$ (resp. $`\epsilon (n)`$) takes a different form which we denote $`ϵ_0(n)`$ (resp. $`\epsilon _0(n)`$). The new form of $`W(n)`$ is denoted $`W_0(n)`$ and the new constant of integration is $`E_{cl}^0`$:

$`{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{dA}{dx}}\right)^2+W_0[n(x)]`$ $`=`$ $`E_{cl}^0`$ (6)
$`\text{with}W_0(n)`$ $`=`$ $`\epsilon _0(n)+(\mu V_0)n+{\displaystyle \frac{J_{\mathrm{}}^2}{2n}}.`$ (7)
It follows from general arguments on the dispersion of elementary excitations of Eq. (2) that the physically acceptable boundary conditions of Eq. (3) correspond to a constant far down-stream density (see ). The asymptotic $`x+\mathrm{}`$ density should thus be equal either to $`n_{1,0}`$ (we denote this as “case A”), or to $`n_{2,0}`$ (case B); $`n_{1,0}`$ and $`n_{2,0}`$ – being the analogous of $`n_1`$ and $`n_2`$ of Fig. 1 – are extrema of $`W_0(n)`$. They are solutions of

$$\mu =ϵ_0(n)+V_0+\frac{J_{\mathrm{}}^2}{2n^2}.$$
(8)
In the non-interacting case the term $`\epsilon _0(n)`$ is absent from $`W_0(n)`$ which has only one minimum ($`n_{10}`$). Case A is therefore the only possible solution in non-interacting systems. Case B describes new non-perturbative effects related to interaction. It corresponds to an asymptotic down-stream density which is part of a gray soliton.

There exists a maximum value $`J_{\mathrm{}}^{max}`$ of $`J_{\mathrm{}}`$ above which Eq. (8) admits no solution: as $`J_{\mathrm{}}`$ is increased (keeping $`\mu `$ and $`V_0`$ fixed), the two extrema of $`W_0(n)`$ move toward each other, until they coalesce and disappear. This marks the onset of a time–dependent flow. If, to be specific, we consider the case $`ϵ_0(n)=g_0n^{\nu _0}`$, then

$$J_{\mathrm{}}^{max}=\left[\frac{2}{\nu _0+2}(\mu V_0)\right]^{\frac{1}{\nu _0}+\frac{1}{2}}\sqrt{\nu _0}g_0^{1/\nu _0}.$$
(9)
The scattering process is now well defined. It corresponds to the matching between two asymptotic densities described by the classical motion of a particle of energy $`E_{cl}`$ in a potential $`W(n)`$ at $`x\mathrm{}`$, and of energy $`E_{cl}^0`$ in a potential $`W_0(n)`$ at $`x+\mathrm{}`$ (with $`E_{cl}^0`$ either equal to $`W_0(n_{1,0})`$ or to $`W_0(n_{2,0})`$). Eq. (3) being non-linear, an important question is how to properly define a transmission and a reflection coefficient; i.e., is it possible to disentangle an incident and a reflected wave in the upstream flow ? We follow here an approach closely related to usual experimental set-ups, and choose to work with an incident and a reflected beam which can be approximated by plane waves. This corresponds to a regime where Eqs. (3,4) can be linearized in the far upstream region. In this regime, for $`x\mathrm{}`$, we write $`n(x)=n_1+\delta n(x)`$ and expand $`nW(n)`$ to second order in $`\delta n`$. Then Eq. (4) leads to

$$\left(\frac{d\delta n}{dx}\right)^2+\kappa _1^2\delta n^2=8(n_1+\delta n)[E_{cl}W(n_1)],$$
(10)
where $`\kappa _1^2=4(v_1^2c_1^2)`$, $`v_1=J_{\mathrm{}}/n_1`$ being the average velocity of the up-stream beam and $`c_1=[n_1(dϵ/dn)_{n_1}]^{1/2}`$ the sound velocity of a beam with constant density $`n_1`$. The linearization (10) is valid provided $`|\delta n(x)/n_1|\kappa _1^2/c_1^2`$. In this regime, if one further imposes $`v_1c_1`$, the up-stream density oscillations can be analyzed in term of incident and reflected particles (and not quasi-particles). This allows to unambiguously define the incident, reflected, and transmitted current as $`J_i=(n_1+\delta n_1/2)v_1`$, $`J_r=\delta n_1v_1/2`$ and $`J_t=n_1v_1=J_{\mathrm{}}`$ (where $`\delta n_1=4[E_{cl}W(n_1)]/\kappa _1^2`$). Hence, once $`E_{cl}`$ is known, the transmission at given incident current $`J_i`$ is determined through

$$E_{cl}=W(n_1)+\frac{\kappa _1^2}{4}\delta n_1=W(n_1)+\frac{\kappa _1^2}{2v_1}J_i(1T).$$
(11)
The linearization procedure explained so far is valid in the case of small upstream interaction (this is the essence of the condition $`v_1c_1`$). However, all interaction effects are fully taken into account in the down-stream region, where they are indeed more important (the constriction acts as a barrier which lowers the velocity of the down-stream flow and, by flux conservation, increases its density ). In the following we solve the exact non-linear equation (3) and use the linearization procedure only to define the transmission coefficient $`T`$. Thus, the results presented below are of very general validity, but their analysis in term of transmission coefficient is only correct in so far as the linearization procedure is valid.

The method is now the following: for a given $`\mu `$ and $`J_i`$, assume a particular value of $`T`$. This determines $`J_{\mathrm{}}=TJ_i`$, fixes the form of the function $`W_0(n)`$, the value of $`n(+\mathrm{})`$ (it is equal to $`n_{1,0}`$ in case A and to $`n_{2,0}`$ in case B) and of $`E_{cl}^0=W_0[n(+\mathrm{})]`$. Integrating (3) backwards from $`x=+\mathrm{}`$ to $`x=\mathrm{}`$ yields $`E_{cl}`$, which should be compatible with (11). If not, the value of $`T`$ has to be modified until self-consistency is achieved.

To understand the physical picture we consider an abrupt step-like constriction. In this geometry, numerical integration of (3) can be bypassed because $`E_{cl}`$ is simply expressed in terms of $`E_{cl}^0`$ (see Eq. (12)). However, the adiabatic approximation (1) is based on the assumption that the typical longitudinal length scale is much larger than the transverse one. This is clearly violated by an abrupt constriction. Nevertheless, in certain parameter ranges, the adiabatic approximation remains valid. In order to illustrate this point we compute the transmission of non-interacting atoms , for which the exact solution can be obtained numerically. We thus consider a linear wave moving in a guide with harmonic confinement whose transverse pulsation changes abruptly (at $`x=0`$ say) from $`\omega _{}^<`$ (up-stream) to $`\omega _{}^>=\alpha \omega _{}^<`$ (down-stream), with $`\alpha >1`$. The incident atoms occupy the transverse ground state of the up-stream potential. The numerical solution of the problem can be worked out by a straightforward 3D generalization of the procedure devised in Ref.  for studying a similar 2D problem. We denote the transmission $`T^L`$, the superscript recalling that we are in a linear (i.e., non-interacting) regime. The result for $`T^L`$ is presented in Fig. 2 for $`\alpha =2`$ and $`\alpha =3`$. The vertical bars indicate the location of the energies of the reflected (thin lines) and transmitted (thick lines) channels. Conservation of angular momentum along the longitudinal axis imposes selection rules between channels. These rules effectively forbid half of the energetically allowed channels, and those henceforth do not play any role in the transmission.

Within the adiabatic approximation, the constriction is described by an abrupt longitudinal step potential of height $`V_0=\omega _{}^>\omega _{}^<`$. The corresponding transmission is $`T_{adia}^L=4[\mu (\mu V_0)]^{1/2}(\sqrt{\mu }+\sqrt{\mu V_0})^2`$ (represented by a dashed curve in Fig. 2). The onset of transmission occurs at $`\mu =V_0`$, i.e., $`\mu /\omega _{}^<=\alpha 1`$. One notices in the figure that when $`\mu `$ is increased from this value, the adiabatic approximation is initially quite accurate. However, at larger values of $`\mu `$ deviations from adiabaticity are clearly visible. For $`\alpha =2`$, a sudden lowering of $`T^L`$ occurs when a new reflected channel opens at $`\mu /\omega _{}^<=2`$ (in all the following we denote as “open channels” those allowed by energy conservation and symmetry rules). From there on, $`T^L`$ diminishes until a new transmission channel opens, at $`\mu /\omega _{}^<=5`$. At this point, a sudden increase of $`T^L`$ is observed. At large values of $`\mu `$, $`T^L`$ tends to unity, as it should. For $`\alpha =3`$, the process is similar, but the breakdown of the adiabatic approximation occurs earlier because the opening of the initial transmitted channel – at $`\mu /\omega _{}^<=2`$ – coincides with that of the first allowed excited reflected channel.

In the following we will concentrate on the region $`\mu /\omega _{}^<\stackrel{>}{}\alpha 1`$ where the adiabatic approximation is well justified and where, as we shall now see, non-linear effects may induce strong modifications of the transmission.

We now turn to the non-linear problem and consider the case $`V_{}(x<0)=0`$, $`ϵ(n;x<0)=ϵ(n)`$ and $`V_{}(x>0)=V_0`$, $`ϵ(n;x>0)=ϵ_0(n)`$. Eq. (3) admits the first integral (4) for all $`x0`$ and (6) for all $`x0`$. $`E_{cl}`$ is determined through $`E_{cl}^0`$ by imposing continuity of $`A`$ and $`A^{}`$ at $`x=0`$:

$$E_{cl}W[n(0)]=E_{cl}^0W_0[n(0)].$$
(12)
Let’s consider case A first. The asymptotic down-stream density is $`n_{1,0}`$ and thus one has, for all $`x0`$, $`n(x)=n_{1,0}`$ (the fictitious classical particle remains at the bottom of the potential well $`W_0`$). In particular, $`n(0)=n_{1,0}`$ and the matching (12) determines $`E_{cl}`$ uniquely. The value of $`T`$ is denoted $`T^A`$ in this case.

Case B is more interesting because the structure of the down-stream solution is richer: $`n(x0)`$ being part of the profile of a gray soliton, $`n(0)`$ can be varied continuously provided the matching (12) is fulfilled at a value acceptable for Eq. (4). Effectively, the only restriction imposed is that $`n_{min}<n(0)<n_{max}`$ ($`n_{min}`$ and $`n_{max}`$ are defined in Fig. 1). As a result, for fixed $`\mu `$ and $`J_i`$, $`E_{cl}`$ is not uniquely determined by $`n(+\mathrm{})`$, and the transmission $`T^B`$ varies between 0 and a value that we denote as $`T_{max}^B`$.

To be specific, we consider a continuous beam of 23Na atoms propagating through a guide with a transverse confinement $`\omega _{}^<=2\pi \times 2`$ kHz (in the region $`x<0`$), to which we impose a narrowing $`\omega _{}^>=2\pi \times 6`$ kHz (in the region $`x>0`$). This represents a barrier of height $`V_0=192`$ nK. In the non-interacting case the transmission $`T^L`$ as a function of $`\mu `$ is plotted in the bottom part of Fig. 2 ($`\alpha =3`$). We take $`\mu =210`$ nK (this corresponds to the kinetic energy of atoms having a velocity of 1.2 cm/s). This value of $`\mu `$ corresponds, in the non-interacting case, to a regime where the adiabatic approximation holds and yields a transmission $`T^LT_{adia}^L0.70`$. The repulsive interaction between atoms introduces in Eq. (3) a non-linear term which, in the region $`x<0`$, reads $`ϵ(n)=gn`$ with $`g=2a_{sc}\omega _{}^<530`$ nK.nm. In the region $`x>0`$, the transverse frequency of the guide is multiplied by 3, and thus $`ϵ_0(n)=g_0n`$ with $`g_0=3\times g`$. An important parameter of the system is the maximum transmitted current $`J_{\mathrm{}}^{max}`$ above which no stationary flow can exist in the down-stream part of the guide. From Eq. (9) (with $`\nu _0=1`$) one obtains $`J_{\mathrm{}}^{max}1.6\times 10^4`$ atom/s.

Fig. 3 summarizes the results obtained. The linearization condition $`v_1c_1`$ is extremely well satisfied in the whole range of incident currents considered (the less favorable case occurs at large $`J_i`$, where $`v_125c_1`$). The horizontal dashed line is the value of the transmission coefficient of non-interacting atoms : $`T^L0.7`$. It is current independent. When the current is increased from zero, $`T^A`$ decreases from this value down to $`T^A0.5`$. At this point (located with a black spot on the figure), $`J_{\mathrm{}}(=T^AJ_i)`$ is equal to $`J_{\mathrm{}}^{max}`$, and a stationary flow of type A is no longer permitted (actually it bifurcates to a type B solution). The prominent feature of the behavior of $`T^A`$ as a function of $`J_i`$ is its decrease compared to the non-interacting value $`T^L`$. The physical reason behind this phenomenon is simple: the available kinetic energy necessary to step over the barrier is reduced when the interaction energy increases, i.e., when the incident current increases. This picture is supported by a perturbative treatment which accurately describes the flow at low incident current ($`J_iJ_{\mathrm{}}^{max}`$) and confirms that the decrease of $`T^A`$ corresponds to an increased fraction of the interaction energy in the chemical potential.

Case B being mediated via interaction, does not exist for low current. It exists only above a critical current $`(2g_0)^1(\mu V_0)^2(8\mu )^{1/2}1.5\times 10^3`$ atom/s $`0.1J_{\mathrm{}}^{max}`$. From this point, $`T_{max}^B`$ increases rapidly up to 1 (reached at $`J_i0.45J_{\mathrm{}}^{max}`$ in the case of Fig. 3), and then decreases down to a point where one can show that it exactly meets the end point of $`T^A`$. From there on, the value of $`T^B`$ is limited by the condition that the flow should be stationary, and one has $`T_{max}^B=J_{\mathrm{}}^{max}/J_i`$, which coincides with the dashed hyperbola in Fig. 3. We emphasize that stationary solutions of type $`B`$ with arbitrary transmission $`0T^BT_{max}^B`$ exists for any current above the critical one.

The nonlinear transport induced by the repulsive two-body interaction has therefore a non-trivial consequence: new solutions – of solitonic character – emerge; they allow for an increased transmission. This contrasts with the behavior of case A where the transmission is lowered by the interaction. One can show that, for any value of $`\mu >V_0`$, there always exists a value of $`J_i`$ such that complete transmission exists in case B. The profile for $`T_{max}^B=1`$ consists of a constant up-stream density $`n(x0)=n_1`$ connected at $`x=0`$ to half a soliton. Fig. 4 displays the density profiles of two stationary flows (case $`A`$ and case $`B`$ at $`T^B=T_{max}^B=1`$) at an incident current $`J_i=0.36J_{\mathrm{}}^{max}`$ (see Fig. 3). Note the significant difference in the densities of the up-stream and down-stream profiles, a purely nonlinear effect mediated by a solitonic profile. We have performed numerical computations that show that the same type of solution also exists for smooth constrictions and that they are dynamically stable, as confirmed by a Bogoliubov analysis. The interactions can thus have two different and, in some sense, opposite consequences on the transport properties of a condensate flow. They diminish the transmission in some instances (case A), but also allow for new stationary flows that can be perfectly transmitted (case B) .

We have restricted our analysis to stationary configurations, but dynamical effects are certainly of interest. Amongst these, one could address the question of the transient that exists before stationarity is reached, or the nature of the flow in parameter regions where stationarity is not possible. An other open problem, that clearly deserves investigation, is the dynamical selection of the different types of stationary profiles discussed above: given an initial low density flow (which, from Fig. 3, is of type A), which branch (A or B) will be followed when the incident flux increases ?

There are different ways to experimentally realize the effect discussed in the present work, namely enhanced solitonic transmission of matter waves. We have studied one possible implementation, where the step–like potential in the longitudinal motion of the condensate is produced by a constriction of the guide. An other possibility is to apply a blue-detuned laser beam on the region $`x>0`$ of a condensate propagating along a guide of constant diameter. In this case, the characteristics of the flow and the barrier should be easily controlled by modifying the laser’s frequency, intensity and waist, thus allowing for a neater experimental observation of the above predicted transport phenomena .

We acknowledge stimulating discussions with D. Guéry-Odelin. L.P.T.M.S. is Unité Mixte de Recherche de l’Université Paris XI et du CNRS, UMR 8626.


text:
# Probing the Edge of the Solar System: Formation of an Unstable Jet-Sheet

## 1 Introduction

Our solar system presents a unique local example of the interaction between a stellar wind and the interstellar medium. As the Sun travels through the interstellar medium it is subject to an interstellar wind. Three different discontinuities are formed due to the interaction between the supersonic solar wind flow and the interstellar wind: termination shock, heliopause and, in the case that interstellar wind is supersonic, a bow shock. At the termination shock, the supersonic solar wind passes to a subsonic regime. The heliopause is a tangential discontinuity that separates the charged components of the two winds. Beyond the termination shock, the solar wind velocity is no longer constant, but decreases with radial distance until it reaches the heliopause. A bow shock may form in the interstellar wind since observations indicate that the solar system travels with a supersonic velocity of approximately $`25km/s`$ (Frisch, 1996). The complex system created by the interaction of the solar wind with the interstellar wind is known as the global heliosphere.

In recent years, there has been increased interest in heliospheric studies in modeling the solar and interstellar interaction, in part due to the two successful Voyager missions, launched in 1977 and in early 2003, at $`87AU`$ and $`69AU`$. There is the expectation that in the near future the Voyager spacecraft, traveling toward the heliopause, will cross the termination shock (Stone et al., 1996). Accurate modeling of this region will allow us to interpret the data taken during the passage of Voyager through (and beyond) the termination shock (for a review on the subject see Suess (1990)). Recent works include elaborate computational models (Baranov & Malama, 1993; Pauls et al., 1995; Pauls & Zank, 1996; Zank et al., 1996; Muller et al., 2000; Liewer et al., 1996; Pogorelov & Matsuda, 2000; Linde et al., 1998; Izmodenov et al., 2001; Washimi & Tanaka, 1996, 2001). Nearly all the models, however, ignored magnetic field effects, especially the solar magnetic field. While the solar magnetic field does not play a major dynamic role in the interplanetary region (Linde et al., 1998), this is not true beyond the termination shock. In the present paper, we focus on the interaction region between the solar wind and the interstellar medium, paying special attention to the region beyond the termination shock.

The solar magnetic field, under quiet solar conditions, is approximately a dipole with opposite polarities in the north and south hemispheres. The supersonic solar wind pulls the solar field out into space creating a heliospheric current sheet separating magnetic field regions of opposite polarities. The heliospheric current sheet corresponds to the region where the azimuthal field, $`B_\varphi =0`$. In the present study we neglect the tilt of the magnetic field axis with respect to the rotation axis. Therefore, near the Sun the current sheet remains at the equatorial plane. An important question overlooked in previous studies was the fate of the current sheet beyond the termination shock. Previous numerical studies by Linde et al. (1998) and Washimi & Tanaka (1996, 2001), that used a 3D MHD simulation which included both the solar and interstellar magnetic field reported results indicating that the current sheet beyond the termination shock remains in the equatorial plane. These studies also neglected the tilt of the magnetic axis with respect to the rotation axis.

The rotation of the Sun coupled to a supersonic solar wind twists the interplanetary magnetic field lines such that the magnetic field, upstream of the termination shock (toward the Sun), takes the shape of an Archimedian spiral (known as the Parker field) (Parker, 1958). In analytic studies of the region beyond the termination shock, Nerney et al. (1991, 1993) predicted the presence of magnetic ridges due to the compression of the azimuthal interplanetary field. Their studies were made in the kinematic approximation where the magnetic field back reaction on the flow was neglected. The flow plasma $`\beta _V`$ (ratio of solar wind ram to magnetic pressures) decreases by a factor of $`10^3`$ near the stagnation point, indicating that a full MHD calculation in the region is necessary. Since there is no natural symmetry in the system, a 3D calculation is required. The enormous range of scales, from hundreds of $`AU`$ to less than $`1AU`$ at the current sheet, requires a time-dependent adaptive grid code, making the calculation computationally very challenging.

## 2 Simulation Results

### 2.1 Description of the Model

We used a full 3D MHD adaptive grid code, BATS-R-US, developed at the University of Michigan. This is a parallel multi-processor code that employs an upwind finite volume scheme based on an approximate Riemann solver for MHD. The single fluid MHD equations are solved in a fully three-dimensional coupled manner. BATS-R-US makes use of a block-based adaptive grid which allows the model to resolve structures of widely varying length scales (see Powell et al. (1999) and Gombosi et al. (2002)).

The inner boundary of our calculation was placed at $`30AU`$ and the outer boundary was on a $`3000\times 3000\times 2400AU`$ box. The grid is designed to resolve in detail the region between the termination shock and the heliopause, especially refining at the current sheet. We used 4.5 million cells ranging from scales of $`1AU`$ at the current sheet to $`36AU`$ at the outer boundary with 6 levels of refinement. The parameters of the solar wind at the inner boundary were: $`n=7.8\times 10^3cm^3`$, $`T=1.6\times 10^3K`$ and a magnetic field of $`B=2\mu G`$ at the equator. For the interstellar medium we used(Frisch, 1996): $`n=0.07cm^3`$, $`v=25km/s`$ and $`T=10^4K`$. The solar wind velocity at the inner boundary was assumed to be isotropic, with $`v=450km/s`$. The magnetic field for the solar system at the inner boundary is taken to be a Parker spiral (with the polarity of the solar magnetic field corresponding to this field being anti-parallel to the axis of solar rotation). In order to eliminate the effect of reconnection at the heliopause in the present study, we did not include an interstellar magnetic field. Only the ionized component was treated.

### 2.2 Formation of a “Jet-Sheet”

Figure 1 is a meridional cut (x-z) at $`t=35`$ years, showing the contours of the magnetic field. The calculation was initiated at $`t=0`$ when the interstellar wind reached the solar wind and the bow shock, heliopause and the termination shock have formed with the current sheet in the equatorial plane. The coordinate system is such that the interstellar wind is flowing from the negative x direction and the solar rotation axis is in the positive z direction. We observe at $`x210AU`$ the presence of magnetic ridges (red) where the magnetic field increases to values of $`2.4\mu G`$.

Beyond the termination shock (TS), we find that a “jet-sheet” forms. As the flow decelerates, beyond the TS, approaching the heliopause (HP), it compresses the azimuthal magnetic field, producing the magnetic ridges. In the current-sheet region, due to the absence of an azimuthal magnetic component, there is no magnetic pressure to slow down the flow and the solar wind streams freely. This leads to the formation of a jet in the meridional plane and sheet in the equatorial plane. Due to the shear between the flow in the “jet-sheet” and the flow in the surrounding medium, the “jet-sheet” (and the current sheet) becomes unstable. We are the first to report this phenomena.

Figure 2a is a meridional cut at $`t=35`$ years showing the velocity contours between TS and the HP. We can see the jet structure: in the equatorial region the wind velocity is much faster than the surrounding medium. The red contours with positive velocities $`(U_x20km/s)`$ indicate solar wind material that was pushed aside by the jet and is flowing towards the TS. This produces the turbulence shown in the black lines that follow the flow streamlines. It can be seen that at $`x=210AU`$, the flow streams at the equator with a velocity $`130km/s`$ while the surrounding medium flows with a velocity of $`40km/s`$. The flow in the current sheet at the TS is subject to the de Laval nozzle effect (Landau & Lifshitz, 1987). As the nozzle widens (at $`x=180AU`$) the flow velocity increases, becoming supersonic, and then decreases approaching the HP. The sonic Mach number $`v/c_S`$, where $`v`$ is the flow velocity and $`c_S`$ is the sound velocity, is 1.1 at $`x=180AU`$ decreasing to 0.2 at $`x=300AU`$. Fig 2c shows a line plot of the velocity in the equatorial plane at $`t=35`$ years. Figure 3 shows the time evolution of the current sheet for nine different times. By $`t=71years`$, the current sheet starts to bend to positive z’s while the nose moves toward negative z. This bending increases the turbulence beneath the current sheet. The current sheet, after an initial bending to the south, reverses its direction and moves to the north. At later times, the HP is highly distorted as can be seen by the flow streamlines in Figure 3g. The flow in the HP is very complicated: more than one stagnation point occurs above the equator and the two fluids appear to intermix. Figure 4 shows a 3D image of the magnetic field contours and lines at $`t=346`$ years.

### 2.3 Kelvin-Helmholtz-like Instability

Due to the velocity shear between the “jet-sheet” and the surrounding medium, a Kelvin-Helmholtz-like (KH) instability develops. The KH growth rate, for an infinite thin layer, is(Chandrasekhar, 1961) $`\mathrm{\Gamma }=0.5\times kU_0[1(2c_A\widehat{k}\widehat{B}_0)^2/(\widehat{k}U_0)^2]^{1/2}`$, where $`U_0`$ is the velocity jump across the shear layer, $`c_A`$ is the Alfvén speed, $`\widehat{k}`$ is the direction of the wavevector (here $`\widehat{x}`$) and $`\widehat{B}_0`$ is the direction of the magnetic field. Maximum growth occurs for $`k1/a`$, $`a`$ is defined as the distance between the maximum velocity and zero (see Figure 2b). For $`t=32years`$ for example, at $`x=210AU`$, $`U_0100km/s`$, and $`a20AU`$. With these values the predicted KH growth rate is $`3.3\times 10^8sec^1`$. The necessary condition for the KH instability is $`c_A<U_0/2B_0`$. At $`x=210AU`$, $`c_A/U_0<0.5`$ for $`z<2AU`$ so the instability can occur near (and at) the current sheet. Measuring the displacement of the current sheet with respect to the equator, at $`x=210AU`$, during the evolution of the instability, gives a linear growth rate of $`1.6\times 10^9sec^1`$ (between $`t=0126`$ years), comparable to the estimate for the KH instability in slab geometry. The wavelength of the perturbation $`\lambda _x2\pi /k63AU`$, using $`k1/a`$ is also consistent with our results.

To determine whether this instability was not observed in prior studies due to the coarseness of the grid, we made an identical calculation, using a less refined grid at the current sheet. For the coarse grid, the cell size at the current sheet was $`46AU`$ compared with the refined grid of $`1AU`$. No instability was seen. The current sheet remained in the equator, consistent with the previous numerical calculations. (e.g,Washimi & Tanaka (1996) used a spherical grid with a resolution varying from $`\mathrm{\Delta }R6AU`$ to $`\mathrm{\Delta }R20AU`$ from the inner to the outer boundary). The line plot in Figure 2b compares the velocity profiles for our coarse and the refined grid cases. For the coarse grid, the velocity of the jet is four times wider ($`(30AU)`$ half width) and smaller by a factor of two. The linear KH growth rate $`\mathrm{\Gamma }kU_0`$ and thus we expect a decrease of about a factor eight. The instability is apparently suppressed in the coarse grid case not only because the growth rate decreases but also because numerical dissipative effects wash it out.

## 3 Discussion

In the present study, we have neglected the effect of the tilt in the solar dipole ($`15^{}`$ from the rotation axis), as in ud-Doula & Owocki (2002) and the previous heliospheric studies by Linde et al. (1998) and Washimi & Tanaka (1996). Time-dependent effects such as solar rotation (27 days) can affect our results. Due to the tilt between the magnetic and the solar rotation axis, the current sheet near the Sun takes the shape of the “ballerina skirt” and can change the region beyond the termination shock. We we plan to investigate how this affects the instability. Past work has shown that the effect of the inclusion of the neutrals coupled to the plasma through charge exchange interactions has several effects such as diminishing the distances to the TS, HP, BS and affecting the heliospheric flow velocity. We intend to include self-consistently the neutrals to study their effect on the jet-sheet structure. The formation of an unstable jet-sheet at the edge of the solar system, can potentially affect, the entry of energetic particles and the intermixing of interstellar and stellar material.

The interaction between stellar winds and the interstellar medium produces several detectable astrophysical structures. Hydrogen walls have been observed in several nearby stars(Wood et al., 2001), including our own solar system(Linsky & Wood, 1996). Stellar wind bow shocks, have been observed on a wide variety of astrophysical scales, associated with pulsars(Cordes et al., 1993), OB-runaway stars(Van Buren & McCray, 1988) and neutron stars(Chatterjee & Cordes, 2002). We predict the formation of an unstable jet-sheet in other stellar systems, as well. The ingredients necessary are a magnetized fast rotating star, possessing a supersonic wind that moves with respect to the interstellar medium. The probable candidates are rotating O stars with hot winds having high velocity in respect to the interstellar wind. The strength of the “jet-sheet” instability is proportional to the shear velocity between the flow at the current sheet and in the surrounding medium that is being slowed by the magnetic ridges. The intensity of the magnetic ridges, depends on the magnetic field, the wind ram pressure and on the rotation rate of the star. Therefore, we predict that the instability will be enhanced for fast rotating stars with strong magnetic fields. Because the compression of the stellar azimuthal field is also dependent on the interstellar pressure, the instability will be further enhanced for fast moving stars relative to the interstellar medium.

We thank M. Velli for helpful discussions. Much of this work was performed at the Jet Propulsion Laboratory of the California Institute of Technology under a contract with NASA. The University of Michigan work was also supported by NASA. GT was partially supported by the Hungarian Science Foundation (OTKA, grant No. T037548).