text:
# The Promiscuous Nature of Stars in Clusters

## 1 Introduction

The rich environment of a star cluster provides an ideal laboratory for the study of self-gravitating systems. Star clusters within our Galaxy range in size from loose associations and open clusters, which contain up to tens of thousands of stars, to the large globular clusters which host a million stars or more. Globular clusters are some of the oldest objects known and dynamically they lie in a very interesting regime. Compared to the solar neighbourhood the central density of stars is high enough (a factor of 10 million or more greater in some globular clusters) that a significant fraction of the cluster stars are likely to experience at least one close encounter with another star during their lifetime. At the opposite end of the scale, the size of a globular cluster is small enough, in comparison to a galaxy, that its dynamical relaxation timescale is less than its age. These combined considerations make a globular cluster an exciting place for a star to reside. As such we expect the stellar populations of a star cluster, the various classes of stars and binaries, to exhibit a dynamical signature. In particular, many of the stars contained in the cluster will deviate strongly from the evolutionary paths predicted by standard stellar and binary evolution theory. This is indeed what we observe: either when pointing the Hubble Space Telescope (HST) at the centres of globular clusters or when generating star-by-star models of star cluster evolution.

## 2 Observations of Exotic Objects

Observational evidence for the modification of stellar populations in star clusters owing to dynamical interactions between the cluster stars is found in many instances (see Bailyn 1995 for a review). A common indicator of this effect is a variation in the radial distribution of a particular stellar sub-population within the cluster. If, for example, a certain type of exotic star is found primarily in the core of a cluster, where the number density of stars is highest and hence interactions between stars are most likely, this is taken as strong evidence that these stars have a dynamical origin.

Blue stragglers (BSs) are stars burning hydrogen in their cores that, according to their stellar mass, should already have evolved off the main-sequence (MS) and past the red giant stage to become white dwarfs (WDs): stellar corpses. Somehow these stars have been rejuvenated. Popular theories of the cause of this phenomenon are mass transfer from a companion star in a close binary system or the collision of two MS stars. HST observations of the dynamically old and metal-poor Galactic globular cluster M30 have revealed a number of BSs (Guhathakurta et al. 1998). When plotted in a cumulative radial distribution these BSs appear concentrated towards the cluster core relative to standard MS stars. Conversely, the population of bright red giants in M30 is depleted in the core. This suggests that collisions are destroying giants, by ripping off their envelopes and reducing them to WDs, and creating BSs. Ferraro et al. (1999) used HST to identify 305 BSs in the high density globular cluster M80. These are also centrally concentrated with respect to the giants in the cluster and the suggestion is that M80 is being observed at a critical phase in its evolution when stellar interactions are delaying the collapse of the core (Ferraro et al. 1999) - leading to the production of collisional BSs.

Blue Stragglers are not the only stars exhibiting enhanced populations in the cores of globular clusters. X-ray binaries are at least $`\mathrm{1\hspace{0.17em}000}`$ times more abundant in clusters than in the Galactic field, and are being found in ever-increasing numbers with the Chandra X-ray observatory (Grindlay et al. 2001). Millisecond pulsars, including some in binaries or with planets (D’Amico et al. 2001) are prevalent in cluster cores. Cataclysmic binaries in clusters have been predicted and are also being found (Shara et al. 1996; Grindlay et al. 2001). A few sdB stars have been located and characterized (Moehler et al. 1997).

Even though the stellar neighbourhood within an open cluster is less dense than that of a globular cluster it is still capable of producing exotic objects. The number of BSs found in the open cluster M67 is much greater than we would expect if they are simply produced via mass transfer in binaries exhibiting the same distribution of orbital characteristics as binary stars found in the field. Furthermore, these BSs have a variety of living arrangements (Leonard 1996, and references within). Some are single while others are found with a companion. In some cases the BS and the companion star interact in an intimate and regular manner (short-period circular orbit), in other cases the relationship is distant (long-period orbit) and may even appear eccentric. One of the BSs is so massive, a super-BS, that it must represent the merger of three stars, and another is observed in an active triple-system, a menage-a-trois (van den Berg et al. 2001). All of this suggests that the BSs have varied formation histories and that no one mechanism is responsible for their production (Leonard 1996). In particular, the presence of a super-BS and of BSs in eccentric binaries is not predicted by standard binary evolution and points to dynamical interactions tampering with the destinies of stars.

## 3 Simulation Methods

On a basic dynamical level a star cluster is composed of $`N`$ bodies interacting with each other due to the gravitational force of every other body in the system. The ideal method for following the evolution of such a system is to integrate directly the $`N`$ individual equations of motion, the $`N`$-body approach. However, the cost of integrating the cluster for a relaxation-time scales as $`N^3`$$`N`$ per force calculation, $`N^2`$ to integrate $`N`$ bodies for one crossing (dynamical) time, and there are of order $`N`$ crossing times per relaxation time – so the method is computationally expensive. As a result most $`N`$-body simulations performed until recently, although state-of-the-art at the time, have involved a varying number of simplified and unrealistic conditions, such as including only single stars, using only equal-mass stars, neglecting stellar evolution or assuming no external tidal field (McMillan, Hut & Makino 1991; Heggie & Aarseth 1992). Performance has been enhanced by the development of improved computational algorithms. For example, the use of individual time-steps (Aarseth 1963) enables each star to evolve on its own natural dynamical timescale rather rather than forcing a small time-step on the entire system when only a fraction of the stars may need frequent attention. Quantizing these time-steps into a series of hierarchical levels increases the efficiency of the integration (McMillan 1986). Also, regularization techniques have been developed to remove the singularities involved with compact binaries and few-body configurations, as well as to improve the accuracy of the treatment (Mikkola & Aarseth 1998).

The very nature of the large-$`N`$ effects which make the $`N`$-body method so computationally expensive can be used to justify the use of statistical algorithms such as the orbit-averaged Fokker-Planck equation (Hénon 1961; Giersz 2001; Joshi, Nave, & Rasio 2001). In theory the $`N`$-body and Fokker-Planck methods, when applied to the same initial conditions, should give similar results in the limit of large $`N`$ (Takahashi & Portegies Zwart 2000). Ultimately the $`N`$-body approach remains the method of choice for creating dynamical models of star clusters as a minimum number of simplifying assumptions are required and it is relatively easy to implement additional realistic features (see Meylan & Heggie 1997 for a full discussion). It is also possible to follow the journey made by a particular star as the cluster evolves.

A major impact on the suitability of $`N`$-body methods has been recent advances in computing power, in particular the development of special-purpose hardware. The GRAPE-4 system became commercially available in 1995 and was designed specifically for high-accuracy simulations of dense stellar systems (Makino & Taiji 1998). It basically acts as a Newtonian force accelerator: the host workstation supplies the GRAPE (GRAvity PipE) with the positions, velocities and masses of a list of particles and the GRAPE returns the force and its time derivative for each particle. A 96-chip GRAPE-4 board operates at a peak speed of $`50`$Gflops (50 billion floating-point operations per second) and has enabled production of realistic open cluster models in a reasonable time. To be precise, a simulation of $`\mathrm{10\hspace{0.17em}000}`$ stars with a moderate binary fraction of 10% could be completed within a week. Dealing with binary systems complicates matters because as far as the GRAPE is concerned it sees binaries solely as the centre-of-mass particle and the integration of the orbits of the stars within each binary system must be dealt with on the host machine. Therefore the inclusion of a large proportion of binaries in a simulation is a potential bottle-neck to the computational performance. These concerns are non-trivial owing to the fact that upwards of 50% of the stars in any one globular cluster may be binaries and that the majority of these will be primordial (Hut et al. 1992). On the upside it has been shown by Makino & Hut (1990) that the computational cost of primordial binaries becomes (relatively) less of a problem when $`N`$ is increased.

The next generation machine, the GRAPE-6 (Makino 2001), became available in early 2001 and has provided a further leap forward for $`N`$-body methods. A single GRAPE-6 chip has an operating speed of $`30`$Gflops so that a 32-chip board represents $`1`$Tflops of computing power. The peak speed of GRAPE-6 is up to 100 times faster than GRAPE-4. This exciting development means that models of $`\mathrm{50\hspace{0.17em}000}`$ stars can be simulated from formation to death, i.e. $`10`$Gyr, in a week of wall-clock time. As such, it is now possible to produce realistic $`N`$-body models of small globular clusters for the first time. The largest Galactic globular clusters are still out of reach for a single GRAPE-6 board. However, at the University of Tokyo, where the GRAPE project is based, work is underway to link many GRAPE-6 boards in parallel and simulate systems containing upwards of a hundred-thousand stars111We refer the interested reader to astrogrape.org for further information on the GRAPE project.. To tackle the million body problem we will have to wait for the GRAPE-8 which is expected to run at Petaflops speed.

The Aarseth NBODY4 code (Aarseth 1999) has been developed to investigate the evolution of star clusters utilising the GRAPE hardware. Importantly NBODY4 incorporates algorithms for the detailed treatment of stellar evolution and the full range of possible interactions within binary stars (Hurley et al. 2001, and references within). In the dense environment of a star cluster it is possible for the orbital parameters of a binary to be significantly perturbed owing to close encounters with nearby stars. As a result the orbit may even become chaotic (Mardling & Aarseth 2001). Collisions may occur and binaries can form as a result of tidal capture (Fabian, Pringle & Rees 1975). Gravitational encounters between single stars and binaries, or binary-binary encounters, can lead to an exchange interaction, where one star finds it energetically preferable to leave its current partner and become bound to another (Heggie 1975; Heggie & Hut 2002). Other possible outcomes of these scattering incidents include the destruction of the one (or two) binary system(s) to leave only single stars or formation of a triple system. Rather than being instantaneous events, exchange interactions often involve a period of indecisiveness where stars swap between partners in quick succession before choosing a companion to remain with (see McMillan, Hut & Makino 1991 for an early example of such behaviour). All these processes are accounted for in NBODY4. This allows the generation and interaction of the full range of stellar populations possible. A similar effort is represented by the kira code (Portegies Zwart et al. 2001b) which also takes advantage of the GRAPE hardware.

## 4 Promiscuous Stars

### 4.1 Cluster Models

The use of $`N`$-body codes on the GRAPE hardware has already given rise to a number of important results regarding star cluster evolution. Aspects studied include the effect of dynamical evolution on the mass-functions of globular clusters (Vesperini & Heggie 1997), the effect of the Galactic tidal field on the evolution of star clusters (Giersz & Heggie 1997), and the evolution of cluster morphology (Boily, Clarke & Murray 1999). The dependence of cluster lifetimes on their size and the tidal field in which they reside was recently investigated by Baumgardt (2001) who found that the lifetimes do not scale with the relaxation timescale, as had previously been assumed. Portegies Zwart et al (2001a) simulated the evolution and disruption of young compact star clusters near the Galactic centre and discussed the likelihood of their being observed.

Hurley et al. (2001) demonstrated the effectiveness of the cluster environment in modifying the nature of the stars it contains by modelling the BS population of M67. Not only was the number of BSs produced twice that expected from standard binary evolution but formation paths for all the various types observed in M67 were created. The models of young Galactic open clusters performed by Portegies Zwart et al. (2001b) also exhibited strong dynamical activity. These studies have highlighted the importance of combining dynamics with stellar evolution.

Within each of these open cluster size simulations there exist numerous interactions between stars that are often glossed over when discussing the simulation result. However, each of these interactions is interesting in its own right and we now take the opportunity to illustrate some of these fascinating tales in detail.

### 4.2 Lurid Examples

The first documentation of the behaviour of a particular star in a $`N`$-body simulation involving stellar evolution, and the motivation for this paper, was presented by Hurley et al. (2001). This star began life as a $`1.33M_{}`$ single star. Four exchange interactions and one tidal capture event later, involving a total of four collisions within eccentric binaries, its mass (after $`\mathrm{4\hspace{0.17em}100}`$Myr of cluster evolution) had grown to $`7.7M_{}`$, an amazing factor of six greater than the MS turn-off mass at that time. At various intervening points in the evolution this star would have been observed as a single super-BS, a super-BS in a binary, a BS in an eccentric binary, or a BS in a triple system. What we want to emphasise here is that such promiscuous behaviour is not at all rare.

Consider a similar case involving BS formation that arose in a recent simulation on GRAPE-6 aimed at modelling the evolution of planetary systems in star clusters (Hurley & Shara 2002). At the beginning of the simulation the star of interest had a mass of $`0.99M_{}`$ and was in a primordial binary with a period of $`P=\mathrm{6\hspace{0.17em}513}`$yr and eccentricity $`e=0.73`$. The companion star had a mass of $`0.54M_{}`$ and the binary was situated just beyond the cluster half-mass radius. The metallicity of the stars in this simulation was $`Z=0.004`$. After $`40`$Myr the orbital parameters of this wide binary had been reduced to $`P=63`$yr and $`e=0.13`$ owing to perturbations resulting from weak gravitational encounters with nearby stars. At $`T490`$Myr the binary had drifted inside the cluster half-mass radius and was involved in an exchange interaction with a $`0.29M_{}`$ star. Subsequently the low-mass interloper was ejected from the 3-body system leaving the original binary intact. Following this tumultuous but short-lived relationship, the binary partners remained quietly monogamous until becoming involved in a 5-body interaction at $`T=\mathrm{3\hspace{0.17em}463}`$Myr, by which time the binary had sunk further towards the cluster centre, through mass segregation, lying slightly outside of the core. The other participants in this encounter were a $`1.12M_{}`$ single star and a binary with component masses of $`1.17`$ and $`0.32M_{}`$ and $`P=58`$yr. Initially the $`0.54M_{}`$ star was ejected, destroying the primordial binary, and a quasi-stable 4-body system remained. Shortly afterwards the $`0.32M_{}`$ star was also ejected and the $`0.99M_{}`$ and $`1.12M_{}`$ stars formed a binary with $`P=51`$yr and $`e=0.95`$ in a triple system with the $`1.17M_{}`$ star. The presence of the third body gradually drove the eccentricity of the inner binary up to 0.99 at which point ($`T=\mathrm{3\hspace{0.17em}515}`$Myr) the two MS stars collided to form a $`2.11M_{}`$ BS which did not remain bound to the third star. Figure 1 documents this interaction, and all subsequent interactions involving the initially $`0.99M_{}`$ star, in terms of the masses of the stars involved. We note that the $`0.54M_{}`$ star actually escaped from the cluster, with a velocity exceeding the stellar velocity dispersion by a factor of three, approximately $`10`$Myr after it was removed from its primordial companion \- but as a direct result of the energy exchanged in that interaction.

The BS (now significantly more massive than the average cluster star) quickly sank inside the core of the cluster where, at $`T\mathrm{4\hspace{0.17em}000}`$Myr, it was involved in an exchange interaction with a $`2.2M_{}`$ single star and a primordial binary. The single star, itself a BS formed via the coalescence of two MS stars in a semi-detached binary, formed a bound pair with the original BS. This BS-BS binary had a period of $`487`$yr and an eccentricity of 0.72. Its orbit was perturbed by interactions with a $`1.2M_{}`$ star at $`T=\mathrm{4\hspace{0.17em}049}`$Myr and with a $`0.9M_{}`$ star at $`T=\mathrm{4\hspace{0.17em}092}`$Myr which increased the eccentricity to 0.98. The orbit then became chaotic and at $`T=\mathrm{4\hspace{0.17em}162}`$Myr the BSs collided to form a $`4.31M_{}`$ super-BS. At $`T=\mathrm{4\hspace{0.17em}180}`$Myr the super-BS exchanged itself into a binary comprising a $`1.11`$ and a $`1.52M_{}`$ star and left bound to the more massive of the two, also a BS. Then, in a final showdown, this binary formed a 4-body system with another binary and at $`T=\mathrm{4\hspace{0.17em}210}`$Myr the $`4.31M_{}`$ super-BS and the $`1.52M_{}`$ BS collided to form a $`5.78M_{}`$ super-BS. Being as massive as it was, almost five times greater than the MS turn-off mass at the time, this promiscuous star quickly evolved off the MS and shortly afterwards lost its envelope on the asymptotic giant branch to become an oxygen-neon-WD.

This lurid example, and the one presented by Hurley et al. (2001), highlight an important point. Though the possibility of any particular star becoming a BS as a result of a dynamical encounter is quite random, once a star does become a BS and hence one of the most massive stars in the cluster, subsequent interactions are almost inevitable. They also show that if a binary emerges from the indiscriminate and short-lived relationship that is an exchange interaction it will more often than not comprise the two most massive stars involved in the interaction, i.e. in this society it is desirable to be heavy.

In a typical GRAPE-6 simulation with $`N=\mathrm{20\hspace{0.17em}000}`$, initially comprised of $`\mathrm{18\hspace{0.17em}000}`$ single stars and $`\mathrm{2\hspace{0.17em}000}`$ binaries, where the evolution was followed for $`5`$Gyr, the number of exchange interactions observed was $`500`$. These involved 730 different stars with some stars taking part in multiple interactions. The number of stars that swapped partner once was 494, twice was 105, three times was 48, four times was 27, and 14 stars swapped partner on five occasions. Amongst these were a number of re-marriages where a “star changed its mind” and returned to its original partner. The component stars of one particularly flirtatious binary were actually involved in 22 exchange interactions, including 10 re-marriages.

The total number of merger events observed during the entire simulation was 104. Of these, mass transfer in a primordial binary accounted for 46 cases while 13 mergers came from mass transfer in a binary formed via an exchange interaction. The remainder were the result of collisions in eccentric binaries: 21 in primordial systems where the orbit was strongly perturbed by nearby stars.

BS formation as a result of a hyperbolic collision is rare in simulations of open clusters but does occur. One example involved a $`0.32M_{}`$ star that began life in the core of the cluster and slowly drifted out to the half-mass radius, owing to mass-segregation, where at $`\mathrm{1\hspace{0.17em}300}`$Myr it collided with a $`0.93M_{}`$ star. The relative velocity of the two stars at infinity was $`2.8\mathrm{km}\mathrm{s}^1`$ and the collison product was assumed to be a fully-mixed $`1.25M_{}`$ MS star. When the cluster was $`\mathrm{3\hspace{0.17em}850}`$Myr old this star first appeared as a BS and by $`\mathrm{4\hspace{0.17em}500}`$Myr, when the simulation ended, it had sunk inside the cluster core. The incidence of direct collisions will be greater in the higher density conditions of a globular cluster simulation. However, these are abrupt encounters and much less interesting than the sociable nature of exchange interactions and the binary systems they produce. The presence of binaries also acts to magnify the chance of collisions because in the case of a binary it is the semi-major axis that sets the relevant cross-section for collision, rather than the stellar radius which is used in the case of single stars. Here we have to be careful with the terminology used. If a binary is involved in a hyperbolic interaction any collision that results will occur in a two-step process and will not be direct: firstly a (resonant) capture may produce a hierarchical system and subsequently a collision can occur within this system.

Interesting cases of stellar interactions are not limited to BS formation alone. Consider the case of the primordial binary comprised of MS stars of 0.6 and $`1.5M_{}`$ with an eccentricity of 0.42 and an orbital period of $`270`$yr that was part of the same GRAPE-6 simulation mentioned above. Residing in the core of the cluster this binary suffered a series of weak perturbations to its orbit which drove the eccentricity up to 0.95. After the $`1.5M_{}`$ star became a sub-giant it was involved in two exchange interactions, finally settling into a $`33`$yr orbit with eccentricity of 0.9 about a $`1.4M_{}`$ MS star after $`\mathrm{2\hspace{0.17em}130}`$Myr of evolution. The sub-giant then evolved onto the giant branch (GB) and by $`T=\mathrm{2\hspace{0.17em}160}`$Myr tidal forces within the binary had circularized the orbit resulting in a separation of $`250R_{\mathrm{}}`$. While on the GB the $`1.5M_{}`$ star filled its Roche-lobe and a phase of common-envelope evolution began. This stripped the envelope of the giant and left a helium-WD and a MS star separated by $`52R_{\mathrm{}}`$. The MS star subsequently evolved to the GB, filled its Roche-lobe, and another common-envelope event ensued, resulting in a pair of 0.4 and $`0.3M_{\mathrm{}}`$ helium-WDs with an orbital period of $`0.7`$d. Owing to gravitational radiation this system would easily merge within $`10^{10}`$yr to possibly form a blue sub-dwarf star (Iben 1990), i.e. a helium-burning object with a thin hydrogen envelope.

The final case that we choose to highlight involves the formation of a Thorne-Żytkow object (TZO, Thorne & Żytkow 1977). Although TZOs have not been directly observed, they are thought to result from stellar mergers involving either a neutron star or a black hole, where the merger product is unstable and rapidly ejects all the material involved other than the neutron star or black hole, which remains. The system of interest began as a primordial binary with component masses $`10.8`$ and $`5.3M_{}`$, an eccentricity of 0.87 and a period of $`\mathrm{1\hspace{0.17em}870}`$d. After $`20`$Myr of evolution the orbital period had been reduced to $`\mathrm{1\hspace{0.17em}350}`$d owing to perturbations that had hardened this already hard222For a binary to be termed hard the magnitude of its binding energy must be greater than the mean kinetic energy of the cluster stars, otherwise it is soft. In stellar dynamics it is found that hard binaries will become harder, and soft binaries will be broken-up, as a result of encounters with other stars. binary. A few Myr later the more massive of the two stars evolved onto the sub-giant branch and, as the region of convection within its envelope grew, tidal forces began to circularize the orbit. At $`T=23`$Myr the primary star filled its Roche-lobe and began transferring mass to the companion. When the primary evolved onto the giant branch the rate of mass transfer accelerated so that the giant quickly overfilled the Roche-lobes of both stars to leave the $`2.42M_{}`$ helium core of the giant and the $`5.31M_{}`$ MS star contained within a common-envelope. Orbital friction then caused these two objects to spiral towards each other and the energy released was enough to drive off the envelope before they coalesced. This left a naked helium star and a MS star in a circular orbit with a $`4.3`$d period. The helium star then evolved onto the helium GB and filled its Roche-lobe. A period of stable mass transfer ensued until the envelope of the helium star had been removed and it became a $`1.42M_{}`$ WD of oxygen-neon composition. The orbital period at that stage was $`10`$d and the companion star was now a $`6.07M_{}`$ BS. At $`T=74`$Myr the BS evolved off the MS onto the sub-giant branch and filled its Roche-lobe which lead to another common-envelope event and the formation of a naked helium star in an orbit of $`P=0.1`$d with the oxygen-neon WD. When this helium star evolved onto the helium GB it also filled its Roche-lobe and began transferring mass to the WD. It was assumed that the helium-rich material would be accepted by the WD and steadily burnt on its surface. This quickly caused the mass of the WD to reach the critical Chandrasekhar mass, $`1.44M_{}`$, at which point it underwent an accretion-induced collapse to form a neutron star. The neutron star remained bound to the helium giant with an orbital period of $`0.04`$d. At $`T=94`$Myr the helium giant again filled its Roche-lobe which lead to a third phase of common-envelope evolution and the formation of a TZO.

In this particular instance dynamical interactions were not crucial to the outcome of the evolution: without perturbations to its orbit this binary would still have formed a TZO. However, we do see examples of similar systems being disrupted by encounters with other stars in the cluster so that the ultimate fate of the system is altered. Such occurences include the orbit being so strongly perturbed that the stars merge while both are on the MS. Alternatively a perturbation takes place just prior to the onset of a common-envelope phase, in a way that leads to coalescence of the two stars, rather than the formation of a close binary.

The examples given here represent only a small fraction of the wealth of information regarding stellar interactions and populations that is generated by realistic $`N`$-body simulations of star clusters. Much of this information is quickly discarded if it is not of direct relevance to the initial purpose of that particular simulation. With GRAPE-6 increasing the $`N`$ that can be simulated in a reasonable timeframe there exists a distinct possibility that those working in the $`N`$-body field will rapidly be overwhelmed with data. Therefore, the creation of an archive for the results of $`N`$-body simulations that is accessable to all members of the scientific community is crucial (see Teuben et al. 2002 for a related discussion). This database will complement the HST stellar populations archive (Zurek: HST Proposal AR-9225) currently under construction at the American Museum of Natural History.

It should also be noted that these interactions are not simply interesting: they have wide reaching implications. Take, for example, the case of type Ia supernovae which are important standard candles in cosmology and have recently been used to demonstrate that the expansion of the Universe is accelerating (Perlmutter et al. 1999). Short-period double-degenerate binaries in which the component WDs have a combined mass in excess of the critical Chandrasekhar mass and will merge within a Hubble time are thought to be possible progenitors of type Ia supernovae. Close inspection of the results of $`N`$-body simulations performed on GRAPE-6 have shown that the production rate of these systems is remarkably enhanced - over an order of magnitude in open clusters and likely much more in globulars - relative to the field (Shara & Hurley 2002).

## 5 Summary

The introduction of the GRAPE-6 hardware, coupled with $`N`$-body codes that have the capability to produce realistic cluster models, places the astrophysicist in the position of knowledgeable voyeur. An exciting aspect of this is the capability to investigate and understand the range of stellar populations that are observed, and how these are affected by dynamical interactions between cluster stars. As we move to full utilization of the GRAPE-6, and the simulation of globular clusters, the fascinating tales of promiscuous stars in clusters will become very graphic indeed!

We are extremely grateful to Piet Hut and Jun Makino for their efforts in instigating and maintaining the GRAPE program. Immense gratitude is extended to Sverre Aarseth for his tireless work in producing state-of-the-art $`N`$-body codes. We also thank Peter Eggleton, Rosemary Mardling, Steve McMillan, Seppo Mikkola, Onno Pols and Christopher Tout for important contributions made to the development of realistic $`N`$-body codes. Thanks to Simon Portegies Zwart and David Zurek for many stimulating discussions relating to the interactions of stars in clusters, and to the referee (Piet Hut) for a number of highly constructive suggestions. Doug Ellis’ generous support has enabled AMNH to purchase new GRAPE-6 boards and we are most grateful to him.


format: text

text:
# 1 Introduction

## 1 Introduction

The objective of this simple consideration is to attract slightly more careful attention to the interplay between electromagnetic and gravitational interactions, than usually is required for cross sections calculation. Total cross section for the typical elementary electrogravitational process like that considered in this note, $`\mathrm{}^++\mathrm{}^{}\gamma +`$ graviton $`𝒢,`$ is constant proportional to $`\alpha G_N`$ with $`\alpha 1/128`$ and $`G_N2.610^{20}yb,`$ while typical elementary cross section for electroweak process (with gauge boson instead of $`𝒢`$) is of the order $`\alpha ^2/s,`$ where $`s`$ is the energy squared of collision. It follows, that both cross sections should equate at energies of the order $`\sqrt{s}10^{18}`$ GeV. This huge gap is the reflection of the hierarchy problem, one of the ways how to circumvent this problem was proposed recently  and outlined in a cursory manner below.

First attempt of unified description of Electromagnetism and Gravity was undertaken by Theodor Kaluza . He achieved his goal by adding an extra, fifth dimension to the visible four-dimenshional space-time and interpreted the $`\mu 5`$ $`(\mu =1,2,3,4)`$ components as the electromagnetic vector potential. Afterwards Oscar Klein  suggested that this fifth dimension has a peculiar, periodic topology, and, therefore, he ”compactified” it to a circles of small radii, attached to each point of our visible world.

The fruitful Kaluza’s idea revivaled in 1980’s  and, more recently, with the advent of theories, known under the generic name ”M theory” , which inevitabely lead to consideration of ten or eleven (space) dimensions. At last, at the end of preceeding century, interest to the extra dimentions greatly arose, because, on the one hand these theories may provide a natural solution of hierarchy problem and on the other hand give rich phenomenology for the collider physics . The idea is, that while all the usual particles live in our 3-dimensional world (brane), gravitons are ubiquitous and live in additional dimensions (bulk) also. This explains why the Newtonian constant is so small - gravitons waste their strength smearing it onto the extra dimensions. Potentially, gravitational interaction can be enhanced by Kaluza-Klein excitations of gravitons up to order of electroweak one. How does this mechanism work see, for instance, in the recent reviews . According to concept outlined above the gravitational and electromagnetic interactions can be considered on the equal foot and influence one of them to another might be perceptible. Therefore, we can now speak about electrogravitational processes hoping to observe them at the laboratory. Many such processes are considered already , but in what follows the particular accent will be done on the mass dependence of the corresponding cross sections. Two types of particles (both are yet to be discovered) are related closely to the masses: they are Higgs bosons and gravitons. So, it is not unexpected, that some similarity might be loomed between behaviours of these two particle types. Let us then first recall some features of the Higgs boson production amplitudes.

## 2 Preliminaries

Consider the so called Björken process, $`\mathrm{}^++\mathrm{}^{}Z+`$ higgs $`H^0.`$

In calculation of the cross section for it one usually restricts to the consideration of the graph (a) from Fig.1 only. So do we, but let us take into account masses of the initial fermions (muons, to be definite). Asymptotics of this process at $`\sqrt{s}\mathrm{}`$ is as follows

$`\sigma _{as}^{(c)}(\mu ^+\mu ^{}ZH^0)|_{m_\mu 0}={\displaystyle \frac{2\pi \alpha ^2g_A^2}{\mathrm{sin}^4(2\theta _𝒲)}}{\displaystyle \frac{m_\mu ^2}{m_Z^4}},`$ (1)
where $`g_A^2`$ is the constant of weak axial-vector coupling of the Z boson to muon, $`\theta _𝒲`$ is the angle of electroweak mixing, and all other notations are selfevident. It is easily to see that this result contradicts the unitarity condition. Indeed, the cross-section obtained is azimuthal angle independent. It means, that the scattering process occurs in the s-wave. But the s-wave amplitudes must satisfy the inequality

$$\sigma _{J=0}1/s,$$

while one sees from Eq.(1), that it is constant (equal to $`1.210^2fb`$). In the expression above $`J`$ is the angular momentum of scattering,

In order to avoid the contradiction we include two more amplitudes in the calculation corresponding to the diagrams (b) and (c) of Fig.1. Cross section due to the sum of these two amplitudes is angle independent and is equal to $`1.210^2fb`$ again, but the interference between them and that of (a) results in $`2.410^2fb`$, completely eliminating confusion . This remarkable fine tuning reflects fundamental features of the underlying Standart Model . Namely, it is related to the renormalizibility of the latter, because the asymptotic behaviour of the tree amplitudes reflects the behaviour of the loop integrals with respect to their limits.

## 3 Graviton-photon annihilation of $`\mathrm{}^+\mathrm{}^{}`$ pair

Now, let’s turn to the ”analogous” process $`\mathrm{}^++\mathrm{}^{}\gamma +𝒢`$, to which diagrams on the Fig. 2 correspond.

They are labeled by capital letters: (A) + (B) are Compton like, (C) is the photon exchange and (D) is the contact one, for future references.

In what follows we consider several cases with respect to the mass values of participating particles.

### 3.1 Massless graviton, the lepton mass values are neglected

The cross section is given in the differential form as

$`{\displaystyle \frac{d\sigma }{d\mathrm{cos}\vartheta _\gamma }}={\displaystyle \frac{\pi }{2}}\alpha G_N^{enh}\left(1+\mathrm{cos}^2\vartheta _\gamma \right),`$ (2)
where $`\vartheta _\gamma `$ is the emitting angle of the photon with respect to the negative lepton beam and $`G_N^{enh}`$ is now the constant of graviton coupling to the matter fields enhanced by the virtue of the Kaluza-Klein mechanism.

If, however, we calculate the cross section, corresponding to the diagram (C) of Fig.2 only, we obtain the same result as that of Eq.(2). This is a consequence of the Ward identity. For the illustration, let us write down only the contributions of the Compton like diagrams (A) and (B). Squaring the sum of (A) and (B) amplitudes, we obtain the following expression

$`{\displaystyle \frac{\mathrm{"}d\sigma \mathrm{"}^{Comptonlike}}{d\mathrm{cos}\vartheta _\gamma }}={\displaystyle \frac{5}{2}}\pi \alpha G_N^{enh}`$ (3)
Note, that this expression is constant and scattering angle independent. Therefore, in the full expression for the cross section the interference term appears:

$`2\mathrm{}e\left(_A+_B\right)\left(_C+_D\right)^{}+2\mathrm{}e_C_D^{}.`$ (4)
This term completely compensates the contribution of the Compton like diagrams (A) and (B). Note that after summing over the lepton polarizations the squared (D) diagram vanishes:

$`{\displaystyle \underset{spins}{}}|_D|^2=0.`$ (5)
It is expedient to mention the paper , where the fact of mutual cancellation between contributions of diagrams, describing the gravity in the strong electromagnetic field, was also noted.

### 3.2 Massless graviton, the lepton mass values are not neglected

First, let us calculate that piece of cross section, which corresponds to the photon exchange diagram only in order to see if the case of preceeding subsection is repeated. The result is as follows

$`{\displaystyle \frac{d\sigma ^{(C)}}{d\mathrm{cos}\theta _\gamma }}={\displaystyle \frac{\pi }{2}}\alpha G_N^{enh}\left(2\beta ^2\mathrm{sin}^2\theta _\gamma \right),`$ (6)
while the cross section corresponding to the full set of diagrams is given in the differential form as

$`{\displaystyle \frac{d\sigma }{d\mathrm{cos}\vartheta _\gamma }}={\displaystyle \frac{1}{2}}\pi \alpha G_N^{enh}\beta {\displaystyle \frac{1+2\beta ^2\mathrm{sin}^2\vartheta _\gamma \beta ^4\left(1+\mathrm{sin}^4\vartheta _\gamma \right)}{(1\beta ^2\mathrm{cos}^2\vartheta _\gamma )^2}}\mathrm{sin}^2\vartheta _\gamma ,`$ (7)
where $`\beta =\sqrt{1\frac{4m_{\mathrm{}}^2}{s}}`$, $`m_{\mathrm{}}`$ is the mass of initial lepton.

Evidently, the expression of Eq.(6) and that of Eq.(7) look utterly unalike - the marvellous cancellation of the preceeding subsection disappeared due to presence of lepton masses. Explicitly its behaviour with respect to the scattering angle and to the center of mass energy of collision is depicted on the Fig.3. Being zero at the threshold it growths with energy and quickly goes to plateau. On this picture both variable change in the following limits: $`1cos\vartheta _\gamma +1,`$ and $`2m_\mu 0.212GeV\sqrt{s}3GeV.`$ Units for the cross section values are arbitrary.

Integrating the expression (7) over $`cos\vartheta _\gamma `$ in the limits above we obtain

$`\sigma _{tot}={\displaystyle \frac{4}{3}}\pi \alpha G_N^{enh}\left\{(1+8{\displaystyle \frac{m_{\mathrm{}}^2}{s}})(13)+12{\displaystyle \frac{m_{\mathrm{}}^4}{s^2}}(1+2)\right\},`$ (8)
where $`=\frac{m_{\mathrm{}}^2}{s}\frac{1}{\beta }\mathrm{ln}\frac{1+\beta }{1\beta }`$.

### 3.3 Massive graviton, the lepton mass values are neglected

When the lepton mass is negligible and the graviton mass $`m_𝒢`$ is kept finite, the expression for the cross section looks as

$`{\displaystyle \frac{d\sigma }{d\mathrm{cos}\vartheta _\gamma }}={\displaystyle \frac{\pi }{2}}{\displaystyle \frac{\alpha G_N^{enh}}{(1\frac{m_𝒢^2}{s})^2}}\left\{1({\displaystyle \frac{m_𝒢^2}{s}})^5(1{\displaystyle \frac{m_𝒢^2}{s}})^3\mathrm{cos}^2\vartheta _\gamma \left[3{\displaystyle \frac{m_𝒢^2}{s}}+(1{\displaystyle \frac{m_𝒢^2}{s}})^2\mathrm{cos}^2\vartheta _\gamma \right]\right\}{\displaystyle \frac{1}{\mathrm{sin}^2\vartheta _\gamma }}`$ (9)
The remarkable feature of this expression is the fact, that it has smooth limit at $`m_𝒢0.`$ Namely, in this limit we obtain the expression of the Eq.(2). Thus, the so called van Dam - Veltman - Zakharov discontinuty  is absent in the physical quantity - cross section, although in deriving expressions of Eq.(2) and that of Eq.(9) we have used projection operators for the massless and massive gravitons, for which such a discontinuty presents explicitly.

Now, the cross-section of the Eq.(9) is evidently singular with respect to the scattering angle. This is consequence of neglecting the initial state masses. The analytical expression for the cross section with the initial muon masses taken into account is too cumbersome, so we will treat this case numerically in the next subsection.

### 3.4 Massive graviton, the lepton mass values are not neglected

In this subsection we present the results of the numerical study of the process under investigation with masses of participating particles kept finite ( apart from the zero mass value for the photon, obviously).

Results of this analytical computation of the matrix element squared with the aid of HECAS program  and Monte Carlo phase space integration are very different from that of the massless case and are as follows.

We consider only the one particular case, when $`m_𝒢^2/s1,`$ i.e. $`E_\gamma 0.`$ It appears, that in this case the cross section behaviour is completely determined by the Compton like diagrams only,

$$\sigma ^{alldiagrams}\sigma ^{tot}\sigma ^{\mathrm{"}Comptonlike\mathrm{"}}$$

and there is no infrared cut-off by gravitation, as in the case of massless graviton, i.e.

$`\sigma ^{tot}|_{m_𝒢^2s}\mathrm{}.`$ (10)
At the same time contributions from the photon exchange diagram and the contact one tend to go to the zero in this limit, $`\sigma ^{(C+D)}|_{m_𝒢^2s}0,`$ while the interference between two subsets, namely, between Compton like and the sum of photon exchange and contact one remains constant and negative, $`\sigma ^{(A+B)(C+D)}|_{m_𝒢^2s}const,`$ which , however, has no strong influence on the whole cross section behaviour in view of Eq.(10). So, we were not able to draw any definite conclusion from the case under consideration as opposite to the case of massless graviton. Nevertheless, questions raised in this talk, deserve future study, because not all possibilities are yet considered. To this end we need the fully analytical expression for the cross section in the most general case, i.e. when all masses of the participating particles except photon, are taken into account. We hope to do this job in the nearest future.

## 4 Conclusion

By a simple example of the process $`\mathrm{}^++\mathrm{}^{}\gamma +𝒢`$ the interplay between electromagnetic and gravitational interactions is considered. Some striking features of the amplitudes behaviour are revealed:

* in the case of massles graviton and in the limit of neglecting the initial state mass values the cross section is completely determined by the photon exchange diagram only (diagram (C) of Fig.(2));
* the phenomenon disappears, when masses of initial states are taken into account;
* in the case of massive graviton the behaviour of amplitudes is very different from that considered above.

## 5 Acknoledgments

I would like to thank Prof. Y.S.Kim for invitation me to participate in the work of Symposium and Mrs. Mary Ridgell for excellent job. Also, I am indebted to Drs. O.S.Pavlova and S.K.Stepanyan for useful information they delivered to me and to Prof. D.V.Gal’tsov for making me aware of Ref..


text:
# Structure of Sn isotopes beyond 𝑵=𝟖𝟐

## ACKNOWLEDGMENTS

This work was supported in part by the Italian Ministero dell’Università e della Ricerca Scientifica e Tecnologica (MURST). NI thanks the European Social Fund for financial support.


text:
# Ballistic Spin Injection from Fe(001) into ZnSe and GaAs

## Abstract

We consider the spin injection from Fe into ZnSe and GaAs in the ballistic limit. By means of the ab initio SKKR method we calculate the ground state properties of epitaxial Fe$`|`$ZnSe(001) and Fe$`|`$GaAs(001) heterostructures. Three injection processes are considered: injection of hot electrons and injection of “thermal” electrons with and without an interface barrier. The calculation of the conductance by the Landauer formula shows, that these interfaces act like a nearly ideal spin filter, with spin polarization as high as 99%. This can be traced back to the symmetry of the band structure of Fe for normal incidence.

The controlled injection of a spin polarized current into a semiconductor (SC) is one of the central problems in the new field of spin electronics, since it is a prerequisite for the development of new spin dependent devices Datta . Recently some important successes have been achieved. Fiederling et al. Fiederling  have demonstrated the injection from the paramagnetic II-VI SC BexMnyZn1-x-ySe into GaAs with very high spin polarization using an external magnetic field, while Ohno et al. Ohno  were able to show the injection from the ferromagnetic SC Ga1-xMnxAs into GaAs with an efficiency of 1%. However both methods have the disadvantage that they require a low temperature. Therefore the injection from a ferromagnet with large Curie temperature such as Fe would have strong advantages. Such attempts, though, have not been very successful in the past, i.e., the reported spin injection efficiency was low Hammar ; Filip . Schmidt et al. Schmidt  revealed that a basic obstacle for spin injection from a ferromagnetic metal into a SC exists, being represented by the large conductivity mismatch between both materials. Nevertheless Rashba Rashba  as well as Fert and Jaffrès Fert  have recently pointed out, that this obstacle can be overcome by introducing a tunneling barrier. Meanwhile, and independently, Zhu et al. Zhu  were successful in demonstrating the spin injection at room-temperature from Fe(001) into GaAs with an efficiency of 2% which they attributed to tunneling through a Schottky barrier.

Kirczenow Kirczenow  has lately pointed out that contrary to the ferromagnet$`|`$metal interface the interface between a ferromagnet and a SC could act as an ideal spin filter, if e.g. the Fermi surface of the majority or minority spin bands has a hole at the $`\overline{\mathrm{\Gamma }}`$-point of the two-dimensional Brillouin zone, so that only electrons of the other spin band can scatter into the conduction band states of the SC at the $`\overline{\mathrm{\Gamma }}`$-point. However relevant hybrid systems like Fe$`|`$ZnSe(001) and Fe$`|`$GaAs(001) for which epitaxial growth has been demonstrated, do not show this simple property.

Recently two ballistic calculations Grundler ; Hu  for the spin injection process have been published, which basically rely on a free-electron description of the majority and minority spin bands. Grundler Grundler  could argue in this way that the Fe$`|`$SC interface can act as a spin filter with an efficiency of a few percent. Motivated by this work we present here an ab initio calculation of the ground state properties and the ballistic transport through the Fe$`|`$ZnSe(001) and Fe$`|`$GaAs(001) interfaces. In contrast to the above mentioned methods our calculations include the whole complexity of the band structures of the ferromagnet and the SCs as well as the even more complex properties of the interface. The important result of our calculation is, that the considered Fe$`|`$SC interfaces act like nearly ideal spin filters, with spin injection ratios as high as 99%. We can attribute this to the different symmetries of the majority and minority $`d`$-bands of Fe at the Fermi level, a behavior which cannot be described in the free-electron model. Taken together with the results of Zhu et al. Zhu  our calculations give a bright outlook for the spin injection from ferromagnetic Fe into SCs.

Our method is based on the local density approximation of density functional theory and apply the screened KKR-method Wildberger . The heterostructure consists of a Fe halfspace and a SC (either ZnSe or GaAs) halfspace, both oriented in the (001) direction and being epitaxially bonded at the interface, so that the SC lattice constant is double the Fe constant (a$`{}_{}{}^{\mathrm{exp}}{}_{\mathrm{Fe}}{}^{}`$ = 5.425 a.u. is used in the calculation). The two halfspace Green’s functions are determined by the decimation technique Turek . In the interface region the potentials of 4 monolayers (ML) of Fe and 2 ML of SC are determined selfconsistently. The potentials of all other ML are identified with the asymptotic bulk values. In all calculations we use a cut-off of $`\mathrm{}_{\mathrm{max}}=2`$ for the wavefunctions and an atomic-sphere-approximation for the potentials, but include the full charge density. The ballistic conductance $`G`$ is calculated by the Landauer-Büttiker formalism for $`T=0`$. Here we use an expression similar to the one derived by Baranger and Stone Baranger , but adjusted to the asymptotic Bloch character of the wavefunctions and the two-dimensional translation symmetry of the system. The in-plane component $`𝐤_{}`$ of the $`𝐤`$-vector enumerates then the scattering channels, and we can express the $`𝐤_{}`$-dependent conductance $`G(𝐤_{})`$ wholly in terms of the Green’s function of the system. Spin-orbit coupling is neglected in the calculation.

As we will demonstrate in this paper, the spin injection process is to a large extent determined by the symmetries of the bulk band structures. For this reason we show in Fig. 1 the band structure of Fe and of the SCs ZnSe (a) and GaAs (b), for Bloch vectors $`𝐤=(0,0,k_z)`$ normal to the interface. These are the states relevant for the injection process, since in the SC only states close to the conduction band minimum $`E_C`$ will be populated, having $`𝐤_{}0`$. The left panel shows the spin split majority and minority bands of Fe in the region of the Fermi level $`E_F`$. As usual the different bands in (001) direction are indexed by $`\mathrm{\Delta }_1,\mathrm{\Delta }_2`$, etc. indicating the symmetries of the wavefunctions Callaway . On the right side, the SC bands are shown with $`E_F`$ assumed to be located in the middle of the gap. Most important is here that the lowest conduction states have $`\mathrm{\Delta }_1^{\mathrm{SC}}`$-symmetry; they are invariant under all symmetry operations of the zink-blende lattice, that transform the Bloch vector $`𝐤=(0,0,k_z)`$ in itself. These operations form the symmetry group $`C_{2v}`$, which is at the same time identical with the symmetry of the whole Fe$`|`$SC(001) interface. It is now important to single out those Fe states, which are compatible with this $`C_{2v}`$ symmetry. In Fe, the $`\mathrm{\Delta }`$-nomenclature refers to the $`C_{4v}`$ symmetry group, since, contrary to the zink-blende lattice, in the bcc lattice the (001) direction is a fourfold axis. Thus, not only the $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$-states, consisting locally of $`s,p_z`$ and $`d_{z^2}`$ orbitals, can couple to the $`\mathrm{\Delta }_1^{\mathrm{SC}}`$-band states, but also the $`\mathrm{\Delta }_2^{}^{\mathrm{Fe}}`$-states consisting locally of in-plane $`d_{xy}`$ orbitals. Here we assume that the $`x`$ and $`y`$ directions point along the cubic axes. On the other hand the Fe states of $`\mathrm{\Delta }_2^{\mathrm{Fe}}`$-symmetry (with $`d_{x^2y^2}`$ character) as well as the Fe states with $`\mathrm{\Delta }_5^{\mathrm{Fe}}`$-symmetry (with $`p_x`$ and $`d_{xz}`$ or $`p_y`$ and $`d_{yz}`$ character) cannot couple to the $`\mathrm{\Delta }_1^{\mathrm{SC}}`$-states, since they do not show the full symmetry $`C_{2v}`$ of the heterostructure. For the spin injection it is now important, that in the majority band at $`E_F`$ and above there exists only a $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$-band (below $`E_F`$ also a $`\mathrm{\Delta }_2^{}`$-band is available) while in the minority band around $`E_F`$ only a $`\mathrm{\Delta }_2^{}^{\mathrm{Fe}}`$-band exists that can couple to the $`\mathrm{\Delta }_1^{\mathrm{SC}}`$-states, since the $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$-band appears here at about 1.3 eV above $`E_F`$ (see $`\mathrm{\Gamma }_{12}`$ in Fig. 1).

Not shown in Fig. 1 is the lower $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$-band separated from the upper $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$-band by the so-called $`s`$-$`d`$ hybridization gap. This gap is characteristic for the transition metals and arises from the hybridization of the $`s`$ with the $`d_{z^2}`$ orbitals. For the (001) orientation this gap is so large that for the minority band $`E_F`$ lies in the gap, giving rise to the spin filtering effect discussed in this paper. This effect is also important in magnetic tunnel junctions Mavropoulos .

Firstly we discuss the injection process of hot electrons with Fe states well above $`E_F`$. Although for hot spin injection states with non-zero $`𝐤_{}`$ values also play a role, we consider here for simplicity only states with normal incidence. The calculated transmission probabilities for injection into ZnSe are shown in Fig. 2 for both spin directions, with Fig. 2(a) referring to a Zn terminated interface and Fig. 2(b) to a Se terminated one. The transmission starts at the energy $`E_C`$ of the SC conduction band minimum. In the majority band the conductance strongly increases to values of around 0.6 or 0.7 (in units of $`e^2/h`$), while the conductance in the minority band is much smaller. As a result, the spin polarization of the injected current is very large, for Zn termination always larger than 97%, for Se termination larger than 75%. However, the situation completely changes, if the energy of the injected Fe electrons exceeds the value $`E_{\mathrm{\Gamma }_{12}}`$ of the minimum of the minority $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$-band. There the transmission in the minority band increases very sharply and even overcomes the majority transmission, so that the spin polarization changes sign. This clearly illustrates, that the absence of the $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$-state in the minority band leads for lower energies to the very large spin polarization of the current. Similar results are also obtained for the hot spin injection into GaAs(001), resulting, for lower energies, even in polarizations extremely close to 100%.

The strong spin polarization can be understood from the different spatial orientation and extent of the $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$ and $`\mathrm{\Delta }_2^{}^{\mathrm{Fe}}`$-states. The $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$-states have $`s`$, $`p_z`$ and $`d_{z^2}`$ admixtures. In particular the $`s`$ and $`p_z`$ components have large spatial extent and a strong overlap with the SC states. Moreover the $`d_{z^2}`$ and in particular the $`p_z`$ orbitals point directly into the SC, so that a large transmission is possible. In contrast to this the minority $`\mathrm{\Delta }_2^{}^{\mathrm{Fe}}`$-states consist of in-plane $`d_{xy}`$ orbitals which are much less extended and point in the wrong direction.

To model the injection of electrons at $`E_F`$ we lower the potential in the SC halfspace such that the Fermi level falls slightly above the conduction band energy $`E_C`$. Here we consider two situations, by simulating the injection process both without and with a tunneling barrier. In the first case, we lower the potentials of the 3rd, 4th and all further away SC ML by the same constant value, so that $`E_FE_C`$ becomes positive. We do not change the potentials of the two SC ML closest to the interface, since they are important for the interface characteristics. By continuously varying the potential step, we calculate then the conductance as a function of $`E_FE_C`$. Fig. 3 shows the resulting conductance at the $`\overline{\mathrm{\Gamma }}`$-point for an Fe$`|`$ZnSe(001) junction with Zn interface termination (Fig. 3(a)) and for an Fe$`|`$GaAs(001) junction with Ga termination (Fig. 3(b)). The energy scale in the order of 10 meV refers to typical carrier concentrations in a two-dimensional electron gas Grundler . The insets show the results over a larger energy region. The minority intensities are enhanced by a factor of 10 for ZnSe and by a factor of $`10^4`$ for GaAs. Thus the spin polarizations are larger than 97% for ZnSe and practically 100% for GaAs. Very similar results are also obtained for the other terminations not shown here, i.e. the Se termination of ZnSe and the As termination of GaAs.

All the calculated results in Fig. 2 and Fig. 3 suggest, that near the energy $`E_C`$ of the SC conduction band minimum, the transmission probability varies for both spin directions as $`\sqrt{E_FE_C}`$ broadening . Since the square root-like behavior is the same for both the majority and the minority electrons, the spin polarization remains constant for $`EE_C`$. Moreover, in the interesting energy region of about 10 meV, the reduction of the conductance is rather modest. The square-root like behavior of the transition probability can be understood from a simple picture where a potential step in one dimension is assumed. For a constant potential of height $`V_B`$ in the right halfspace and a vanishing potential in the left halfspace, the transition probability for an incident electron with energy $`E=k^2`$ into a transmitted state with the same energy $`E=V_B+k^2`$ is given by $`T=\frac{4kk^{}}{(k+k^{})^2}4\frac{k^{}}{k}\sqrt{EV_B}`$ for $`k^{}0`$.

To simulate the effect of a Schottky barrier, we modify the above model by smearing out the potential step, i.e. by lowering the external potential continuously over a distance of $`N`$ SC-ML. Within this barrier of $`N`$ ML thickness effectively the Fermi level slowly increases with respect to the local potential, from the ground-state value deep in the gap to an energy value slightly above $`E_C`$. Assuming for this final position of $`E_F`$ a typical energy value $`E_FE_C`$ = 10 meV, we list in Table 1 the resulting spin polarizations $`P`$ at the $`\overline{\mathrm{\Gamma }}`$-point obtained for Fe$`|`$ZnSe(001) and Fe$`|`$GaAs(001) junctions with four different barrier thicknesses of $`N`$ = 8, 32, 80 and 144 ML. As an example for the polarization obtained by integrating over the two-dimensional Brillouin zone, in the last row the polarization is given for a 32 ML thick barrier. Since the integration affects both spin channels approximately equally, the polarization is changed only slightly. While in the case of Se, Ga and As termination the polarization of the spin current is equally high ($`97\%`$) as in the barrier-free case, we see a gradually lowering of the spin polarization for the Zn termination, which however levels off at a value of about 77% for large barrier thicknesses. This effect arises from the existence of minority interface states at the Fe$`|`$SC(001) interface. These states of $`\mathrm{\Delta }_1`$-symmetry lie within the $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$-bulk gap and become resonant due to the coupling with the $`\mathrm{\Delta }_2^{}^{\mathrm{Fe}}`$-band. In the case of Zn-termination, the interface state at $`\overline{\mathrm{\Gamma }}`$ lies relatively close to $`E_F`$, i.e. 0.15eV below. Its effect is to reduce the (positive) spin injection ratio. If this state would coincide with $`E_F`$, its effect would be much bigger and could even lead to a strong negative polarization.

In summary, we have performed ab initio calculations to investigate the ballistic spin injection from a Fe half-crystal into ZnSe and GaAs SCs. Three processes of injection have been considered: the injection of hot electrons as well as the injection of electrons at the Fermi level with and without an interface tunneling barrier. The calculations demonstrate that the Fe$`|`$ZnSe and Fe$`|`$GaAs (001) interfaces act as highly spin-polarizing filters yielding polarizations as high as 99%. This behavior can be traced back to some simple properties of the band structure of Fe for normal incidence: the majority states at the Fermi level have $`\mathrm{\Delta }_1^{\mathrm{Fe}}`$-symmetry and a strong $`s`$ and $`p_z`$ admixture, so that they can couple well to the conduction band states of the SC, while the Fe minority states at $`E_F`$ have a different symmetry and can either couple only weakly or not at all to the SC states. This picture becomes clearer, the more ordered the interface is, since interface disorder breaks the $`𝐤_{}`$ conservation and can reduce the spin polarization of the current. Our calculations and the recent successful observation of a 2% spin injection in the Fe$`|`$GaAs(001) system Zhu  suggest, that much larger spin injection efficiencies should be achievable.

###### Acknowledgements.

The authors thank G. Schmidt for helpful discussions. This work was supported by the RT Network Computational Magnetoelectronics (Contract RTN1-1999-00145) of the European Commission.


text:
# Mechanical detection of FMR spectrum in a normally magnetized YIG disk

thanks: Corresponding author

## Abstract

The ferromagnetic resonance spectrum of a normally magnetized YIG disk, with thickness of 4.75$`\mu `$m and radius of 80$`\mu `$m, is measured at room temperature both by magnetic resonance force microscopy and by standard detection of the microwave susceptibility. The comparison indicates that MRFM represents one of the most potent means of obtaining the *complete* FMR spectra of micron-size samples. In the weak coupling regime, the measured data can be quantitatively understood within the framework of the Damon and Eshbach model.

preprint: To appear in J. Appl. Phys. (manuscript CE-09)

There has been recently a renewed effort, both theoretical Albuquerque et al. (2001) and experimental Hoffmann et al. (1991), to describe the spin dynamics of patterned magnetic thin films in anticipation of new fast magnetic devices. The time scale of interest is the nanosecond because it corresponds to the precession frequency of the magnetization about the internal field direction. The classical method of exciting spin waves is ferromagnetic resonance (FMR). These experiments study the splitting of the ground state through the observation at microwave frequencies of continuous wave absorption spectrum Wigen (1984). But standard cavity methods usually require sample dimensions that are much larger than the typical lateral size of actual micro-fabricated device.

For micron-size samples Zhang et al. (1997), a new sensitive technique, called ferro-Magnetic Resonance Force Microscopy (fMRFM) Zhang et al. (1996), has been developed and its principle is to mechanically detect the net change of dipole moment induced by FMR resonances. It uses a permanent micro-magnet placed in the stray field of the sample, which couples to the longitudinal magnetization (component along the static field, $`M_\textcolor[rgb]{0,0,0}{z}`$) through the dipolar interaction Suh et al. (1998). The force and torque acting on the probe produces a measurable elastic bending of a cantilever onto which the magnet is affixed.

The sensitivity of the mechanical detection will be illustrated with a comparison of MRFM data and standard measurements of the microwave susceptibility. In this paper, the experiments will be restricted to the weak coupling regime, where a large separation ($`h=100\mu `$m) is set between the probe and the sample.

Yttrium Iron Garnet (YIG) is the marvel material of FMR because of its remarkably low loss in magnetic propagation Damon and van de Vaart (1965). An important simplification of the MRFM spectroscopic signal can be obtained by using disk samples, since a probe magnet placed on the disk axis will preserve the axial symmetry. Cylindrical geometries have been extensively studied in the sixties on millimeter-size rods and disks Dillon (1960); Eshbach (1963); Yukawa and Abe (1974). The measured resonances are ascribed to magnetostatic waves propagating radially across the sample. Their separation is determined by the cylinder aspect ratio Dillon (1960). They are labeled by $`(n,m)`$, the number of nodes respectively in the diametrical and circumferential directions. Magneto-exchange modes have been found to be negligible for film thickness above 5$`\mu `$m Wigen (1984).

Our sample is micro-fabricated (ion milled) from a single crystal film (thickness $`S=4.75\mu `$m) oriented along the  direction (easy axis) into a disk of radius $`R80\mu `$m as shown in Fig.1b. The dimensions are large enough so that standard FMR experiments can be carried out on the sample. Fig.1c shows the microwave susceptibility of the disk as a function of the dc magnetic field applied along the disk axis ($`Oz`$). The transverse absorption spectrum ($`M_\textcolor[rgb]{0,0,0}{t}`$) is measured at 10.46GHz by a wide band spectrometer using a standard crystal diode detector. The incident power is 2mW which is about the maximum excitation possible while keeping nonlinear effects reasonable Seagle et al. (1985). Four magnetostatic modes are resolved. They correspond to the longest wavelength spin-wave modes allowed.

The same disk has then been used to test our mechanical FMR detection. The setup is schematically represented in Fig.1a. The MRFM fits between the poles of an iron core electromagnet which applies a static field $`H_{\text{ext}}`$ along the disk axis. A proton NMR gaussmeter is used for the calibration of $`H_{\text{ext}}`$. The microwave field $`H_\textcolor[rgb]{0,0,0}{1}`$ is generated by a 10.46GHz synthesizer and fed into an impedance matched strip-line resonator of length $`5`$mm and width $`0.5`$mm. The YIG disk is placed at the center of this half-wavelength resonator, with the 0.19mm thick GGG substrate intercalated between the YIG and the strip-line. The microwave field $`H_\textcolor[rgb]{0,0,0}{1}`$ can be considered *homogeneous* (within 4%) over the volume of the sample. The sample temperature is fixed at $`T=285`$K and the saturation magnetization is $`4\pi M_\textcolor[rgb]{0,0,0}{s}(T)=1815`$G. A magnetic bar, 18 $`\mu `$m in diameter and 40 $`\mu `$m in length, Malliavin et al. (1999) is glued (see Fig.1b) at the extremity of a soft cantilever (spring constant $`k=0.01`$N/m). SQUID measurements indicate that the room temperature saturation magnetization of our bar is around $`5\times 10^\textcolor[rgb]{0,0,0}{2}`$emu/cm3. The bar is lifted 110$`\mu `$m above the YIG in zero field. A large separation is purposely chosen, so that the coupling between the probe magnet and the sample is weak (the magnetic field gradient produced on the sample is less than $`0.16\text{G}/\mu `$m). This situation is obviously not optimal for imaging purposes, but provides a meaningful comparison with the data measured without the bar. The dipolar field acting on the bar produces an Å-scale displacement measured by a laser beam deflection on a photodiode. When $`H_{\text{ext}}`$ reaches the 0.5 tesla field range, the cantilever bends by 10$`\mu `$m towards the sample surface and the effective spring constant stiffens to $`k=0.2`$N/m. The MRFM signal is proportional to the changes of the longitudinal magnetization $`\mathrm{\Delta }M_\textcolor[rgb]{0,0,0}{z}`$ and thus it increases linearly with microwave power ($`H_\textcolor[rgb]{0,0,0}{1}^\textcolor[rgb]{0,0,0}{2}`$) below saturation (the transverse component $`M_\textcolor[rgb]{0,0,0}{t}`$ being proportional to $`H_\textcolor[rgb]{0,0,0}{1}`$).

Fig.2 shows the field dependence of the mechanical signal when the bar is placed on the symmetry axis of the disk and the amplitude of $`H_\textcolor[rgb]{0,0,0}{1}`$ is fully modulated at the resonance frequency of the cantilever $`f_\textcolor[rgb]{0,0,0}{c}2.8`$kHz. The microwave peak power is increased gradually during the sweep, from 25$`\mu `$W for the longest wavelength modes up to 2.5mW for $`H_{\text{ext}}<4.7`$kOe. The normalized result is shown on a *logarithmic* scale. A series of 50 absorption peak is resolved by our instrument thereby demonstrating the excellent sensitivity of mechanical detection. The linewidth of the peaks is of the order of 1.5G, a typical value for YIG disks Eshbach (1963). In FMR, the linewidth is much smaller than $`\gamma \mathrm{\Delta }H_\textcolor[rgb]{0,0,0}{i}`$, the field distribution inside the sample, because of the propagating character of the spin waves throughout the sample. $`\mathrm{\Delta }H_\textcolor[rgb]{0,0,0}{i}`$ in our disk is set by the dipolar field ($`2\pi M_\textcolor[rgb]{0,0,0}{s}`$) and the additional broadening introduced by the probe magnet is comparatively negligible (less than $`6.2`$G). These findings are in sharp contrast to results obtained with paramagnetic resonance imaging where the excitation is localized to the sheet satisfying the resonance condition which leads to inhomogeneously broaden linewidth.

The theory of magnetostatic modes in thin films has been established by Damon and Eshbach in 1961 Damon and Eshbach (1961) and modified by Damon and van de Vaart Damon and van de Vaart (1965). The transverse mode dispersion relation of forward volume waves can be expressed as:

$$k_\textcolor[rgb]{0,0,0}{t}=\frac{2}{S}\frac{1}{\sqrt{p}}\mathrm{tan}^\textcolor[rgb]{0,0,0}{}\textcolor[rgb]{0,0,0}{1}\left(\frac{1}{\sqrt{p}}\right)$$
(1)
with $`p`$ a parameter that equals to $`p_\textcolor[rgb]{0,0,0}{0}=\left\{B_\textcolor[rgb]{0,0,0}{i}H_\textcolor[rgb]{0,0,0}{i}\omega ^\textcolor[rgb]{0,0,0}{2}/\gamma ^\textcolor[rgb]{0,0,0}{2}\right\}/\left\{\omega ^\textcolor[rgb]{0,0,0}{2}/\gamma ^\textcolor[rgb]{0,0,0}{2}H_\textcolor[rgb]{0,0,0}{i}^\textcolor[rgb]{0,0,0}{2}\right\}`$ in the absence of exchange ($`\omega `$ is the frequency of excitation and $`\gamma `$ is the gyromagnetic ratio). $`p`$ depends on the inhomogeneous internal fields, $`H_\textcolor[rgb]{0,0,0}{i}=H_{\text{ext}}+H_\textcolor[rgb]{0,0,0}{a}4\pi M_\textcolor[rgb]{0,0,0}{s}n_{\textcolor[rgb]{0,0,0}{z}\textcolor[rgb]{0,0,0}{z}}`$ and $`B_\textcolor[rgb]{0,0,0}{i}=H_\textcolor[rgb]{0,0,0}{i}+4\pi M_\textcolor[rgb]{0,0,0}{s}`$, where $`n_{\textcolor[rgb]{0,0,0}{z}\textcolor[rgb]{0,0,0}{z}}`$ is the depolarization factor and $`H_\textcolor[rgb]{0,0,0}{a}=58`$G the magneto-crystalline anisotropy field along the  direction. An analytical expression for $`n_{\textcolor[rgb]{0,0,0}{z}\textcolor[rgb]{0,0,0}{z}}(r,z)`$ exists in the case of uniformly magnetized cylinders Joseph and Schlömann (1965). This expression, however, neglects the rotation of the magnetization close to the edge, a correction of the order $`M_\textcolor[rgb]{0,0,0}{s}/H_{\text{ext}}`$. From Eq.1, we infer that magnetostatic modes are propagating in the region where $`p`$ is positive and are otherwise evanescent waves. One window of special interest (Fig.2a) is when $`4.577\text{kOe}<H_{\text{ext}}<5.375`$kOe which corresponds to $`H_\textcolor[rgb]{0,0,0}{i}<\omega /\gamma <\sqrt{H_\textcolor[rgb]{0,0,0}{i}B_\textcolor[rgb]{0,0,0}{i}}`$ at $`r=0`$. There, the propagating region is the central part of the disk. Excited spin waves experience a force due to the internal field gradient and are accelerated radially towards the center. The wave reflects at the circle $`r=r_\textcolor[rgb]{0,0,0}{1}`$, defined by $`p=\mathrm{}`$ ($`\omega /\gamma =H_\textcolor[rgb]{0,0,0}{i}`$). Standing waves occur when the phase shift over a period $`4_\textcolor[rgb]{0,0,0}{0}^{\textcolor[rgb]{0,0,0}{r}_\textcolor[rgb]{0,0,0}{1}}k(r)𝑑r`$ is equal to $`n\times 2\pi `$ (the integral form incorporates the spatial nonuniformity). In this picture, the internal field distribution affects the phase delay, i.e. the locus of the resonance, but not the linewidth. In Fig.2b, we observe another broad absorption regime for $`H_{\text{ext}}<4.577`$kOe which has been reported by EshbachEshbach (1963) and corresponds to magneto-elastic modes excited by spin waves localized at the edge.

We have used Eq.1 to calculate $`n`$ as a function of $`H_{\text{ext}}`$ for our sample. The resonance condition is calculated on the median plan of the disk ($`z=S/2`$) and the result is shown as a long dashed line on Fig.3. There is no fitting parameter or relative corrections. Our theoretical model, however, assumes a uniform magnetization in the disk and thus underestimates the resonating field of each mode (because it underestimates $`n_{\textcolor[rgb]{0,0,0}{z}\textcolor[rgb]{0,0,0}{z}}`$ Joseph and Schlömann (1965)). Magnetization non-uniformities at the periphery of the disk introduce a diminution of the radial decay of $`H_\textcolor[rgb]{0,0,0}{i}`$ near the disk center. In our picture, this is equivalent to an increase of the disk radius. We have adjusted the radius $`R=85\mu `$m to fit the data at $`n=30`$ (the optimal fitting range for the magnetostatic model). For large $`k`$ ($`n>50`$), Eq.1 needs to be modified to include exchange effects. In our case, this effect can be calculated as a first order correction of the characteristic equation Damon and van de Vaart (1965):

$$k_\textcolor[rgb]{0,0,0}{t}=k_\textcolor[rgb]{0,0,0}{t}(p_\textcolor[rgb]{0,0,0}{0})\frac{D}{\gamma \mathrm{}}\frac{B_\textcolor[rgb]{0,0,0}{i}+H_\textcolor[rgb]{0,0,0}{i}(1+2p_\textcolor[rgb]{0,0,0}{0})}{B_\textcolor[rgb]{0,0,0}{i}H_\textcolor[rgb]{0,0,0}{i}(\omega /\gamma )^\textcolor[rgb]{0,0,0}{2}}\frac{k_\textcolor[rgb]{0,0,0}{t}^\textcolor[rgb]{0,0,0}{2}(p_\textcolor[rgb]{0,0,0}{0})}{S}\left\{\frac{S}{2}k_\textcolor[rgb]{0,0,0}{t}(p_\textcolor[rgb]{0,0,0}{0})+\frac{1}{1+p_\textcolor[rgb]{0,0,0}{0}}\right\}$$
(2)
with $`D=0.93\times 10^{\textcolor[rgb]{0,0,0}{}\textcolor[rgb]{0,0,0}{28}}`$erg.cm2 the exchange parameter. The solid line includes exchange effects and the dashed line in the inset shows the behavior if $`D`$ is omitted. We have also calculated the alteration of the spectrum due to the presence of the magnetic bar. For $`h=100\mu `$m, the correction is small (less than 0.1%), but the changes become more important if the magnet is brought closer to the surface. The variation of the spectrum with $`h`$ will be published elsewhere along with the details of the calculation. For the sake of completeness, the solid line displayed is the full result including the presence of the magnetic probe. The MRFM data agree quantitatively with the model over the full range Damon and Eshbach (1961).

We are greatly indebted to C. Fermon, H. LeGall, O. Acher and A.L. Adenot-Engeluin for their help and support in this work


text:
# Scaling at the selective withdrawal transition

## Abstract

In the selective withdrawal experiment fluid is withdrawn through a tube with its tip suspended a distance $`S`$ above an unperturbed two-fluid interface. At low withdrawal rates, $`Q`$, the interface forms a steady state hump and only the upper fluid is withdrawn. When $`Q`$ is increased (or $`S`$ decreased), the interface undergoes a transition so that the lower fluid is entrained with the upper one, forming a thin steady-state spout. Near this discontinuous transition the hump curvature becomes very large and displays power-law scaling behavior. This scaling is used to show that steady-state profiles for humps at different flow rates and tube heights can all be scaled onto a single similarity profile.

Is it possible to classify topological transitions in nonlinear fluid systems  in the same manner as one classifies thermodynamic transitions? When the topological transition involves formation of a singularity in the fluid flows or interface shapes, a similarity solution can provide a simplified description of the flows and help make such a classification . A crucial component to this approach involves determining how characteristic physical quantities and lengths describing the fluid system scale near the singularity. In many cases these singularities manifest themselves in the transition dynamics  and do not appear in the steady state flows. Here, we report on steady-state interface profiles near the topological transition associated with the selective withdrawal experiment. Despite the transition being discontinuous, scaling of the interface is observed as the transition is approached.

In the selective withdrawal experiment a tube is immersed in a filled container so that its tip is suspended a height $`S`$ above an unperturbed interface separating two immiscible fluids. When fluid is pumped out through the tube at low flow rates, $`Q`$, only the upper fluid is withdrawn and the interface is deformed into an axi-symmetric steady-state hump (Fig. 1) due to the flows in the upper fluid. The hump grows in height and curvature as $`Q`$ increases or $`S`$ decreases until the flows undergo a transition where the lower fluid becomes entrained in a thin axi-symmetric spout along with the upper fluid. The two-fluid interface becomes unbounded in the vertical direction thus changing the topology of the steady state. Once the spout has formed, an increase in $`Q`$ or decrease in $`S`$ thickens the spout.

The interfacial profiles at different flow rates and tube heights are recorded. Near the transition, the steady-state radius of curvature of the hump tip is orders of magnitude smaller than the length scales characterizing the boundary conditions (e.g. the tube diameter, $`D`$). This separation of length scales suggests that a similarity analysis of the steady-state hump profiles might be possible. However, for the range of parameters explored thus far, even when the system is arbitrarily close to the transition from hump to spout, the mean curvature of the hump tip, $`\kappa `$, while large, remains finite. Nevertheless, by fixing $`S`$ and looking at the steady-state profiles as $`Q`$ is increased, we observe that both the hump curvature and height display scaling behavior characteristic of systems approaching a singularity. Since the divergence is cut off before a singularity is reached this transition appears to be “weakly-first-order.”

As shown in Fig. 1, the parameters important for this experiment are the upper and lower fluid viscosities and densities ($`\eta _u,\eta _l,\rho _u,\rho _l`$), interfacial tension ($`\gamma `$), orifice diameter ($`D`$), tube height ($`S`$), and flow rate ($`Q`$). In looking for scaling of the steady-state profiles, care must be taken to design an experimental apparatus capable of isolating the profiles near the transition. Experiments were performed in large tanks (30cm $`\times `$ 30cm $`\times `$ 30cm) capable of holding fluid layers that were each about 12 cm in height. To ensure that the upper fluid level remained constant, the withdrawn fluid was recycled back into the container (the bottom fluid layer thickness remains constant when the system is in the hump state). Steady withdrawal rates were achieved by using a gear pump. We verified that for the tube diameter ($`D=0.16`$ cm), tube heights ($`0.1`$ cm $`S1.1`$ cm), and flow rates ($`Q10`$ ml/sec) used in the experiments, the container walls were sufficiently distant and the fluid layers sufficiently thick so as not to affect the flows . We measured the upper(Heavy Mineral Oil) and lower(Glycerin-Water mixture) fluid viscosities to be $`\eta _u=2.29`$ St and $`\eta _l=1.90`$ St, the upper and lower fluid densities to be $`\rho _u=0.88`$ g/ml and $`\rho _l=1.24`$ g/ml, and the surface tension to be $`\gamma =31`$ dynes/cm. Attempts to increase $`\kappa `$ near the transition by decreasing $`\gamma `$ cause fluid mixing and result in a diffuse interface at high shear rates so that some fraction of the lower fluid is always being withdrawn .

While many of the parameters mentioned influence the flows, our understanding of the scaling behavior can be conveyed by focusing on $`S`$ and $`Q`$. We can fix $`Q`$ and track the development of the hump profiles as a function of $`S`$. Below the tube height, $`S_u`$, the hump is unstable and undergoes a transition to a spout. Figure 2 shows that $`S_uQ^{0.30\pm 0.05}`$ . At low $`Q`$ the transition is hysteretic: the value of $`S`$ where the spout becomes unstable and decays back into the hump is larger than $`S_u`$. We define the difference of the two heights or hysteresis as $`\mathrm{\Delta }S`$. Figure 2 indicates that the data is consistent with an exponential decrease: $`\mathrm{\Delta }S=0.04\mathrm{exp}^{Q/0.032}`$. For $`Q>0.1`$ ml/sec, $`\mathrm{\Delta }S`$ was too small to measure.

In order to measure the mean curvature of the hump tip, $`\kappa `$, we first fit the tip of the recorded profile with a Gaussian and then calculate the curvature of the fitting function at the hump tip. Figure 2 also shows the hump height, $`h_u`$, and mean radius of curvature, $`1/\kappa _u`$, at the transition as a function of $`Q`$. The dramatic decrease in $`\mathrm{\Delta }S`$ coincides with the onset of a flat asymptotic dependence for $`1/\kappa _u`$ at $`Q>0.1`$ ml/sec. We quantify this correlation, by fitting the curvature data with the form $`1/\kappa _u=0.02+0.32\mathrm{exp}^{Q/0.032}`$ which has the same exponential decay with $`Q`$ as does $`\mathrm{\Delta }S`$. For $`Q>0.1`$, ml/sec we find both $`\kappa _u`$, and $`h_u`$ to be independent of the orifice diameter $`D`$ . We restrict our scaling analysis to this regime.

Figure 3a plots $`\kappa `$ (where $`0\kappa \kappa _u`$), versus $`Q`$ for six data sets corresponding to different values of $`S`$. As shown in the inset of Fig. 3a, all fifteen data sets display a power-law divergence for $`\kappa `$, as $`Q`$ approaches $`Q_c`$ (a fitting parameter). While the power-law exponents remain constant as $`S`$ is varied, the prefactors to the power laws, $`c_\kappa (S)`$, vary slightly with $`S`$ and are scaled out in the inset. Figure 3b plots the hump height, $`h_{max}`$, versus $`Q`$. The inset to Fig. 3b shows that as $`Q`$ approaches $`Q_c`$ (obtained from Fig. 3a inset), the hump height approaches $`h_c`$ (a fitting parameter) as a power law. Once again, the power-law exponents for these data sets are the same over this range of $`S`$. The prefactors, $`c_h(S)`$, are scaled out in the inset. Note that $`Q_c`$ changes with $`S`$ indicating that the system can approach a continuous line of divergences. Combining the two scaling dependencies in Fig. 3c, we plot $`(h_ch_{max})/h_{max}`$ versus the normalized curvature, $`\kappa /n`$. We find that $`(h_ch_{max})/h_{max}`$ scales as $`(\kappa /n)^{0.85\pm 0.09}`$ indicating that even though both $`h_c`$ and the power-law prefactor, $`n`$, change with $`S`$, the power-law exponents are independent of $`S`$ for this range of tube heights. Note that $`n(S)=c_h(S)[c_\kappa (S)^{0.85}]`$. The transition cuts off the evolution of the hump states making it impossible for the system to approach arbitrarily close to the singularity and limiting the precision with which we can determine the exponents.

The scaling observed for $`h_{max}`$ and $`\kappa `$ suggests that the hump profiles should display universal behavior as $`h_{max}`$ nears $`h_c`$. The quantities $`1/(\kappa /n)`$ and $`(h_ch_{max})/h_{max}`$ track how quickly the radial and axial length scales decrease as the system approaches the singularity and are therefore used to scale the profiles. We define the scaled variables:

$$H(R)=\frac{h_ch(r)}{h_ch_{max}}andR=\frac{r\kappa }{n},$$
(1)
where $`h(r)`$ is the hump profile and $`h_c`$ is taken from Fig. 3. The transformation shifts the profiles so that under scaling the singularity occurs at the origin and the maximum hump heights occurs at $`H=1`$ and $`R=0`$. Figure 4 shows eight scaled profiles for the $`S=0.830`$ cm data set. In the bottom inset we overlay the hump profiles scaled in the main figure. We find excellent collapse of the profiles. The solid line in Fig. 4 is a power law that fits the data in the region beyond the parabolic tip. The picture that emerges is of a parabolic tip region which decreases in its radial scale and is simultaneously pulled towards the singularity in the axial direction leaving behind it a power-law profile with an exponent of $`0.72\pm 0.08`$. This exponent is within error (although slightly smaller) of the exponent observed in the scaling relation of Fig. 3c which predicts an exponent of $`0.85\pm 0.09`$.

Typically, the observed scaling dependencies in these types of problems result from the local stress balance. A scaling analysis where the viscous stresses of the upper and lower fluids balance the stress arising from the interfacial curvature predicts linear scaling dependencies and conical profile shapes. The non-linearity of the observed dependencies indicates that either a different stress balance governs the flows (e.g. only viscous stress due to upper fluid balances stress due to the interface curvature) or that non-local effects are coupling into the solution. A more detailed discussion can be found in .

Finally, we compare the similarity curves for five different tube heights in the upper right inset of Fig. 4. The profiles all display the same power-law dependence. Within error, the normalized curvature $`\kappa /n`$ (taken from Fig. 3c) can be used to scale the radial components of these profiles and obtain good collapse. In Fig. 3c we find that the normalization prefactors, $`n`$, decrease as $`\mathrm{exp}^{2.5S}`$. Here, we correlate this decrease with the observation that the profiles become shallower at larger $`S`$. The points of deviation for the $`S=0.255`$ cm and $`0.381`$ cm profiles mark the transition from the similarity regime to the matching regime beyond which the profiles become asymptotically flat. At large enough radii all of the scaled profiles display these deviations.

We have shown that in the $`Q>0.1`$, ml/sec regime, a similarity analysis can be used to describe the flows near the selective-withdrawal transition . We have observed power-law scaling of the hump height and curvature (Fig. 3) and used these scaling relations to collapse the hump profiles at different flow rates and tube heights onto a single universal curve (Fig. 4). However, the origin of the saturation of $`\kappa `$ at large $`Q`$ remains an important unexplained problem. Further insight into this cutoff behavior may be gained by comparing with an analogous two-dimensional (2-D) problem which roughly corresponds to replacing the tube with a line sink. Jeong and Moffatt  showed that in an idealized case where the bottom fluid is inviscid while the top fluid is very viscous, the 2-D hump interface forms a 2-D cusp singularity above sufficiently high withdrawal rates. Recently, Eggers  showed that the solution changes when the lower fluid has a finite viscosity; the system no longer manifests a singularity . Instead, the approach to the singularity is cut off and the system undergoes a transition to a different steady state. In this new state, a sheet of the lower fluid is entrained along with the upper fluid into the line sink. However, the finite lower fluid viscosity prevents the hump profiles from scaling onto a similarity solution.

Here, we have shown that for the three-dimensional selective withdrawal system even when both fluids are viscous ($`\eta >1`$ St for both fluids) the effects of a singularity manifest themselves in the scaling of the hump profiles. Furthermore, preliminary experiments show that a reduction of the lower fluid viscosity to 0.01 St has little effect in determining the final curvature of the hump tip or equivalently, how close the system is to forming a cusp. This suggests that for our 3-D problem, either the effects of the lower fluid viscosity enter as a higher-order perturbation to the profile shapes, or a different mechanism underlies the avoidance of the cusp formation. If the latter scenario is correct, it may be possible for the system to manifest the singularity at a finite lower fluid viscosity. In either case, determining which variables affect how close the system approaches the singularity would allow for control of the maximum hump curvature and minimum spout diameter. This control could then be used to advance technologies such as coating microparticles , creating mono-dispersed micro-spheres , and emulsification through tip streaming  which take advantage of the selective withdrawal geometry.

We are grateful to W. W. Zhang, S. Venkataramani, J. Eggers, H. A. Stone, T. J. Singler, J. N. Israelachvili, C. C. Park, S. Chaieb, S. N. Coppersmith, T. A. Witten, L. P. Kadanoff, R. Rosner, P. Constantin, R. Scott, T. Dupont H. Diamant, and V. C. Prabhakar for sharing their insights. This research is supported by the University of Chicago (MRSEC) NSF DMR-0089081 grant