text:
# Hot fragmentation of nuclei11footnote 1 Invited Talk at Seventh International Conference on Nucleus-Nucleus Collisions, Strasbourg, France, July 3 - 7, 2000

## 1 INTRODUCTION

Fragmentation is a term commonly used to specify a nuclear disassembly by force. Hot fragmentation is meant to indicate the most violent of these processes, following excitations beyond the limits of nuclear binding, but still ending with bound nuclear fragments of different sizes in the final channels. The formation mechanism of these fragments, whether they are the remnants of an incomplete destruction or the products of a condensation (’selforganization’, cf. Ref. ) out of the disordered matter, has continued to be the topic of very active research in recent years .

A two-stage scenario has proven to be fruitful for the interpretation and modelling of hot-fragmentation reactions. It is motivated by the differences of the wave lengths and time scales governing the entrance and exit channels and justified by the remarkable success of statistical approaches for the second stage . The intermediate states are not necessarily equivalent to hot nuclei but should be, more generally, viewed as systems of highly excited nuclear matter, populating a phase space characterized by global quantities like mass, charge, energy, density or temperature. Molecular-dynamics calculations can now continuously follow the reaction process from the first encounter to the final disassembly stages without the need to specifically assume equilibration . In that sense, they present a challenge to the statistical two-stage picture. On the other hand, the extraction of thermodynamical parameters from transport-model calculations has demonstrated that a connection to the statistical approach can be established .

There is little doubt that the densities are low and the temperatures high in the intermediate state, coinciding with the values predicted for the coexistence region of liquid and gaseous nuclear matter. To search for observable links of the multifragmentation phenomenon to the predicted phase transition has therefore been a major motivation for many experimental and theoretical activities. It is evident from the titles and abstracts submitted to this conference that there is growing confidence that signals of this phase transition are actually being observed.

This talk divides into three main parts: a brief review of the reactions that are being used to produce hot nuclear systems, a discussion of the caloric curve of nuclei as a potential signature of the liquid-gas phase transition in finite systems, and new experimental results for the conditions at breakup in spectator reactions. This last part will mainly reflect recent activities of the ALADIN collaboration.

## 2 USEFUL REACTIONS

Thermal fragmentation denotes a concept of studying the breakup of thermally excited nuclear systems formed in collisions of relativistic hadrons with heavy target nuclei. It is based on the expectation, tested with intranuclear-cascade calculations, that light energetic projectiles will generate a statistical (thermal) disorder without exciting collective modes such as compression, rotation or shape deformations. The latter are known to break nuclei very efficiently , and thus would mask the nuclear response to the thermal excitation.

The potential and the limits of this approach have been explored by the ISiS collaboration in their recent experiments at the AGS in Brookhaven with a variety of primary and secondary beams including pions and antiprotons in the momentum range up to 14.6 GeV/c . The properties of intermediate systems produced in these reactions are summarized in Fig. 1. Target-like residues with several hundreds of MeV excitation have the largest cross sections, but excitation energies exceeding 1 GeV can be reached, even though with rapidly dropping cross sections. There is also a loss of mass caused by the heating mechanism at relativistic energies. Some of the nucleons participating in the cascading processes are too energetic to remain part of the intermediate system.

The multi-fragment channels open up at the higher excitation energies, as illustrated in Fig. 2. The transition from residue production to hot fragmentation, highlighted by the dotted vertical lines, is, most notably, associated with a striking decrease of the emission time scales to values in the vicinity of 50 fm/c.

According to Fig. 1, the most efficient projectiles are antiprotons of high momentum which only are difficult to use because of their low abundance in the secondary beam . Some extra heating is generated by the pion cascades from the $`\overline{\mathrm{p}}`$ annihilation which leads to considerable excitation energies already at much lower $`\overline{\mathrm{p}}`$ momenta. The reported fragment multiplicities, however, are significantly smaller for the less energetic antiproton beams even if identical bins of excitation energy are selected . It is therefore an open question whether the excitation energy by itself is the only parameter that governs the decay properties and fragment production. Similar questions have recently been raised by other authors who find fewer fragments in the experiment than are predicted by statistical models . While this may be partly connected to the difficulties inherent with experimentally determining the thermalized excitation energy , it is nevertheless obvious that the solution to this problem will help us to better understand how fragments are formed .

The limits of generating excitation energy can be overcome with composite projectiles. Systematic sets of data with projectiles of different mass and energy have been collected by the FASA , EOS  and KEK/HIMAC  collaborations, and new results are reported in contributions to this conference. The continuing rise of the cross section for high excitations with increasing mass of the collision partner is demonstrated in Fig. 3 for the fragmentation of 197Au projectiles studied by the ALADIN collaboration. These data are shown as a function of the variable $`Z_{\mathrm{bound}}`$, representing the sum of the atomic numbers $`Z_\mathrm{i}`$ of all projectile fragments with $`Z_\mathrm{i}`$ 2, which is inversely correlated with the excitation energy . As an example of the $`Z_{\mathrm{bound}}`$ scaling, a prominent feature of spectator reactions up to very high energies , the isotope temperatures $`T_{\mathrm{HeLi}}`$ are also given in the figure. They are close to 6 MeV for the major part of the $`Z_{\mathrm{bound}}`$ range, but tend to higher values at low $`Z_{\mathrm{bound}}`$, i.e. at the highest excitation energies reached in these reactions. The temperatures depend only on $`Z_{\mathrm{bound}}`$ but not on the specific target that is used to fragment the 197Au projectile.

Large systems at even higher excitations can be produced in central collisions of heavy systems. In this case, the idea of excitation as a simple heating process has to be abandoned, however. Collective modes, compression as well as the directed outward motion of particles and clusters from primary collisions generate an explosive pattern, quantified as collective radial flow . The production of large clusters is rare and inversely correlated with the observed amount of flow. An extreme case, observed at the AGS, has been reported very recently. For 197Au beams of 11.5 GeV/c per nucleon, centrally colliding with heavy targets, the fragment mass yields are steeply exponential with a penalty factor of about 50 for each additional mass unit . The characteristic transition from power law to exponential spectra, as radial flow sets in , has been reproduced with molecular dynamics calculations, very recently e.g. with quantum molecular dynamics , but to fully understand the clusterization mechanism in the dynamical environment still remains an interesting problem for future research .

Heavy symmetric systems below the threshold of collective radial flow (about 50 MeV per nucleon) have been extensively studied by the Miniball/Multics and INDRA collaborations. High-statistics data permit the selection of single-source formation which occurs with small cross sections in central collisions . In more peripheral encounters, several sources contribute to the fragment yields with a clear enhancement in the mid-rapidity domain . The potential mechanisms of mid-rapidity emissions, as e.g. neck formation, and the isospin effects associated with it  are closely connected to several dynamical and statistical aspects of hot fragmentation, among them the isotopic separation expected in the liquid-gas coexistence zone . These topics are the subjects of other plenary talks at this conference.

## 3 LIQUID-GAS PHASE TRANSITION

It is commonly accepted that extended nuclear matter should exhibit a liquid-gas phase transition, following from the Van-der-Waals-like range dependence of the nuclear force . Sharp discontinuities in the infinite system are expected to broaden as the system size decreases . This, however, does not contradict the existence or prevent the identification of phase transitions in small systems with constituent numbers on the nuclear scale .

The appropriate experiment for identifying the liquid-gas phase transition in finite nuclei has recently been done, theoretically. Fermionic as well as antisymmetrized molecular dynamics (FMD, AMD) models were used to study the equilibrium dynamics of small nuclear systems . None of the necessary ingredients were missing in these studies, a container to confine the system, a controllable heating technique, long propagation times to allow the system to settle into equilibrium, and a suitable technique of measuring the temperature. As a remarkable result it was shown that the nuclear dynamics, as represented in these advanced transport-type models, generate a phase transition without any further assumptions.

Caloric curves obtained in the FMD study are shown in Fig. 4. The extracted temperatures exhibit a plateau that extends over about 8 MeV per nucleon of excitation energy and has its origin in the coexistence of liquid and nuclear phases in the system. The densities reflect the transition from a liquid phase in equilibrium with its surrounding vapor to the pure vapor phase. The properties of the system in these asymptotic states were identified as those of a Fermi liquid and a Fermi gas (Van-der-Waals gas in the AMD). Both groups have also demonstrated that the external conditions, the container which controls the pressure, have a strong influence on the properties of the transition. Changing the confinement gives the possibility to map out the phase diagram. The latent heat of the 16O system approaches zero at $`T`$ 10 MeV which may be associated with the critical temperature for that system.

There is no external confinement in the real experiment, and the time scales of hot-fragmentation reactions are rather short. It is therefore even more surprising that the temperature-energy correlation measured for the breakup states exhibits such a similar behavior. The data for 197Au fragmentation shown in Fig. 5 represent the results for 600 MeV per nucleon , with small modifications due to additional experimental information and corrections, and the results for 1000 MeV per nucleon obtained more recently . The temperature of the transition region is close to that obtained with the dynamical  or statistical models  and does not change with the bombarding energy. In contrast to it, the energy associated with the spectator source increases by, on the average, 30% over the range 600 to 1000 MeV per nucleon, a behavior inconsistent with the universality of the spectator decay that so clearly appears in other variables . It is caused by the energy dependence of the mean kinetic energies of nucleons in the spectator frame and most likely indicates that contributions from the early stages of the reaction have been included in the calorimetry for the spectator source. As a consequence, the apparent latent heat deduced even for the lower energy of 600 MeV per nucleon should be considered as an upper limit.

Breakup temperatures and energy contents of the fragmenting system were measured and correlated by several other groups for a variety of different reactions . The basic methods were identical with some variation in the approximations that had to be made. Differences exist, e.g. in whether and how the effects of sequential decays on the temperature were taken into account, whether the neutron multiplicities and kinetic energies were measured or estimated, and whether and how preequilibrium components were identified and explicitly excluded. The last point is part of the bigger problem of identifying and properly selecting the fragmenting source which is crucial.

The resulting caloric curves have several features in common but are considerably different in detail. A deviation from the behavior of a Fermi liquid is observed for all reactions but at different temperatures between 5 and 7 MeV. The slopes are somewhat different and the upbend at high excitations is only seen in the 197Au spectator decay (Fig. 5). It is not observed for the other reactions which, however, do not easily lead to comparable excitation energies (Fig. 3).

Apart from the experimental differences and imperfections, it is the transient nature of the reaction process which most likely prevents a single universal curve to emerge from these studies. To the extent that the equilibrated breakup state is an idealization, measured fragment yields represent integrals over finite emission times , with pre-breakup and post-breakup contributions varying among different reactions. Furthermore, the expansion dynamics may generate transient pressures that are different for different types of hot-fragmentation reactions. The prominent role of the pressure, however, is known from the model experiments .

## 4 PARAMETERS OF THE BREAKUP STATE

Temperatures and excitation energies are only two of the quantities of interest that characterize the breakup states. Techniques for determining other breakup parameters have been successfully developed and applied. Among them, the density is of particular importance because an expansion to low density is a basic ingredient of the multifragmentation scenario.

Low densities in agreement with model expectations were recently reported for spectator decays following 197Au on 197Au reactions at 1000 MeV per nucleon . They were deduced from measured correlation functions for proton pairs and for unlike pairs of protons or deuterons in coincidence with $`\alpha `$ particles. The correlation functions were found to exhibit the surprising property that their variation with $`Z_{\mathrm{bound}}`$ is rather small, indicating source extensions that do not change dramatically with impact parameter in these reactions. The observed variation of the density with impact parameter is mainly caused by the variation of the mass of the intermediate spectator system (Fig. 6).

Associated time scales were deduced from the same proton-proton data with the technique of directional analysis . They are rather short, of the order of 20 fm/c, and comparable with the collision time in the entrance channel. While this may indicate that a majority of these protons, selected to have energies of $`E`$ 20 MeV in the spectator frame, may originate from early stages of the collision, it is still obvious that hot fragmentation initiated by relativistic projectiles is a fast process. The emission times obtained from fragment-fragment correlations (Fig. 2) are larger by only a factor of two.

With experimental values for temperatures and densities at hand, the breakup conditions can be displayed in the temperature-versus-density plane, commonly used to characterize the nuclear phase behavior. The three data points (errors are omitted) shown in Fig. 7 are obtained by correlating the measured temperatures $`T_{\mathrm{HeLi}}`$ (Figs. 3 and 5) with the p-p densities of Fig. 6. The resulting trajectory in the phase diagram tends to higher temperatures and lower densities in the more violent collisions. It is important to note that this breakup (equilibrium) trajectory is practically orthogonal to the actual reaction trajectory followed by the system as it expands from higher to lower temperatures .

The measured breakup conditions are at comparable densities but at temperatures far below the critical point of nuclear matter as obtained in the calculations of Refs. . They are also below the critical temperature obtained for the finite 16O nucleus from the FMD experiment . The conclusion that the breakup occurs in the coexistence region thus seems justified also from this perspective. Alternatively, it may be derived from the observation that multifragmentation populates the partition space predicted for the coexistence region by the statistical models.

The answer to the question what causes fragmentation is less obvious. As sketched in the figure, the experimental breakup points are barely inside the adiabatic spinodal of the infinite system. In finite nuclei, this region of volume instability is probably limited to much lower temperatures. The expanding systems therefore seem to fragment before coming even close to it. Such a conclusion may appear speculative at present, but the need for systematic studies of the instabilities of finite nuclei is rather obvious . Surface modes are an alternative to bulk instabilities. While they are expected to be slow in homogeneous systems , they may be rapidly excited during the early stages of the collision. Nuclear systems are predicted to be resilient to gentle surface excitations but not to major distortions following more violent encounters .

Additional insight in this direction may come from an improved understanding of the kinetic energy spectra of light particles and fragments. In cases of bulk fragmentation, the slope temperatures extracted from fragment spectra are inconsistent with the chemical temperatures obtained with the double-ratio method. However, for spectator decays at relativistic bombarding energies, it has been shown that these temperature values can be consistently understood if the slope temperatures are assumed to reflect the intrinsic Fermi motion, as assumed in the Goldhaber model . Recent calculations with transport models, which incorporate Fermi motion, support this interpretation . The energies of spectator fragments are well reproduced, and the coexistence of qualitatively different internal (or local) temperatures and fragment slope temperatures has been demonstrated (Fig. 8). The experimental and theoretical findings suggest that fragments are preformed at an early stage in these collisions ($``$ 50 fm/c) before they may expand to typical breakup densities. To resolve this apparent contradiction to the filling of the phase space at breakup, the basis of the statistical approaches, is an interesting task for the future.

## 5 CONCLUDING REMARKS

There are many different ways for transforming cold into hot nuclear matter, and we have learned to exploit their different features in our efforts to understand the fragmentation processes. The question of how fragments are formed, the problem of identifying the dominant mechanisms, has been a recurrent theme during this talk. It continues to be a challenge, even though remarkable progress has been made within the last few years.

While this requires a realistic modelling of the dynamics, it does not reduce the statistical approaches in their role and importance. Not only is it demonstrated that the instabilities governing the fragmentation fill the phase space, but it is also the statistics and thermodynamics of the process that allow us to establish the connection to the nuclear phase transition.

The liquid-gas phase transition continues to act as a major motivation and has inspired wider investigations of phase transitions in small systems, extending beyond the nuclear domain. This more general approach seems extremely interesting and rewarding. The caloric curve of nuclei has been emphasized in this talk, but there are other potential observables, not necessarily less challenging on the experimental side, which will be discussed at this conference and should be further pursued and exploited in the future.

The uncertainties associated with the measurement of temperatures have been extensively investigated. The temperature differences seen with different methods have led to new insights into the fragmentation mechanism. Experimentally, the bigger challenge probably is the identification of the source and the measurement of its excitation energy at breakup. There will be limitations in how far we can go with the assumption of a well defined breakup configuration. The continuous evolution of the reaction and emission processes may not allow a precise distinction between the equilibrated emissions and those preceding it.

Finally, the majority of the data included here has come from elaborate experiments with approximately 4-$`\pi `$ coverage in the respective source frames. It is a pleasure to just look at such data, and it is gratifying to see the big investments in funds and manpower being justified. Moreover, it provides encouragement for those who are presently designing or completing potentially even more powerful detection devices for future research.

Illuminating and inspiring discussions with my colleagues at the GSI, notably A.S. Botvina, H. Feldmeier, U. Lynen, W.F.J. Müller, W. Nörenberg, W. Reisdorf and C. Schwarz are gratefully acknowledged. In addition, I like to thank L. Beaulieu, D. Durand, T. Gaitanos, J. Łukasik, K.H. Schmidt and V.E. Viola for providing me with graphics for the talk and for this manuscript.


format: text

text:
# Untitled Document

Particle propagation upstream of CIR shocks

M. Savopulos
Information Dept, Education and Technological Institute of Thessaloniki, Greece
J. J. Quenby and M.K. Joshi
Physics Department, Imperial College, London, SW7 2BZ, UK
M. Fraenz
Queen Mary and Westfield College, Astronomy Unit, Mile End Rd., London E1 4NS, UK.
Abstract The first high solar latitude pass of the Ulysses spacecraft revealed the presence of MeV particle increases up to latitudes above 60 degrees, well outside the CIR belt but associated in time with the regular passage of these plasma interfaces at more equatorial latitudes. The particle increases have been explained variously as due to diffusion from a field line connection with a CIR at a greater distance, perpendicular diffusion in latitude, propagation from the inner Corona or an acceleration process on a connecting field line. Numerical solutions to the 1-D propagation equation upstream of a CIR shock for reasonable diffusion mean free paths are shown here to limit the source of the increases to within about 2 AU of the CIR’s, asssuming the CIR’s either trap or accelerate energetic particles and they diffuse from the nearby CIR interface. The problem of allowing propagation from the CIR’s is eased if the CIRs are located according to the predictions for the current sheet with a solar source surface at $`2.5R_s`$, rather than at $`3.25R_s`$. Energetic electron observations showing delays with respect to the CIR-associated ions suggest possibly an acceleration in the inner corona, rather than in the CIR-associated shock. A more likely explanation, however, is the ability of electrons to take a longer a longer route, from beyond the point of observation, than for that of the ions and yet arrive within a few days of CIR closest approach time.
1. Introduction
In the course of the first Ulysses southern latitude polar pass it was found that whereas the interaction regions between high speed and low speed solar wind driving CIR’s were limited to solar latitudes between $`13^{}`$ and $`40^{}`$, associated energetic particle increases appeared up to $`64^{}`$ in the case of 0.5-1.0 MeV protons and to $`75^{}`$ in the case of 50 keV electrons. Despite the well known longitudinal spreading of such particle increases in equatorial regions (O’Gallagher and Simpson, 1966 ), this observation constituted a major surprise to those working on interplanetary energetic particle transport. It clearly provides crucial information on the relative importance of various competing plasma and energetic particle transport processes (Fisk, 1997). The particles may have originated in the CIR’s due to interplanetary diffusive shock acceleration (Palmer and Gosling, 1978), or they may have been mainly accelerated close to the sun and carried out in a trapping region (Lim et al., 1996). In either case, cross-field diffusion from the CIR to Ulysses (Kota and Jokipii , 1995) or sunward diffusion along a field line connecting to an expanded CIR beyond the Ulysses orbit (Keppler et al., 1995, Lanzerotti et al., 1995) or sun-ward diffusion along a field line constrained to move in latitude under the combined influence of differential solar rotation and differing wind speeds which eventually meet the streamer belt ( Fisk, 1996) are all possible modes of propagation for the high latitude energetic particles. A third alternative is that particles, originally accelerated at the sun, diffuse in the inner corona before they move out along connecting field lines to Ulysses ( Quenby et al., 1996). In the following we describe data relevant to the high latitude increases and discuss possible models.
An additional complication in the observational data lies in the delayed onset of the CIR-related energetic electron increases, relative to ions, above the low latitude streamer belt, Simnett and Roelof (1995).
2. Data
The helium data to be employed has already been presented by Keppler et al., (1995) and was obtained from the Ulysses EPAC experiment. 4 identical, three-element semiconductor telescopes mounted at different angles to the spin axis yielded 8 sector information over $`80\%`$ of the solid angle The geometry factor to measure 0.3 to 1.5 MeV protons and 0.4 to 6 MeV/N heavy ions was 0.08 cm2ster.
We confine ourselves to the use of 0.4-1.0 MeV/N omnidirectional helium data obtained on the first Ulysses southern solar pass. Figure 1 is reproduced from Keppler et al., (1995) in order to clarify subsequent discussion. The regular series of energetic alpha particle increases, 26 days apart are shown plotted against time, spacecraft-sun distance and heliolatitude. In general, the greater part of these particle increases lie between well-identified forward-reverse shock pairs in the low latitude ’streamer belt’ but cease to be enclosed by significant plasma discontinuities at high latitudes (Phillips et al, 1995). $`40^{}`$ is taken as the transition point between ’low’ and ’high’ in our subsequent analysis. A mean spectrum for energetic He ions within a CIR is established from data obtained close to the reverse shocks identified between latitudes -32 and-38. We will also use data at higher latitudes, again coinciding with the regular particle increases which seem CIR-associated, at the distance (AU) and southern latitude (degrees) pairs; (4.2,41),(4.0,44),(3.8,47),(3.6,52),(3.4,56),(3.2,60).
Electron and ion data are obtained from the Ulysses HI-Scale experiment, (Lanzerotti et al., 1992). We use the archived data from the thin, single solid-state detector system with active anti-coincidence and thin foil separation of electrons in the energy range 30 keV to 300 keV and ions in the energy range 300 keV to 5 MeV. The telescopes involved have solid angles of 0.48 cm-2 sr.
3.Shock Expansion Model
Since it will turn out that the models for the high latitude increases most stringently constrained by the above data are those involving su-ward diffusion from a CIR source located well beyond Ulysses, we will centre our numerical analysis on one of these, namely the expanded CIR beyond the spacecraft orbit. In order to determine the likely range over which diffusion takes place, it is thus necessary to establish plausible estimates of the expansion history of a CIR as it moves beyond Ulysses. Gosling et al., (1993) pointed out that the reverse shock moved poleward. Figure 2 is a cartoon representing a meridian plane cross-section of an interaction region in the southern hemisphere with the equator-ward forward shock expansion and the pol-ward reverse shock expansion. Ulysses is located at a latitude such that it does not directly see the CIR interaction region but lies on a radial line from the sun which subsequently intersects the expanded CIR. A simple estimate of the expansion is obtained following Quenby et al., (1995), who assumed a radial speed of 500 km/s and found the pole-ward or latitude expansion speed arising from the excess pressure within the CIR to be the Alfven speed of 150 km/s. In this calculation, it is further assumed that the extent of the streamer belt remains constant in it’s upper latitude boundary at a given radial distance throughout the southern latitude pass. Starting from the observed extent of the streamer belt as reaching 36S at 4.4 AU, the CIR is then found to lie on the same radial vector as Ulysses at the following pairs of distances and latitudes; (5.8,41),(6.9,44),(8.2,47),(10.8,52), (13.6,56),(17.1,60). This calculation of the CIR expansion is to be regarded as the maximum possible and indeed, the simulations by Pizzo, (1991) suggest a reduction in the rate of latitude increase beyond $``$ 10 AU.
It is not easily established how the intensity increase inside the CIR varies with position. We have tried correlations with radial distance, Ulysses latitude and the maximum extent of the neutral sheet at the time. There is only a poor correlation with the later two parameters but while the radial increase yields the dominant correlation, we cannot from our data and epoch extrapolate beyond 5 AU and so we take a favourable case and assume that the intensity increase stays constant. Thus we take the mean intensity found previously as applying both as the source intensity inside the CIR for the (4.2,41) observation and to the subsequent points where Ulysses is connected on a solar radius vector to the CIR which has propagated beyond and above this position. Infact, if the Van Hollebeke et al., (1978) data on the equatorial plane radial dependence of the peak CIR particle intensity applies, the source flux within the CIR will diminish beyond 4-5 AU. All the intensities corresponding to the above pairs of positions are assumed to refer to observations upstream (or sun-ward) of the expanded reverse shock of the CIR. These intensities, normalised to the intensity inside the CIR (shock and observation at 5.8 AU) are plotted in figures 3 and 4, together with the estimated upstream distances on the Quenby et al., (1996) model.
If the connection between the CIR from beyond Ulysses to the spacecraft is via a field line moving in latitude according to the Fisk, (1996) model, the radial distance to be moved is typically 10 AU for a line seen at 70 in the inner heliosphere. Thus calculations on the flux expected upstream of shock on the expanded CIR model are applicable to geometries based upon the Fisk, (1996) model provided in each case, the dominant mode of particle diffusion is radial.
4.Semi-quantitative consideration of the possible propagation
The ability of MeV particles to propagate $`10`$ AU to be seen with the observed intensity and to be reasonably in phase with the established low latitude solar rotation periodicity are critically dependent on the radial or perpendicular diffusion mean free path adopted. It is difficult to better the ’Palmer Consensus’ value of 0.1 AU (Palmer, 1982) as applying in a wide rigidity and distance range for $`\lambda _{rr}`$. However a detailed study of the time variability in the magnetic fluctuation spectum on an hour to hour basis and of the anisotropy injection profiles of long lived particle events by Wanner and Wibberenz, (1991) found $`\lambda _{}`$ varying between 0.01 and 1.0 AU. Reames (1999) maintains that both prompt and gradual SEP events are consistent with $`\lambda _{}`$ 1 AU. The ’Palmer’ value is in accord with recent realistic, magnetometer based computations of the parallel mean free path (Drolias et al., 1997) over the Ulysses orbit which yield $`\lambda _{}0.1`$ AU provided we are at $`1`$ AU. There is, however, already a puzzle beyond $`5`$ AU in adopting this relatively low value of $`\lambda _{}`$ where $`\lambda _{rr}=\lambda _{}cos^2\chi `$ where $`\chi 70^{}`$ or $`\lambda _{rr}0.01`$ AU , unless $`\lambda _{}\lambda _{}`$.
Using the expression for the diffusion time, $`\tau `$, to reach peak flux over distance, L, with velocity, v, $`\tau (3/4)L^2/\lambda _{rr}v`$, we find it takes 94 days to diffuse 10 AU, 23 days to diffuse 5 AU and 90 hours to diffuse 2 AU. These numbers suggest that unless the ’source’ of the high latitude increases is within a few AU of Ulysses, the 27-day periodicity in the enhancements will be difficult to explain.
5.Quantitative radial back diffusion model
In this section, we numerically solve the Fokker-Planck transport equation for the case of spherically symmetric propagation back towards the sun from the reverse shock of a CIR. The Fokker-Planck with only radial coordinate dependence and assuming a steady state is

$$\frac{1}{r^2}\frac{}{r}(r^2k\frac{f}{r})(VV_s)\frac{f}{r}+\frac{2}{3}\frac{Vv}{r}\frac{f}{v}=0$$
(1)
Here k is the radial diffusion coefficient where $`k=\lambda _{rr}v/3`$ , f the distribution function, V the wind speed, $`V_s`$ the shock speed, r the solar radial distance and all quantities measured in a solar reference frame. The equation expresses the facts that while the diffusive and adiabatic deceleration divergence terms depend on distance from the Sun, the convective motion must be considered relative to the shock. The boundary condition derived from continuity of streaming at the CIR reverse shock is

$$k(\frac{f}{r})_1+\frac{v}{3}(\frac{f}{v})_1(V_sV_1)=\frac{v}{3}(\frac{f}{v})_2(V_sV_2)$$
(2)
where 1 and 2 refer to upstream and downstream conditions. Using the inverse compression ratio $`\beta `$, $`V_sV_2=\beta (V_sV_1)`$. The boundary condition equates upstream diffusive and convective flows with downstream convection with negligible diffusion. Following Savopulos et al., (1995), equations 1 and 2 are combined and numerically solved by inverting a quindiagonal matrix allowing $`k=k_{}v^\alpha r^ϵ`$ . $`k_{}`$ will be quoted with v measured in units of $`10^8`$ m s-1 and r in units of $`10^{11}`$ m.
We show in figures 3 and 4 as continuous lines solutions to equations 1 and 2 for two relations for k. In figure 3, $`\alpha =ϵ=1`$ and $`\lambda _{rr}`$=0.008AU at 1 MeV/N and 1AU. In figure 4, $`\alpha =1.6,ϵ=0.0`$ and the $`k_{}`$ value of 5.1435 in the above units corresponds to $`\lambda _{rr}=0.44`$ AU at 1MeV/N independent of distance. The top curve is the fit to the He spectrum inside the CIR while the others correspond to the computed upstream spectra for the assumed shock-Ulysses distances on the model of section 3 for the latitudinal spread of a CIR. These distances can similarly be employed in the model of the latitudinal wandering of an interplanetary field line. The figures show that a variety of mean free path models can fit the data (as demonstrated by Savopulos and Quenby, 1998) although independent propagation evidence clearly favours an intermediate mean free path value. However, both figures confirm that there is little prospect of fitting the intensity fall-off apparently observed upstream of the reverse shock except within about 1 or 2 AU of the CIR, if either back diffusion model is adopted.
The fall in intensity due to convective sweeping and adiabatic expansion over $``$ 5AU also calls into question models requiring coronal plus interplanetary propagation.
Some significant reduction in the requirement on the distance travelled in latitude is gained by following Sanderson et al.,(1999) who employ the Wilcox Solar Observatory curent sheet locations based on the $`2.5R_s`$ source surface model, rather than that with the source surface at $`3.25R_s`$, which is often used. Both models assume potential fields and a radial source surface field, but the former adds a sharply peaked polar field to the surface line-of-sight data while the latter assumes a purely radial surface field. Both models are reasonably, but not entirely successful in predicting to the IMF (eg Hoeksema,1995). Using the $`2.5R_s`$ model the curent sheet is found to have passed within about $`5^{}`$ of Ulysses up to about day 60, 1994, although there still remains a gap opening to $`30^{}`$ at the last ion increase seen around day 160, 1994, so the problem is only possibly partly solved with the model with the lower source surface. The active interaction region thought to be operative at these times is not considered to extend much south of the maximum southward extent of the neutral sheet at the Ulysses radius.
6. Interpretation Of ’Delayed Electrons’ Figures 5a and 5b illustrate the delayed electron arrival relative to the ions above the streamer belt. Figure 5a runs from day 183 to day 273, 1993 and shows the electron channels 30-50 keV(E1), 50-90 keV(E2), 90-165 keV(E3), 165-300 keV(E4) and the ion channels 300-550keV (FP4), 550-keV-1MeV (FP6) and 1-5MeV (FP6). The electron and ion increases up to day 225, corresponding to a South latitude of about $`37^{}`$, tend to peak at the same time, although the electron increase often extends to later times or further upstream from the reverse shock, behaviour confirmed by similar plots at lower latitudes. By day 261, corresponding to a latitude of about $`41^{}`$, the electron increase is delayed 2 days in the lowest energy channel at least. Figure 5b for the time period day 274 to 365, 1993 shows a delay of about 2 days developing at all energies shown, again a number confirmed by early 1994 data.
In seeking to understand the electron data delay we note the conclusion of Palmer (1982) that the parallel electron mean free path between 10 keV and 10 MeV lies in the interval 0.1-1 AU with a likely value $`0.1`$U AU. Scatter-free events with $`\lambda 1`$ AU are judged to be relatively rare. Assuming perpendicular diffusion over the shortest possible mean free path, let us attempt an interpretation of the extra electron time delay around day 343, 1994, when Ulysses was at $`3.9`$ AU at latitude $`46^{}`$ south. According to Sanderson et al,(1999), the current sheet was at $`40^{}`$ when it passed closest to Ulysses. Hence the distance over which perpendicular diffusion is to take place from the CIR to Ulysses is $`0.4`$ AU. A way of estimating the MeV proton propagation time is to take a $`\lambda _{}0.1\lambda _{}`$. This is within an acceptable range and is favourable to explaining the electron data by the same propagation mode. With$`\lambda _{}0.1`$AU, $`\lambda _{}`$=0.01AU. This yields a 1.7 day delay. Suppose a total diffusion time across the field above the CIR of 1.7 day (suggested by the proton diffusion data) plus the 2 days extra electron delay at 230 keV. Then we require an electron $`\lambda _{}=0.00030.0004`$ AU in the range 40-230 keV. There is thus no consistency with the ’Palmer consensus’ value of the cross field $`k_{}/\beta 10^{21}cm^2s^1`$ or $`\lambda _{}7.10^3`$AU. Neither is the theoretical field line wandering estimate (Forman et al.,1974) of $`k_{}/\beta 4\times 10^{20}cm^2s^1`$ satisfied, especially the velocity independence of this quantity.
Computations by Kota and Jokippi, (1998), suggest that the 100keV electrons, in contrast to the MeV ions, can infact arrive from a CIR position further out than Ulysses. These authors employ a tilted dipole model of the heliospheric field, with purely radial flow, that incorporates the forward and reverse shocks at the interfaces between the slow and fast streams. The 3-D test particle transport equation is solved, including already accelerated electron and ion spectra injected at the shock fronts. The simulations reveal a broader time peak in the high latitude arrival with several days delay with respect to the ions, suggesting propagation predominantly from further out.
Using our ’semi-quantitative approach’ as before, we estimate a three day delay for 100 keV electrons with $`\lambda _{}`$ =0.1 AU would enable a distance of 1.8 AU to be traversed along a field line whereas this distance increases to 5.7 AU if $`\lambda _{}`$ =1.0 AU. Since these relativistic electrons are in a propagation parameter parameter regime where diffusion dominates convection in the transport equation (1), the Fisk and Lee (1980) approximation to the remaining, diffusion and adiabatic deceleration terms should yield the radial intensity dependence upstream of the shock. This solution gives $`f=F(v)(r/r_s)^{2\beta /(1\beta )}`$ provided $`k=k_{}vr`$ with $`r_{}`$ as shock location. Hence with $`\beta =1/4,fr^{2/3}`$ and therefore the intensity reduction over 2-6 AU is not drastic. This fact can also be inferred from fig 4 which in the high velocity regime and for $`\lambda _{}`$ =0.44 AU corresponds to the diffusion dominated solution regime. Note that the theoretical curves approach each other at the highest velocities, denoting a low spatial gradient in the parameter space where the radial dependence is expected to be particle velocity and species independent according to the Fisk and Lee (1980) approximation. Clearly with an appropriate choice of the $`k_{}/k_{}`$ ratio, one can arrange for the cross field diffusion electron intensity to be less than the along the field diffusion electron intensity. Note this 1-D discussion of an essentially 3D situation depends on $`k_{}<<k_{}`$.
An alternative approach to the electron delay is to assume most acceleration takes place close to the sun within the streamer belt confines or with particle release mainly situated within this region. It is unlikely that electrons appear at higher latitudes relative to protons because of a basic asymmetry in the output of reconnection event aceleration because such events are orientated at random. In the IMF at the epoch of observation, the gradient and curvature drifts will produce a drift velocity in a direction such that the electrons appear later, but the drift speed is only is only $`3\times 10^5`$ cm/s. Hence this is unlikely to be a cause of the measured delay as insufficient latitude separation occurs. A possibility is that at $``$ 100 keV or $``$ 0.3MV, $`\lambda _{}`$1 AU. Then electrons accelerated within a few $`R_s`$ propagate easily to Ulysses, outside the CIR belt, after coronal diffusion while protons $``$0.4 MeV or $``$30 MV are much attenuated in intensity in the coronal diffusion region or even in the IMF. This idea presupposes that favourable acceleration conditions involving reconnection or shock propagation across closed magnetic loops occur mainly within the streamer belt. The observation by Van Hollebeke et al, (1978), that within CIR’s the proton intensity increases more than 100% from 0.4 AU to 1 AU simply states that some interplanetary acceleration occurs within the CIR but does not invalidate the idea of an important coronal component to acceleration. A coronal diffusion model explains the reduction in appearence rate with increased latitude for electrons noticed by Lanzerotti et al, (1995). CIR particle abundances are mainly similar to solar wind abundances (Keppler, 1998), rather than flare abundances. However, the likely source of the solar wind is small-scale reconnection events, so there can be no basic incompatibility in the source abundances. Large-scale flares are clearly potentially different in plasma constitution.
7.Conclusions
With accepted values of the radial mean free path in the few AU, few MeV/N region, the models for high solar latitude CIR associated particle increases which involve diffusion over 5-10 AU, either from behind Ulysses or directly from the Sun, have difficulty in fitting the relatively high intensities seen up to 60 when spherically symmetric solutions of the Fokker-Planck equation are employed. The problem is enhanced by the fact that favourable assumptions were made for the intensity of the source particles within CIR’s beyond Ulysses and for the height of the CIR expansion in latitude. There is an additional problem in the long time required for such diffusion to take place rendering the observed phase coherence of the particle increases difficult to understand. An alternative is to consider a perpendicular diffusion mean free path not much smaller than the parallel $`\lambda `$, allowing the CIR associated fluxes to propagate only 1-2 AU to locations directly above the shock in latitude. Use of the neutral sheet model with solar source at $`2.5R_s`$ reduces the amount of latitude propagtion required although it remains significant.
The relatively large value of the perpendicular diffusion coefficient required in the above analysis fits with the relatively small latitudinal gradients in low energy cosmic ray intensity seen by Ulysses which appear significantly less than expected on full drift models with $`K_{}=0.05K_{}`$ (see Drolias et al., 1997 and references therein).
100 keV electron increases observed above the streamer belt are not easily explained on the purely perpendicular diffusion transport model with reasonable parameters. Instead, it is possible that a coronal source is the basic accelerator and the scattering mean free path outside CIR’s at sub MV rigidities is of the order of 1 AU. Most ion transport would then be effected by trapping within the CIR’s while most electron transport would be first via the corona and then in the quiet IMF. However, back diffusion from the expanding CIR beyond Ulysses is a likely explanation for the delayed and extended electron arrival as fewer special assumptions arise in this last model.
8.Acknowledgements
Helpful discusions with Bern Blake, Erhard Keppler, Bob Forsyth, Trevor Sanderson and Bernd Heber are acknowledged. PPARC visitors’ support for M. Savopulos was gratefully received.
9.References
Drolias,B.,Garaud,P.,Quenby,J.J. and Smith,E.J.:1997, Proc.25th. ICRC,1,229.
Fisk,L.A.:1996, JGR,101,A7,15,547.
Fisk,L.A.:1997, Proc.25th. ICRC,8,27.
Fisk,L.A and Lee,M.A.:1980, Astrophys. J. 237,620.
Forman,M.A.,Jokippi,J.R. and Owens.:A.J.,1974, Astrophys.J.,192,535.
Gosling,J.T.,Baker,D.N.,Bame,S.J.,Feldman,W.C.,Zwickel.R.D.and Smith,E.J.:1993, Geophys. Res. Lett.,20, 2789.
Hoeksema, J.T.: 1995, Space Sci., Rev.,72/1-2,137.
Keppler,E.,Franz,M.,Korth,A.,Reuss, M.K.,Blake,J.B.,Seidel,R,Quenby,J.J. and
Witte,M.:1995,Science,268,1013.
Keppler,E.:1998, Surveys in Geophysics, 19,(3 SIS),211.
Kota,J. and Jokipii,J.R.:Science, 268(5213), 1024.
Kota,J. and Jokipii,J.R.:1998, Space Sci.Rev.,83,137.
Lanzerotti,L.J.,Gold,R.E.,Anderson,K.A.,Armstrong,T.P.,Lin,R.P.,Krimigis,S.M.,Pick,M,. Roelof,E.C.,Sarris,E.T., Simnett,G.M. and Frain.,W.E.:1992, Astron. Astro. Phys. Suppl. Series,92,349.
Lanzerotti,L.J.,Armstrong,T.P.,Gold,R.E.,Maclennan,C.E.,Roelof,E.C.,Simnett,G.M.,
Thomson,D.J.,Anderson,K.A., Hawkins III,S.E.,Krimigis,S.M.,Lin,R.P.,Pick,M.,Sarris,E.T. and
Tappin,S.J.:1995,Science,268,1010.
Lim,T.L,Quenby,J.J.,Reuss,M.K.,Keppler,E.,Kunow,H.,Heber,B. and Forsyth,R.J.:1996, Ann. Geophys.-Atmos. Hydrospheres Space Sci. 14(4),400.
O’Gallagher J.T. and Simpson, J.A.: 1966, Phys. Rev. Lett., 16, 1212,
Palmer,I.D. and Gosling,J.T.:1978,JGR, 83,2037.
Palmer,I.D.:1982, Revs. Geophys. Space Physics, 20, 335
Phillips,J.L., Bame,S.J.,Feldman,W.C.,Goldstein,B.E.,Gosling,J.T.,Hammond,
C.M.,McComas,D.J.,Neubauer,M., Scime,E.E.,and Suess,S.T.,1995.: Science, 268, 1030.
Pizzo,V.J.,1991.: J. Geophys. Res., 96, 5405.
Quenby,J.J.,Drolias,B.,Keppler,E.,et al.: 1995, Geophys. Res. Lett.,22,3345.
Quenby,J.J.,Witcombe,A.,Drolias,B.,Fraenz,M. and Keppler,E.:1996,Astron. Astrophys. 316,506.
Reames,D.V.:1999, Space Sci. Rev.,90,413.
Sanderson,T.R.,Lario,D.,Maksimovic,M.,Marsden,R.G.,Tranquille,C.,Balogh,A.,Forsyth,R.J. and Goldstein,B.E.: 1999, Geophys. Res. Lett.,26,1785.
Savopulos,M.,Quenby,J.J. and Bell,A.R.:1995, Solar Phys.,157,349.
Savopulos,M. and Quenby,J.J.:1998, Solar Phys.,180,479
Simnett, G.M. and Roelof,E.C.:1995, Space Science Reviews,72,308
Van Hollebeke, M,A,I., McDonald,F.B.,Trainor,J.H. and Von Rosenvinge,T.T.: 1978, JGR, 83, 4723,
Wanner,W. and Wibberenzz,G.: 1991, 22nd Int. Cosmic Ray Conf., Dublin, 3,221.
Figure Captions
Figure 1 Overview of the first Ulysses south high latitude pass showing the solar wind, magnetic field strength and energetic helium flux as a function of heliolatitude.
Figure 2 Representation of a heliospheric meridian plane showing the oscillating neutral current sheet, the developing forward and reverse shocks and the radial line through the Ulysses spacecraft.
Figure 3 Observed high latitude He intensities as a function of particle velocity measured in MeV/N units. Positions of observation and modelled positions where the expanding CIR’s reach the radial line from the sun through Ulysses beyond the spacecraft are indicated. Computational solutions corresponding to these upstream of the shock positions are shown as solid lines. The diffusion parameters employed are stated at the bottom.
Figure 4 Similar to figure 3 for alternative diffusion parameters.
Figure 5 Ulysses electron CIR observations above the streamer belt, 1993, days 183-273. Comparison ion channels are also shown.
Figure 6 Ulysses electron CIR observations above the streamer belt 1993, days 274 to 365. Comparison ion channels are also show.


text:
# Superluminal Recession Velocities

## Two Kinds of Velocity

Despite the efforts of Murdoch C:tdavis:murdoch77 , Harrison C:tdavis:harrison81 ; C:tdavis:harrison91 ; C:tdavis:harrison93 , Stuckey C:tdavis:stuckey92 , Ellis & Rothman C:tdavis:ellis93  and Kiang C:tdavis:kiang97  there is much confusion about superluminal recession in the literature C:tdavis:davis00 . Special relativistic calculations are often incorrectly invoked in a cosmological context and recession velocities with an upper limit of c are obtained. A common, but mistaken line of reasoning goes as follows: 1. The most distant object we can observe has an infinite redshift. 2. An infinite redshift implies a velocity equal to the speed of light. 3. Therefore the most distant object we can see has a velocity of $`c`$. This logic falters because the second premise is invalid: (2.) would only be true if the cosmological redshift were a special relativistic rather than a general relativistic effect. Since redshifts and recession velocities are central to the idea of an expanding universe, the distinction between recession and peculiar velocity is fundamental. We define the velocities and redshifts mentioned above and present spacetime diagrams showing how light can reach us from regions of space which are receding superluminally. We have generalised the solutions to arbitrary $`(\mathrm{\Omega }_M,\mathrm{\Omega }_\mathrm{\Lambda })`$.

The Friedmann-Robertson-Walker (FRW) metric C:tdavis:peacock99  describes an homogeneous, isotropic universe. We are interested in only the radial components: $`ds^2=c^2dt^2+a(t)^2d\chi ^2`$. Radial distance along a constant time slice ($`dt=0`$) is given by: $`D=a\chi `$. Differentiating this physical or proper distance gives a two component velocity:

$`v_{tot}`$ $`=`$ $`\dot{a}\chi +a\dot{\chi }`$ (1)
$`=`$ $`v_{rec}+v_{pec}`$
Recession velocity ($`v_{rec}`$) is the velocity in Hubble’s law since Hubble’s constant $`H=\dot{a}/a`$. Peculiar velocity ($`v_{pec}`$) is unassociated with the expansion of the universe and corresponds more closely to common usage of the word ‘velocity’. For example setting $`ds=0`$ for a photon in the FRW metric gives a radial peculiar velocity $`a\dot{\chi }_\gamma =c`$. The conversion from cosmological redshift to recession velocity is different from the conversion from an SR Doppler shift to peculiar velocity (Fig. 1):

$`v_{pec}(z)=`$ $`c{\displaystyle \frac{(1+z)^21}{(1+z)^2+1}}`$  SR Doppler shift (2)
$`v_{rec}(z)=`$ $`\dot{a}(z){\displaystyle _o^z}{\displaystyle \frac{cdz^{}}{H(z^{})}}`$ GR cosmological redshift (3)
Equation (3) converts high redshifts into superluminal recession velocities but does not contradict special relativity because nothing ever overtakes a photon and all observers measure local photons to be traveling at $`c`$.

In Figure 2 the boundary between the light and dark grey regions is the Hubble sphere, the distance, $`D_{HS}(t)=c/H(t)`$, at which galaxies are receding at the speed of light. The comoving distance to the Hubble sphere increases when the universe decelerates and decreases when the universe accelerates. The Hubble sphere is not an horizon of any kind; it passes over particles and photons in both directions.

How can photons move from the light grey superluminal region into the dark grey subluminal region? How can a swimmer make headway against a current that is faster than she can swim? It seems impossible. First consider the problem in comoving coordinates. Photons propagate towards us along our trumpet-shaped past lightcone with a comoving velocity $`\dot{\chi }_\gamma =c/a`$. Therefore, photons move inexorably through comoving space, irrespective of their position in relation to our Hubble sphere. Their comoving velocity is always negative. Their comoving distance always decreases.

Now consider the problem in physical space. The tear-drop shape of our past lightcone means that the distance between us and the most distant photons we now see, was once increasing. The relevant quantity for understanding this behaviour is the total velocity of a photon that is heading towards us: $`v_{tot}=v_{rec}c=HDc=\dot{a}\chi _\gamma c`$. The total velocity of distant photons is not constant because it is the sum of the distance-dependent recession velocity ($`v_{rec}`$) and the constant peculiar velocity, $`c`$. When $`\dot{a}\chi _\gamma >c`$ the distance between us and the photon increases.

Photons now reaching us from our particle horizon were emitted at the Big Bang. Since $`lim_{t0}\dot{a}=\mathrm{}`$ in all viable cosmological models, $`\dot{a}\chi _\gamma >c`$ would have been satisfied and these photons would have inially receded from us. Similarly, the first photon we receive from any object was emitted from a region with $`\dot{a}\chi _\gamma >c`$, $`v_{tot}>0`$. In the early universe as time progresses, both $`\chi _\gamma `$ and $`\dot{a}`$ decrease. Thus, $`v_{tot}`$ of these photons evolves from positive to negative, and the teardrop shape of our physical past lightcone is ubiquitous to all times.

The swimming analogy fails because, unlike recession velocity which is smaller at smaller comoving distances, the current the swimmer has to face is the same at all comoving distances. Our swimmer has to battle an unrelenting cur rent, while the photon constantly moves into regions with a slower “current” (slower $`v_{rec}`$).


text:
# Pressure determinations for incompressible fluids and magnetofluids

## 1 Introduction

It has long been the case that pressure determinations for incompressible flows, both Navier-Stokes and magnetohydrodynamic (MHD), are known to be highly non-local. Taking the divergence of the equation of motion

$$\frac{𝐯}{t}+𝐯𝐯=\frac{𝐣\times 𝐁}{\rho c}P+\nu ^2𝐯,$$
(1)
and using $`𝐯=0`$ leaves us with a Poisson equation for the pressure $`P`$, which is said to function as an equation of state:

$$^2P=(𝐯𝐯\frac{𝐣\times 𝐁}{\rho c})$$
(2)
Here, $`𝐯=𝐯(x,t)`$ is the fluid velocity field as a function of position and time, $`𝐁`$ is the magnetic field, $`𝐣=c\times 𝐁/4\pi `$ is the electric current density, $`c`$ is the speed of light, $`\nu `$ is the kinematic viscosity, assumed spatially uniform and constant, and $`P`$ is the pressure normalized to $`\rho `$ the mass density, also spatially uniform. (1) and (2) are written for MHD. Their Navier-Stokes equivalents can be obtained simply by dropping the terms containing $`𝐁`$ and $`𝐣`$.

If we are to solve (2) for $`P`$, boundary conditions are required. In the immediate neighborhood of a stationary “no-slip” boundary, both the terms on the left of (1) vanish and we are left with the following equation for $`P`$ as a boundary condition:

$$P=\nu ^2𝐯+\frac{𝐣\times 𝐁}{\rho c}$$
(3)
We now focus on the Navier-Stokes case, where the magnetic terms disappear from (3), for simplicity. All the complications of MHD are illustrated by this simpler case. It is apparent that (3) must apply to all components of $`P`$, and that while the normal component of $`P`$ is enough to determine $`P`$ through Neumann boundary conditions, the tangential components of (3) at the wall equally well determine $`P`$ through Dirichlet boundary conditions. This is a problem which some inventive procedures have been proposed to resolve, usually by some degree of “pre-processing” or various dynamical recipes which seem to lead to approximately no-slip velocity fields after a few time steps (e.g. \[Gresho 1991\], \[Roache 1982\] and \[Canuto et al. 1988\] ). It is not our purpose to review or critique these recipes, but rather to focus on a set of velocity fields, related to Chandrasekhar-Reid functions \[Chandrasekhar 1961\], for which (2) is explicitly soluble at a level where the Neumann or Dirichlet conditions can be exactly implemented. In § 2, we explore the difference between the two pressures so arrived at. Then in § 3, we propose a replacement for the long-standing practice of demanding that all components of a solenoidal $`𝐯`$ vanish at material walls, in favor of a replacement by a wall friction term for which the above mathematical difficulty is no longer present. Of course, similar statements and options will apply to all comparable incompressible MHD problems.

## 2 Pressure determinations

Restricting attention at present to the Navier-Stokes case, we consider two-dimen- sional, solenoidal, velocity fields obtained from the following stream function:

$$\psi (x,y)=C_{k\lambda }\mathrm{cos}(kx)[\mathrm{cos}(\lambda y)+A_{k\lambda }\mathrm{cosh}(ky)]$$
(4)
The hyperbolic cosine term in (4) contributes a potential flow velocity component to $`𝐯`$ which makes it possible to demand that $`𝐯`$ obey two boundary conditions: the vanishing of both components at rigid walls \[Chandrasekhar 1961\]. The function in (4) is even in $`x`$ and $`y`$, but can obviously be converted into an odd or mixed one by the appropriate trigonometric substitutions.

The velocity field $`𝐯=\psi \times \widehat{𝐞}_z`$ has only $`x`$ and $`y`$ components and is periodic in $`x`$, with an arbitrary wavenumber $`k`$. $`C_{k\lambda }`$ is a normalizing constant, and the constants $`\lambda `$ and $`A_{k\lambda }`$ can systematically be found numerically to any desired accuracy so that both components of $`𝐯`$ vanish at symmetrically placed no-slip walls at $`y=a`$ and $`y=a`$. In fact, for given $`k`$, an infinite sequence of such pairs of $`\lambda `$ and $`A_{k\lambda }`$ can be determined straightforwardly. Thus any such $`𝐯`$, or superposition thereof, is not only solenoidal, but has both components zero at $`y=\pm a`$, and all spatial derivatives exist. Moreover, the “source” term, or $`(𝐯𝐯)`$, from the right hand side of (2), is of a relatively simple nature for such a $`𝐯`$, since every term in it can be written as a product of exponentials of $`kx`$, $`\lambda y`$ and $`ky`$. It is straightforward to find an inhomogeneous solution for $`P`$, which then is the same for all boundary conditions for a given $`𝐯`$ of the form stated. To this inhomogeneous part of $`P`$ must be added a solution of Laplace’s equation. This can be chosen so that the total $`P`$ may satisfy either the normal component of (3) at the walls, or the tangential component of it, but not both. The determination involves only simple but tedious algebra.

We illustrate, in figure 1, an arrow plot of the velocity field given by choosing $`k=\pi /2`$, $`\lambda =2.6424`$ and $`A_{k\lambda }=.3499`$, in units of $`a=1`$. The two pressures resulting from the satisfaction of the normal and tangential components of (3) can best be compared by comparing their respective values of $`P`$, since $`P`$ itself is indeterminate up to an additive constant in both cases. In figure 2, we display, as an arrow plot, the difference between the pressure gradients associated with the velocity field shown in figure 1. We have rewritten (1)-(3) in dimensionless units for this purpose, with the kinematic viscosity being replaced by the reciprocal of a Reynolds number, which may be defined as $`\text{Re}=(𝐯^2/(k^2+\lambda ^2))^{1/2}/\nu `$. Here, the angle brackets refer to the mean of $`𝐯^2`$ taken over the 2-D box, containing one period in the $`x`$ direction and from $`y=a`$ to $`y=a`$. The value of Re used to construct figure 2 is $`\text{Re}=2293`$, with the dimensionless version of $`C_{k\lambda }=5000`$ in (4). The two pressures are similar but not identical.

In figure 3, a fractional measure of the difference between the “Neumann pressure” $`P_N`$ and the “Dirichlet pressure” $`P_D`$ is exhibited as a contour plot of the scalar ratio

$$\frac{(P_DP_N)^2}{(P_N)^2}$$
(5)
There is no absolute significance to the numerical value of this ratio. It initially increases with Re approaching a maximum of about $`2\%`$ near the wall for $`\text{Re}\mathrm{}>10`$. It is considered interesting however, that the fractional difference is nearly x-independent where it is largest. That occurs formally because the algebra reveals it to be dominated by a term which varies as $`\mathrm{cosh}(4ky)\mathrm{cos}(4kx)`$ in a region where $`ky\mathrm{}>1`$.

It is amusing but perhaps not significant to superpose the velocity field from (4) with a parabolic plane Poiseuille flow of a larger amplitude. The resulting flow field is shown in figure 4, and it bears a striking but perhaps not significant similarity to the flow patterns seen in two-dimensional plane Poiseuille flow \[Jones & Montgomery 1994\] when linear stability thresholds are approached. The pressure gradient difference for this case will be fractionally smaller than in figure 3, since pure parabolic plane Poiseuille flow is a rare case where the two pressures happen to agree, and it quantitatively dominates the pressures determined from equation (2) in this example.

## 3 Discussion and a possible modification

An alternative to the no-slip condition is the “Navier” boundary condition \[Lamb 1932\]: the slip velocity at the wall surface is taken to be proportional to the rate of shear at the wall. This may be expressed $`\mathrm{\Delta }V=L_s\dot{\gamma }`$ where $`\mathrm{\Delta }V`$ is the slip velocity of the fluid at the wall, $`\dot{\gamma }`$ is the rate of shear at the wall and $`L_s`$ is a constant with the dimensions length. Molecular dynamic simulations of Newtonian liquids under shear \[Thompson & Robbins 1990\] have shown this to be the case under some circumstances. In fact recent work \[Thompson & Troian 1997\] has shown that, in cases where the shear rate is large, there is a nonlinear relationship between $`L_s`$ and $`\dot{\gamma }`$.

We note that the velocity field shown in figure 1 does not lead to one which obeys the Navier boundary condition, after an initial time step, where the fluid has been allowed to slip at the wall. If the velocity field determined by (4) is advanced in time using (1) with the “Neumann pressure”, the proportionality between the slip velocity and the rate of shear at the wall, after the initial time step, varies sinusoidally with $`x`$.

It is to be stressed that we are concerned here only with initial conditions, not with circumstances under which initial slip velocities might be coaxed dynamically into vanishing after some time.

It is difficult to see in what sense the velocity field obtained from (4) might be an unacceptable one from the point of view of the Navier-Stokes or MHD descriptions. It seems to have all the properties that are thought to be relevant. The family of functions of the same x-periodicity in (4) can be shown to be orthogonal, and is a candidate for a complete set, in which any $`𝐯`$ might be expanded, when supplemented by flux-bearing functions of y alone. The mathematical question of which if any velocity fields, which are both solenoidal and vanish at the wall, would lead to Neumann and Dirichlet pressures that were in agreement with each other, must remain open. Indeed, the question of whether there are any, without some degree of “pre-processing,” must remain open. This is an unsatisfactory situation for fluid mechanics and MHD, in our opinion, even if it is a not unfamiliar one. The search for alternatives seems mandatory.

One alternative that may be explored is one that seemed some time ago, in a rather different context \[Shan & Montgomery 1994a,b\], to have worked well enough for MHD. Namely, we may think of replacing the requirement of the vanishing of the tangential velocity at a rigid wall with a wall friction term, added to the right hand side of (1), of the form

$$\frac{𝐯}{\tau (𝐱)}$$
(6)
where the coefficient $`1/\tau (𝐱)`$ vanishes in the interior of the fluid and rises sharply to a large positive value near the wall. The region over which it is allowed to rise should be smaller than the characteristic thickness of any boundary layer that it might be intended to resolve, but seems otherwise not particularly restrictive. Such a term provides a mechanism for momentum loss to the wall and constrains the tangential velocity to small values, but does not force it to zero. The Dirichlet boundary condition disappears in favor of a relation that permits the time evolution of the tangential components of $`𝐯`$, while demanding that $`P`$ be determined solely by the Neumann condition (the normal component of (3) only). In a previous MHD application \[Shan & Montgomery 1994a,b\] dealing with rotating MHD fluids, the scheme seemed to perform acceptably well, but was not intensively tested or benchmarked sharply against any of the better understood Navier-Stokes flows. This comparison seems worthy of future attention.

The work of one of us (D.C.M.) was supported by hospitality in the Fluid Dynamics Laboratory at the Eindhoven University of Technology in the Netherlands. A preliminary account of this work was presented orally at a meeting of the American Physical Society \[Kress & Montgomery 1999\].


text:
# 1 INTRODUCTION

## 1 INTRODUCTION

The classical dynamics in polygonal billiards (in particular equilateral triangles) on a spherical surface was studied in  to investigate the effect of the positive curvature on the classical motion. The structure of phase space for these curved triangular billiards turned out to be regular but very complex. As a consequence of the focusing mechanism on the sphere the phase space is entirely covered by chains of stable islands. Inside each of these islands the motion is elliptic and can be labelled by an infinite repeating code according to the sequence of reflections. This situation that we could call ’piecewise integrable’ is very different from the one corresponding to plane polygonal billiards which are either integrable (or pseudointegrable) or ergodic, according to the rationality of their inner angles (see, for example, ). The natural question we would like to answer in the present work is how this peculiar dynamics is reflected in the statistical properties of the quantum mechanical spectrum.

There is a lot of evidence, mainly numerical, that integrable dynamics leads to Poisson statistics while classically chaotic systems satisfy the statistics given by one of the Random Matrix Ensembles (see, e.g., ). Mixed systems show a combination of these extremes that depend on the relative size of chaotic and regular regions . On the other hand, a different kind of systems neither regular nor chaotic has been shown to conform to intermediate statistics: a simple example is the plane pseudointegrable billiard which has been studied together with other models with similar behaviour by Bogomolny et al .

The peculiar phase space structure of equilateral triangles on the sphere led us to believe that their spectral properties would not conform to any universal statistics. In this article we calculate the quantal energy spectrum and the eigenfunctions of these systems. For generic triangles the calculation will be performed numerically while the spectra of tiling triangles will be derived analytically just by using symmetry arguments, as was the case in the classical counterpart. It was shown in  that tiling triangles are very particular systems that not only are integrable but for which only periodic orbits are present. We will see how this peculiarity affects the quantal spectrum. In section 2 we present the model and the numerical procedure to calculate the spectra. Section 3 is devoted to the spectral properties of tiling triangles. In section 4 we derive the level spacing distributions of several generic triangles and discuss them in connection with the corresponding classical phase space. Finally conclusions are presented in section 5.

## 2 THE MODEL

Given $`T`$ an equilateral triangle with inner angle $`\omega `$ centered on the north pole of a sphere of radius $`R=1`$, we are concerned with the eigenvalue problem

$$\mathrm{}\psi _E+E\psi _E=0inTand\psi _E=0onT$$
(1)
where $`\mathrm{}`$ is the three-dimensional Laplacian in spherical coordinates $`\theta `$ and $`\varphi `$ given by:

$$\mathrm{}=\frac{1}{\mathrm{sin}\theta }\frac{}{\theta }(\mathrm{sin}\theta \frac{}{\theta })+\frac{1}{\mathrm{sin}^2\theta }\frac{^2}{\varphi ^2}.$$
(2)
and the eigenvalues $`E=\frac{2m}{\mathrm{}}`$, $``$ being the energies, are dimensionless.

The wave function $`\psi _E`$ can be expanded as

$$\psi _E=\underset{m}{}c_mF^m(E,\mathrm{cos}\theta )e^{im\varphi }$$
(3)
where $`m`$ is an integer and $`F^m`$ is a solution of the Legendre’s equation:

$$\frac{}{x}[(1x^2)F(x)]+[\frac{m^2}{(1x^2)}+E]F(x)=0$$
(4)
with $`x=\mathrm{cos}\theta `$. The functions $`F^m`$ which have to be regular inside the triangle are then proportional to the Legendre’s functions of the first kind $`P_\nu ^m(x)`$, with real $`\nu `$ such that $`E=\nu (\nu +1)`$ and real argument $`1<x<1`$. These $`P_\nu ^m`$ ( $`m0`$) go to $`0`$ at $`x=1`$, while they diverge for $`x=1`$.

It will be convenient to define :

$$F^m(E,x)=\sqrt{\frac{\mathrm{\Gamma }(\nu m+1)}{\mathrm{\Gamma }(\nu +m+1)}}P_\nu ^m(x)$$
(5)
These functions are real for $`m\nu +1`$ and satisfy the following recursion relation:

$`F^{m+2}`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{(\nu +m+2)(\nu m1)}}}2(m+1){\displaystyle \frac{x}{\sqrt{1x^2}}}F^{m+1}`$ (6)
$`\sqrt{{\displaystyle \frac{(\nu m)(\nu +m+1)}{(\nu +m+2)(\nu m1)}}}F^m`$
which is used in the numerical procedure to compute the $`F^m`$.

In order to obtain a spectrum free from the trivial degeneracies due to symmetry we have to desymmetrize the billiard. The spherical equilateral triangle is invariant under the point group $`C_{3v}`$ and can be cut along the $`3`$ reflection planes into $`6`$ triangular subdomains. Solving the eigenvalue problem in each of these irreducible domains with given boundary conditions corresponds to finding a fraction of the spectrum of a given symmetry class. In the following we will treat the Dirichlet problem in the fundamental triangle. This corresponds to solutions which are odd under reflection across all the symmetry planes. Subspectra corresponding to other symmetry classes can be obtained by imposing Neumann and mixed boundary conditions. In the case we are considering the expansion of Eq.(3) reduces to:

$$\psi _E=\underset{p}{}c_{3p}F^{3p}(E,\mathrm{cos}\theta )\mathrm{sin}3p\varphi $$
(7)
that vanishes automatically at $`\varphi =0`$ and $`\varphi =\frac{\pi }{3}`$.

In order to solve the problem we use an improved version of the ’point matching’ method introduced in . This consists in requiring that the Fourier coefficients of the wave function evaluated on the boundary

$$I_n=𝑑se^{2i\pi ns/}\psi _E$$
(8)
should vanish. $``$ is the billiard perimeter and $`ds`$ the line element. The existence of a non trivial solution leads to the condition:

$$det[J_{n,m}]=0,withn,m=1,N$$
(9)
where

$$J_{n,m}(E)=𝑑se^{2i\pi ns/}F^m(E,\mathrm{cos}\theta (s))e^{im\varphi (s)}$$
(10)
Eq. (9) determines the energy levels $`E`$.

The numerical calculation was performed for equilateral triangles with inner angle $`\frac{\pi }{2}\omega \frac{2\pi }{3}`$, that is fairly large triangles for which no difficulties are expected at the corners. Before presenting these numerical results we analyze in the next section the particular cases of tiling triangles that can be treated analytically.

## 3 TILING TRIANGLES

As already described in  the classical motion in tiling triangles can be studied by following a unique geodesic on the topological surface obtained by sewing together a finite number of replica of the original billiard. Since the geodesics are closed curves on this compact surface, all orbits are periodic. In other words, every trajectory in a tiling triangle is restricted to a one dimensional subspace. The system is over-integrable, the compactness of the sphere playing the role of an additional integral. We therefore expect a non generic behavior of the level distribution. The eigenvalues and eigenfunctions of tiling triangles can be evaluated by making use of symmetry arguments exclusively. Let us consider the equilateral triangle with $`\omega =\frac{\pi }{2}`$, that is, a triangle whose vertices coincide with those of a face of an octahedron. Following the desymmetrisation scheme of the previous section we cut the triangle into $`6`$ triangles with Dirichlet boundary conditions and calculate the ’desymmetrised spectrum’. Since the domain tesselates the sphere the solutions of Legendre’s equation (4) have to be continuous and one-valued at all points $`1\mathrm{cos}\theta 1`$. This restricts $`\nu `$ to integer values $`l`$ and $`|m|l`$. The seeked eigenfunctions will then be linear combinations of spherical harmonics $`Y_l^m`$ that can be determined by requiring that they should vanish on the boundaries. These boundary conditions result in a reduction of the spherical symmetry of the problem: the full spherical symmetry group including all proper and improper rotations in three dimensions will then be reduced to a point group. For this particular triangle it will be the symmetry group of the octahedron . That means that in order to determine the eigenvalues spectrum, that is, the allowed values of $`l`$ and their degeneracy $`\lambda _l`$ we should determine how each of the representations of the full spherical group $`D_l^\pm `$ may be decomposed into irreducible representations of $`O_h`$.

The number of times the $`\alpha ^{th}`$ irreducible representation of the subgroup $`O_h`$ is contained in the representation $`D_l^\pm `$ of the spherical group is given by :

$$c_\alpha =\frac{1}{g_{O_h}}\underset{R}{}\chi _\alpha (R)\chi _l(R)$$
(11)
where the sum extends over the elements of the group $`O_h`$ the order of which is $`g_{O_h}=24`$. The characters $`\chi _\alpha (R)`$ can be extracted from the corresponding character table while the characters for the proper rotations in the full rotation group are given by  :

$$\chi [D_l^\pm (R_\theta )]=\frac{\mathrm{sin}(l+1/2)\theta }{\mathrm{sin}\frac{\theta }{2}}$$
(12)
Here $`R_\theta `$ denotes a rotation through $`\theta `$ about some axis. The character of an improper rotation (that is, proper rotation times an inversion) is the same as the character for the proper rotation for the $`+`$ representations (that is, for even values of $`l`$) and its negative for the $``$ (odd values of $`l`$) representation.

By applying he $`24`$ symmetry operations of $`O_h`$ it can be easily seen that the eigenfunctions satisfying Dirichlet conditions in the reduced triangle transform as the one-dimensional irreducible representation $`\mathrm{\Gamma }_1^{}`$ of $`O_h`$ . Therefore the degeneracy $`\lambda _l`$ of the eigenstate of energy $`E_l=l(l+1)`$ will be the coefficient of the representation $`\mathrm{\Gamma }_1^{}`$ in the decomposition of $`D_l^\pm `$ according to Eq.( 11).

The first levels of the spectrum are shown in Table 1 with their corresponding degeneracies. The eigenfunctions have been obtained by means of diagonalizing the projector of the representation $`\mathrm{\Gamma }_1^{}`$ of $`O_h`$ : they are the eigenfunctions corresponding to eigenvalues equal to 1. Since the representation is one-dimensional the projector has been expressed in terms of the characters as:

$$P_{\mathrm{\Gamma }_1^{}}=\frac{1}{g_{O_h}}\underset{R}{}\chi _{\mathrm{\Gamma }_1^{}}(R)D_l^\pm (R)$$
(13)
Again the sum runs over all the elements of the group, and the D-matrices corresponding to improper rotations are minus the corresponding to proper rotations.

It is easy to understand that only odd values of $`l`$ appear in the spectrum of Table 1 since the representation $`\mathrm{\Gamma }_1^{}`$ appears only in the decomposition of $`D`$ representations with negative (natural) parity. Furthermore, it can be seen that for $`l13`$ all odd values of $`l`$ are present in the sequence. That means that except for the lowest lying levels the level spacing $`S(l)=E_{l+2}E_l`$ is proportional to $`l`$. In order to keep constant the classical mean level spacing $`\overline{S}`$ which according to Weyl’s law is $`\overline{S}=\frac{4\pi }{A}`$, $`A`$ being the area of the reduced triangle, the level degeneracy $`\lambda _l`$ should also increase linearly with $`l`$. The spectrum of the tiling spherical triangle is thus dominated by number-theoretic degeneracies. These ’accidental’ degeneracies have been studied by Itzykson et al . They appear in some simple integrable quantum systems like, for instance, the harmonic oscillator with rational frequency ratios ( where $`E_{m,n}=m+n+1`$ ) and the plane integrable polygonal billiards ( for the equilateral triangle $`E_{m,n}=m^2+n^2mn`$ ). The harmonic potential was extensively investigated by Berry and Tabor  . It is an example of overintegrable systems, in which the additional integrals are the conmensurability relations. The spacing among adjacent levels is constant while the degeneracy increases with the energy . This makes the mean level spacing go to $`0`$ in the classical limit and no spacing distribution can be defined. The case of the plane equilateral triangle presents more analogy with our case in what concerns the mean level spacing $`\overline{S}`$ which is well defined and given by Weyl’s formula. Also there the states are increasingly degenerate and separated by increasingly large gaps.

In order to study the statistics of the nearest neighbor spacing distribution $`P(S)`$ we now unfold the spectrum by defining the sequence $`e_l=\overline{N}(E_l)`$ where $`\overline{N}(E)`$ is the averaged integrated level density given by

$$\overline{N}(E)=\frac{A}{4\pi }E\frac{L}{4\pi }\sqrt{E}+const$$
(14)
$`L`$ is the triangle perimeter in the spherical metric. The level spacing distribution in the scaled spectrum $`P(s)`$, with $`s_l=e_{l+2}e_l`$ is plotted in Fig.1 (b) for a sequence of $`7500`$ levels. The distribution is bimodal with a strong tight peak at $`s=0`$ and a flat component, giving equal probability to all allowed values of $`s_l`$. Since both degeneracies and gaps diverge with the energy this distribution is not defined in the classical limit. Asymptotically the peak tends to a $`\delta `$\- function and the flat mode extends to infinite values of $`s`$ with a height going to $`0`$. For comparison in Fig.1 (a) we show the spacing distribution for a sequence of $`20000`$ energy levels corresponding to a desymmetrized plane equilateral triangle with Dirichlet boundary conditions . Although the distribution looks more poissonian, as in a generic regular system, the situation is analogous to the one in the curved triangle: as shown in  degeneracies and gaps diverge with energy. But, as pointed out by Berry in  for the similar case of the right plane triangle, this divergence is so slow that this nongeneric level structure governed by number-theoretic degeneracies only appears for very high-lying states. Summarizing, the spectrum of tiling spherical triangles is non generic and although a mean level spacing can be defined, no spacing distribution $`P(s)`$ exists in the classical case. This peculiarity is to be attributed to the number-theoretic structure of the spectrum . The spectra of the spherical tiling triangles with angles $`\omega =\frac{2\pi }{3}`$ and $`\omega =\frac{2\pi }{5}`$ can be calculated in the same way, by considering the symmetry group of the tetrahedron and of the icosahedron respectively. The sequence of eigenvalues corresponding to the desymmetrized triangle with $`\omega =\frac{2\pi }{3}`$ and Dirichlet boundary conditions, that transform as the one-dimensional irreducible representation $`\mathrm{\Gamma }_2`$ of $`T_d`$ is presented in Table 1. In this spectrum odd and even values of $`l`$ are allowed . Except for the lowest lying levels all values of $`l`$ are present, the level spacing $`S(l)=E_{l+1}E_l`$ is again proportional to $`l`$ and the spectrum has the general characteristics stated above.

## 4 GENERIC TRIANGLES

The numerical procedure introduced in section 2 was carried out for curved equilateral triangles with inner angle $`\omega `$ ranging from $`\frac{\pi }{2}`$ to $`\frac{2\pi }{3}`$. We checked the stability of the solutions by varying the number of points on the boundary used to evaluate the integral in Eq.(8) and the number of partial waves in the expansion Eq.(3). The completeness of each spectrum was tested by comparing the calculated cumulative level density $`N(E)`$ with the averaged one given by the Weyl formula (14). This test is particularly relevant for triangles close to the tiling ones, for which most of the levels are near-degenerate and therefore some of them could be missed in the numerical calculation.

In Fig. 2 we show the energy spectra up to $`E=2400`$ for a family of equilateral triangles labeled by $`\omega `$. These are plotted as curves $`E_\nu (\omega )`$. The spectra at the left and right boundaries correspond to the integrable triangles with $`\omega =\frac{\pi }{2}`$ and $`\omega =\frac{2\pi }{3}`$ respectively, which have been evaluated analytically as described in the previous section. The levels in both sequences are labelled by an integer $`l`$ and their multiplicity increase with energy.

All the remaining triangles with arbitrary rational or irrational angle $`\omega `$ are non-integrable and therefore present no degeneracies in their spectra. Accidental degeneracies are not expected neither: by varying $`\omega `$, we are moving in a one-parameter family and, according to Berry  and references therein, variations in one parameter are insufficient to produce degeneracies. However, in the region denoted as $`III`$ we observe a great number of quasi-crossings. Their origin will be explained below.

The level spacing distributions corresponding to the unfolded spectra of six triangles with $`\omega `$ ranging from $`\frac{\pi }{2}`$ to $`\frac{2\pi }{3}`$ are shown in Fig. 3. Between $`1300`$ and $`2500`$ levels were considered for each case, except for case (d) ($`\omega =0.589\pi `$) where the huge number of quasidegeneracies prevented us to go beyond a few hundreds levels. We observe that the distributions depend drastically on the parameter $`\omega `$ and do not exhibit any universal behaviour. In order to understand the general features of the distributions and their angle dependence we will have to keep in mind the results for tiling triangles derived in the previous section and also refer to the classical phase space plot which is shown for each case in Fig. 4. As seen in  for generic curved triangles the classical phase space is covered by chains of elliptic islands of regular motion characterized by an infinite repeating code. The multiplicity of these chains increases and their size decreases with the period of the code and phase space takes a fractal structure. As shown in Fig. 4 (a) and (f), for triangles close to the tiling ones phase space is almost dominated by one chain of $`3`$ islands (corresponding to the $`+++`$ code ) in the case $`\omega =\frac{\pi }{2}`$, and by two chains (corresponding to $`+++`$ and $`+`$ ) in the case of $`\omega =\frac{2\pi }{3}`$. The region outside these big domains is entirely covered by ’dust’, that is by chains of islands corresponding to very long codes and therefore with extremely high multiplicity and small area. These islands , although regular and of the same type as the dominant ones, cannot be resolved quantum mechanically in the region of the spectrum we are analyzing. As we go to triangles far apart from the tiling ones (see Fig. 4 (b),(c), (d),(e)) other chains of considerable size appear and we will see that some of them can be resolved in the range of energies considered. In order to see how this is reflected in the level spacing distributions we now go back to Fig. 3.

Histogram $`(a)`$ corresponds to an angle $`\omega =0.511\pi `$ ( region $`I`$ of Fig. 2) and can be seen as a perturbation of the tiling $`\omega =\frac{\pi }{2}`$. The $`O_h`$ symmetry is broken and the fact that degeneracies disappear and gaps between levels do not diverge with energy but tend to some finite value reflects in a modification of the distribution of Fig. 1 (a). The distribution Fig. 3 (a) is still bimodal but the delta peak at $`s=0`$ is shifted (i.e., there is some level repulsion) and the flat component does not extend to large values of $`s`$. Case (f) corresponding to $`\omega =0.649\pi `$ (also in region $`I`$ of Fig. 2 ) presents an analogous situation, with a breaking of the $`T_d`$ symmetry. But in this histogram the low lying peak is broader (in fact, there is no level repulsion) and the high s - component shows more structure suggesting a superposition of two uncorrelated spectra. In both cases the distributions, and more specifically the presence or not of level repulsion, can be better understood by analyzing the corresponding eigenfunctions. At a qualitative level, we expect that eigenstates spanning an integrable region will show rigid and regularly spaced spectra. Therefore, in case where a single integrable region dominates the phase space the spectrum will show repulsion, since all the eigenfunctions are correlated. On the contrary, if several classes of eigenfunctions coexist, the distribution can be thought as a superposition of uncorrelated spectra and there will be no level repulsion. In case (a) the eigenfunctions are of two types: one class lives in the $`+++`$ elliptic islands which mostly fill the phase space, the other class extends over the region left outside the islands, which is a very small fraction of the space. (see Fig. 4 (a) ). Thus, the eigenfunctions of the first class, which are the majority, will be correlated leading the level repulsion observed in Fig. 3 (a). In case (f) the eigenfuctions are of three types. Apart from the extended class which lives in the region left outside the islands, there are two classes of localized functions: one living in the $`+++`$ elliptic islands, the other in the $`+`$ islands (see Fig. 4 (f) ). Examples of the Husimi representation of functions belonging to these three classes are shown in Fig. 5. The two classes of localized eigenfunctions living in domains which are comparable in size will be weakly correlated : this explains the much weaker level repulsion in distribution of Fig. 3 (f).

The distributions corresponding to $`\omega =0.525\pi `$ and $`\omega =0.633\pi `$ are shown in Fig. 3 (b) and (d) respectively. Both correspond to values of the parameter in region $`II`$ of Fig. 2 . We are now fairly far away from the tiling triangle: the splitting among levels is big, and the two modes have shifted yielding a single peaked distribution. As we see in the corresponding classical phase space plot Fig. 4 (b) for $`\omega =0.525\pi `$ the islands corresponding to codes different from the dominating $`+++`$ are still small and cannot be resolved in the range of energies we are considering. Still, there are two classes of eigenstates coexisting: the ones localized in the $`+++`$ domain and the extended ones which occupy now a considerable fraction of phase space. Therefore two classes of weakly correlated eigenfunctions contribute to the spectrum. For $`\omega =0.633\pi `$ three classes of functions are present, since the $`+`$ islands are large, as shown in Fig. 4 (d), and can be resolved.

Summarizing, we can say that the spectra of generic curved triangles are non universal and follow a parameter- dependent intermediate statistics. The analysis of the corresponding classical phase space plots suggest that the characteristics of the distributions depend on the number of domains that can be quantum mechanically resolved in a given energy range. The conjecture is the following. If only one big domain is resolved the eigenfunctions in this domain are of the same type and, therefore, correlated. If more than one domain can be resolved, the coupling among states living in different domains will be weak, and the spectrum will be a superposition of uncorrelated spectra, showing no repulsion.

Another interesting observation is the existence of a high peak at a small value of $`s`$ in the level spacing distributions corresponding to triangles in region $`III`$ of Fig. 2. Two examples are shown in In Fig.3 (c) and (d) for angles $`\omega =0.571\pi `$ and $`\omega =0.589\pi `$ respectively. This peak is due to the presence of numerous quasi-crossings in this region, the lowest of them are seen already in Fig.2. To understand the origin of these quasi-degeneracies we have to refer once again to the classical phase space plots corresponding to these triangles. As pointed out in  generic non tiling triangles might have domains in which all trajectories are periodic (coexisting with the families of elliptic islands). For example, in all rational triangles with $`\omega >\frac{\pi }{2}`$ the islands corresponding to the $`+`$ code are constituted of periodic orbits. On the other hand, in triangles with inner angle satisfying

$$\mathrm{cos}\frac{\omega }{2}=\sqrt{\frac{3}{4}\mathrm{cos}^2\frac{k\pi }{n}}$$
(15)
all orbits of type $`+++`$ are periodic with period $`n`$.

In region $`III`$ we find several triangles with rational inner angle or inner angle satisfying Eq. (15), corresponding to periodic orbits of relatively short period. For example, in case (d) $`\omega `$ satisfies Eq. (15) with period $`n=7`$. Although there is no global symmetry group, as in the case of tiling triangles, the presence of these periodic domains which cover a big fraction of phase space (in particular the corresponding to the $`+++`$ code) lead to quasi degeneracies among the quantum states localized in these domains. This is the origin of the huge peak at $`s=0`$ in Fig. 3 (d). Distribution (d) corresponds to an angle close to $`\frac{4\pi }{7}`$ : now the periodic orbits are of type $`+`$ with period $`n=7`$. Since the area occupied by the $`+`$ islands is smaller, the effect is less pronounced: the peak is lower and located at $`s`$ small but different from $`0`$.

## 5 CONCLUSIONS

In this work we calculated the quantal spectra of generic and tiling equilateral triangles on a spherical surface. The spectra of the tiling triangles could be obtained analytically by using symmetry arguments while for the generic ones we performed a numerical calculation based on an improved version of the ’point matching’ method which allowed us to go up to around $`2000`$ levels.

We found that the spectra of tiling triangles are dominated by arithmetical degeneracies. This is a property shared by other simple quantum systems, in particular by plane square and triangular billiards. In the case of curved triangles the distribution of the degeneracies is easily determined: since the level gap increases linearly with $`l`$, that is, is proportional to $`\sqrt{E}`$ and the mean level spacing in the classical limit is constant, the degeneracies also increase as $`\sqrt{E}`$. But the question asked in  about relating these degeneracies with some ’hidden’ symmetry remains open also in this case.

As expected, the spectra of generic triangles of variable inner angle $`\omega `$ show no universal behavior in this energy range. We first remarked that for triangles close to the tiling ones the level spacing distribution is bimodal : this can be seen as a remnant of the level clustering in the tiling systems. More in general, for an interpretation of the spectral properties we had to refer to the corresponding classical phase space plots. These might be governed by one or more than one island of considerable size, where considerable means that it can be quantum mechanically resolved in the energy range under consideration. The conclusions seem to be the following. If there is only one resolved domain all the eigenfunctions living in this domain will have support on tori and the corresponding energy levels will be regularly spaced, giving rise to a rigid spectrum. This is the case in triangles with inner angle $`\omega `$ close to $`\frac{\pi }{2}`$ for which the $`+++`$ islands covers a great area of the phase phase ((a) and (b) in Fig.4). If several domains of comparable size coexist different eigenfunctions which may correspond to levels close in energy are localized far apart from each other and are weakly correlated. No level repulsion will be present in the spectrum.

In order to check these conjectures and to get a quantitative understanding of the spectra we have to extend the calculation to higher energies. The higher we go, more and more domains in phase space will be resolved and an increasing number of weakly correlated families of eigenfunctions will be present, localized in these domains. It will then be possible to subdivide the spectrum in subspectra corresponding to each type of eigenfunctions (according to their localization in phase space) and study the characteristics of each subspectrum. This subdivision is not possible in the present work since we cannot reach values of $`k`$ for which the small islands will be explored with significant statistics. To remedy this situation and check whether in the high $`k`$ regime the superposition of many uncorrelated spectra recovers asymptotically the Poisson behaviour, we are at present adapting the scaling method  to the geometry of the sphere.

### Acknowledgements

We are grateful to N.N. Scoccola for valuable discussions. This work was partially supported by PICT97 03-00050-01015 and CONICET PIP 0420/98.

Fig.1 – (a) shows the level spacing distribution corresponding to the first $`20000`$ levels of a desymmetrized plane equilateral triangle. In (b) the level spacing distribution for the first $`7500`$ levels of a desymmetrized curved tiling triangle with $`\omega =\frac{\pi }{2}`$.

Fig.2 Energy spectra for a family of equilateral triangles with $`\frac{\pi }{2}\omega \frac{2\pi }{3}`$ . Full points indicate the spectra of the tiling triangles. Region I corresponds to the neighborhood of tiling triangles, while II and III denote intermediate regions (see text).

Fig.3 Level spacing distributions for 6 generic triangles .

Fig.4 Phase space plots for 6 generic triangles.

Fig.8 Semiclassical representation of three types of eigenfunctions in the triangle with $`\omega =0.649\pi `$ .