text:
# Formation of an Evanescent Proto-Neutron Star Binary and the Origin of Pulsar Kicks

## 1. Introduction

It is recognized that radio pulsars have peculiar space velocities between $`30\mathrm{km}\mathrm{s}^1`$ and $`1600\mathrm{km}\mathrm{s}^1`$, significantly greater than those of their progenitor stars. The highest speed ever recorded is from the Guitar Nebula pulsar (Cordes et al. 1993; Cordes & Chernoff 1998), while the lowest has been recently measured for B2016+28 by Brisken et al. (2002) from very accurate VLBA pulsar parallaxes. Statistical studies, aimed at inferring the peculiar velocity at birth from the observed speed, give mean three-dimensional velocities of 100 – 500 $`\mathrm{km}\mathrm{s}^1`$ for the isolated pulsars (Lyne & Lorimer 1994; Lorimer, Bailes, & Harrison 1997; Hansen & Phinney 1997; Cordes & Chernoff 1998; Arzoumanian, Chernoff, & Cordes 2002).

An early explanation for the large space velocities called for recoil in a close binary that becomes unbound at the time of (symmetric) supernova explosion (Blaauw 1961; Iben & Tutukov 1996). Now, a number of observations hint at a natal origin of these high space velocities, or at a combination of orbital disruption (when the progenitor lives in a binary) and internal kick reaction. Evidence for a kick at the time of neutron star birth is now found in a variety of systems: in runaway O/B associations (Leonard & Dewey 1992), in highly eccentric Be/NS binaries (van den Heuvel & Rappaport 1986; Portgies-Zwart & Verbunt 1996), in the binary pulsar J0045-7319 (Kaspi et al. 1996; to explain its current spin-orbit configuration), and in double neutron star binaries such as B1913+16 (Bailes 1988; Weisberg, Romani, & Taylor 1989; Cordes, Wasserman, & Blaskiewicz 1990; Kramer 1998; Wex, Kalogera & Kramer 2000) where misalignment between the spin and orbital angular momentum axes indicates velocity asymmetry in the last supernova. All these observations support the view that the formation of neutron stars is accompanied by anisotropic explosion. This notion is bolstered by evolutionary studies of binary populations (Dewey & Cordes 1987; Fryer & Kalogera 1997; Fryer, Burrows, & Benz 1998) and by studies (e.g. Cordes & Wasserman 1984) on the survival of binary systems into their late evolutionary stages after supernova explosion.

Among the physical processes that have been proposed to account for the kicks are large-scale density asymmetries seeded in the pre-supernova core (leading to anisotropic shock propagation), asymmetric neutrino emission in presence of ultra-strong magnetic fields (see Lai, Chernoff, & Cordes 2001 for a review), or off-centered electromagnetic dipole emission from the young pulsar (Harrison & Tademaru 1975). However, none of these mechanisms can explain kicks as large as $`1600\mathrm{km}\mathrm{s}^1`$.

In this paper we reconsider the idea first put forward by Imshennik & Popov (1998) that, in the collapse of a rotating core, one or more self-gravitating lumps of neutronized matter may form in close orbit around the central nascent neutron star, transfer mass in the short lived binary, and ultimately explode, causing the remaining, massive neutron star to acquire a substantial kick velocity, as high as the highest observed. The light member explodes as mass transfer drives it below the minimum stable mass for a neutron stars. In the light star, stability is lost upon decompression by the $`\beta `$decaying neutrons and nuclear fissions by radio-active neutron-rich nuclei (Colpi, Shapiro & Teukolsky 1989, 1991; Blinnikov et al. 1990), that deposit energy driving matter into rapid expansion (Colpi, Shapiro & Teukolsky 1993; Sumiyoshi et al. 1998). The kick tries its origin from the orbital motion of this evanescent super-close binary, that forms in the collapse of a rapidly rotating (isolated) iron core.

We will study several effects that may modify the magnitude of the kick, such as gravitational bending of the exploding debris, rotational averaging of the momentum impulse, orbit decay, and delayed neutron star cooling.

Formation of a proto-neutron star companion around the main neutron star has never been verified in numerical simulations due to computational limitations. For this reason we elaborate on a study of Bonnell (1994) on the formation of binary/multiple systems in collapsing gas cloud cores and its extension to the stellar core collapse in the aftermath of a supernova explosion (Bonnell & Pringle 1995) to motivate our working hypothesis. The formation of such exotic binaries has been conjectured to occur in many works (Ruffini & Wheeler 1971; Clark & Eardley 1977; Blinnikov, Novikov, Perevodchikova, & Polnarev 1984; Nakamura & Fukugita 1989; Stella & Treves 1987).

## 2. Light fragments around proto-neutron stars

### 2.1. The Scenario

Formation of a light companion around a main body implies breaking of spherical and axial symmetry during collapse and following core bounce. During dynamical collapse, unstable bar modes ($`m=2`$) can grow in a fluid (even non rotating) that may end with fragmentation. However this is known to occur only if the cloud core contracts almost isothermally, as in the case of star’s formation from unstable cold gas clouds (Bonnell 1994). Core collapse in type II supernovae is far from isothermal (it is described by a effective polytropic index $`\gamma 1.3`$) so that instabilities of this type do not have time to grow (Lai 2000) 111Large-scale asymmetries imprinted in the iron core prior to collapse may lead to anisotropic explosions that produce kicks as indicated by Goldreich, Lai & Sahrling (1996). Whether they lead to fragmentation is unknown., and simulations of non-axisymmetric rotating core collapse confirm this trend (Rampp, Muller & Ruffer 1998; Centrella, New, Lowe, & Brown, 2001). Can fragmentation/fission be excited after core bounce ?

Rapid rotation in equilibrium bodies is known to excite non-axisymmetric dynamical instabilities and these instabilities may grow in the proto-neutron star core. Interestingly, core collapse simulations of unstable rotating iron cores (Heger, Langer & Woosley 2000; Fryer & Heger 2000) or polytropes (Zwerger & Muller 1997) indicate that proto-neutron stars, soon after formation, can rotate differentially above the dynamical stability limit set when the rotational to gravitational potential energy ratio $`T_{\mathrm{rot}}/|W|`$ is larger than the value $`\beta _{\mathrm{dyn}}=0.250.26`$ (Saijo et al. 2001). Strong non-linear growth of the dominant bar-like deformation ($`m=2`$) is seen in these cores (described as polytropes by Rampp, Muller, & Ruffert 1998). However there is no sign of fission into separate condensations. The bar evolves, producing two spiral arms that transport the core’s excess angular momentum outwards (see also Shibata, Baumgarte & Shapiro 2000). This reduces the bar’s angular momentum such that a single central body is formed. How can a body develop local condensations when bar-unstable ?

According to Bonnell’s picture, the evolution of the bar instability is more complex, in reality. If the rapidly spinning proto-neutron star core goes bar unstable when surrounded by a fall-back disk, then matter present in the bar-driven spiral arms interacts with this material. The sweeping of a spiral arm into fall-back gas can gather sufficient matter to condense into a fragment of neutronized matter. This occurs because the $`m=1`$ mode grows during the development of the $`m=2`$ mode. The $`m=1`$ mode causes the displacement of the unpinned (free-to-move) core and creates an off-center spiral arm that sweeps up more material on one side than the other during ”continuing” accretion. The condensation eventually collapses into a low mass neutron star.

Detailed simulations confirming or dismissing the occurrence of such instability are still lacking, so further considerations of the lump masses, temperatures and entropy contents are necessarily speculative.

### 2.2. Cooling scenarios and the minimum mass

Here we wish to explore the possibility that a light (proto-)neutron star forms around the main central body, lives for a while, and later explodes, imprinting a kick to the neutron star that remains (due to linear momentum conservation) before gravitational waves or hydrodynamical effects can drive it toward the central star. Can a light (proto)-neutron star form from the condensation of material accumulated in the off-centered spiral arm ? What limits can be imposed on its mass ?

Cooling plays a key role in addressing these questions. Goussard, Haensel & Zdunik (1998), Strobel, Schaab & Weigel (1999), and Strobel & Weigel (2001) have shown that the value of the minimum stable mass $`m_{\mathrm{mmc}}`$ for a neutron star (located at the turning point in the mass-radius relation of equilibria) is a function of the temperature: it varies from $`\mathrm{¿}\mathrm{}1M_{}`$ at 50-100 milliseconds after core bounce (setting the actual value of the mass of the central neutron star), to $`0.7M_{}`$ after $`1`$ second, down to $`0.3M_{}`$ after 30 seconds, reaching the value of $`0.0925M_{}`$ (Baym, Pethick & Sutherland 1971) known for cold catalyzed matter ($`T<1`$ MeV) several hundred seconds later. Thus, the lump of nuclear matter gathered in the spiral arm by the instability may become self-bound if its mass $`m`$ is above the minimum corresponding to that particular temperature. Once formed it is stabilized against expansion by cooling. The actual value of $`m`$ is thus determined by the overall dynamics of collapse after core bounce and varies between 0.0925 and $`1M_{}:`$ It depends on the time at which the instability sets in, on the amount of fall-back material (potential reservoir of matter in the lump) and on the cooling history.

The value of $`m`$ and of the mass ratio $`q=m/M`$ in the binary (with $`M`$ the heavier of the two stars) remains unpredictable to us at this level, so we may just depict three possible scenarios for the formation and evolution of this evanescent binary: case (A) when the instability sets after few hundreds of milliseconds or a second after core bounce and binary formation occurs so that $`m\mathrm{¿}\mathrm{}0.7M_{}`$ (implying a pre-existing iron core of large mass if the primary is as massive as 1.4$`M_{}`$); case (B) when a lighter star forms around a main body with $`m0.20.4M_{}`$ after several tens of seconds from core bounce; case (C) when cooling is sufficiently advanced that the minimum mass approaches its asymptotic value and the binary can have $`m\mathrm{¡}\mathrm{}0.2M_{}`$.

The magnitude of a natal kick in (A) is difficult to estimate and we must await for realistic simulations of core collapse. Values of $`m/M`$ above a given threshold in a binary are known to lead to unstable mass transfer and thus to final coalescence. We note that the evolution in this case might be similar to that described in coalescing neutron star binaries with slightly unequal masses (see Rosswog et al. 2000) where it has been shown that large kicks can be acquired as a consequence of mass loss via a wind. Case (B) and (C) will be explored in this paper. Imshennik & Popov (1998) estimated neutron star recoil speeds of order $`15002000\mathrm{km}\mathrm{s}^1`$ for scenario (C).

In cases (B) and (C), one factor that could lower the final neutron star speed is the finite velocity of the ejecta. Colpi, Shapiro & Teukolsky (1993) have shown that in the dynamical phase of the explosion the ejecta can attain speeds varying from 10,000 to $`50,000\mathrm{km}\mathrm{s}^1`$ (the upper bound being related to the value of the binding energy of the star at the minimum mass relative to a dispersed state of iron). These speeds are close to the escape velocity from the binary and therefore the final kick imparted to the remaining neutron star may be influenced substantially by gravitational deflection of the ejecta, and velocity phase averaging during the explosion. A further threat would be a rapid decay of the orbital separation due to unstable mass transfer that may lead to final coalescence, and emission of gravitational waves.

In $`\mathrm{\S }3`$ we will describe the gravitational bending of the ejecta, while in $`\mathrm{\S }4`$ we will study initial conditions and, subsequently, orbit decay and mass exchange, for cases (B) and (C). We will then give an estimate of the kick speed including phase-velocity averaging.

## 3. Gravitational bending of the exploding debris and the final kick

Consider now the instant at which the light secondary of mass $`m`$ explodes, having reached the dynamical instability point ($`mm_{\mathrm{mmc}}`$), while orbiting the heavier, primary neutron star of mass $`M`$.

Obtaining the final kick speed imparted to the heavy neutron star is complicated by several factors. First, some of the material ejected from the exploding star can remain bound to the system and eventually accretes onto the remaining star. Second, material that escapes may not be moving at extremely large speeds, even right at ejection, and therefore emerges with a different momentum than it has at the point of explosion. Third, the mass of ejected material may be of order 10% of the total mass of the system, and self gravity could play a role. Fourth, the ejecta are a fluid, and, depending on how much time elapses between neutron star formation and the explosion of the low mass star, move through ambient gas left over from the original supernova explosion. If the ejecta are slowed significantly by the surrounding gas, some material that would be judged unbound upon ejection might actually end up falling back toward the neutron star, possibly to accrete onto it.

Suppose that relative to the remaining heavy neutron star, the orbital velocity of the exploding star is $`𝑽`$ at the instant of explosion ($`t=0`$), and the position of the center of mass of the exploding star is at $`𝒓.`$ Let us write the final velocity of the remaining neutron star as

$$𝑽_{\mathrm{kick}}=𝜼\frac{mV}{(M+m)},$$
(1)
where $`V=|𝑽|`$. If the ejecta are expelled isotropically relative to $`𝒓`$, and at speeds large compared with $`V`$, then we expect $`\eta =|𝜼|1`$, in which case the remaining neutron star recoils with the maximum possible kick speed

$$V_{\mathrm{kick},\mathrm{max}}\frac{mV}{(m+M)}$$
(2)
as implied by momentum conservation in the center of mass of the binary.

To get a more realistic, but still approximate, idea of the size of $`\eta =V_{\mathrm{kick}}/V_{\mathrm{kick},\mathrm{max}}`$, we work to zeroth order in the ejected mass, and also treat the trajectories of the ejecta ballistically; this is equivalent to ignoring the third and fourth complications listed above entirely. By working to zeroth order in the mass of the ejecta, we can also assume that all particles are ejected virtually from a single point, $`𝒓`$. Let the velocity of a particle relative to the exploding star be $`𝒘`$, and assume that $`𝒘`$ is distributed isotropically about the position of the exploding star with probability $`P(𝒘)`$. Since we work to zeroth order in the mass of the ejecta, we may take the heavy neutron star to remain at fixed position in the calculation.

The problem that remains is the orbital mechanics for each particle relative to the massive star. The orbits are characterized by the conserved quantities

$`E`$ $`=`$ $`{\displaystyle \frac{1}{2}}|𝑽+𝒘|^2{\displaystyle \frac{GM}{r}}`$
$`𝑱`$ $`=`$ $`𝒓\mathbf{\times }(𝑽+𝒘)`$
$`𝑨`$ $`=`$ $`{\displaystyle \frac{GM𝒓}{r}}+(𝑽+𝒘)\mathbf{\times }𝑱,`$ (3)
which are the orbital energy and angular momentum per unit mass, and the Runge-Lenz vector, respectively. We note that

$$A=|𝑨|=GMe=\sqrt{(GM)^2+2EJ^2};$$
(4)
$`e<1`$ for bound orbits, $`e>1`$ for unbound orbits.

To calculate $`\eta `$, we need to find the momentum per unit mass carried away by the ejecta that escape to infinite distance. First, we need to identify the initial velocities of particles that escape. For these, we must have $`E>0`$, and from equation (3) we find the condition

$$𝑽\mathbf{}𝒘=Vw\mu >\frac{GM}{r}\frac{(V^2+w^2)}{2}$$
(5)
where $`V=|𝑽|`$ and $`w=|𝒘|`$, and $`\mu =𝑽\mathbf{}𝒘/Vw`$. Second, we need to find the mapping between initial and final velocity. To this aim, it is useful to define a coordinate system in which $`\widehat{𝒆}_\mathrm{𝟑}=\widehat{𝑱},\widehat{𝒆}_\mathrm{𝟏}=\widehat{𝑨},`$ and $`\widehat{𝒆}_\mathrm{𝟐}=\widehat{𝑱}\mathbf{\times }\widehat{𝑨},`$ where $`\widehat{𝑱}=𝑱/J`$ and $`\widehat{𝑨}=𝑨/A`$. The orbit is then confined to the $`12`$ plane, and we can read off the final velocity in this system from equation (15.12) in Landau and Lifshitz, Mechanics:

$$𝒗_{\mathrm{}}=[𝑨\sqrt{2E}+(2E/GM)𝑱\mathbf{\times }𝑨].$$
(6)
We wish to express the final velocity $`𝒗_{\mathrm{}}`$ in a coordinate system fixed in the binary at the point of explosion. We thus choose, explicitly

$$𝒓=r\widehat{𝒙}\mathrm{and}𝑽=\mathrm{V}\widehat{𝒚}$$
(7)
so that the orbital angular momentum of the light star relative to $`M`$ just prior to explosion points along $`\widehat{𝒛}`$. We have explicitly selected the point of explosion to occur at either peri- or apocenter of the orbit, and for simplicity, we specialize to circular orbits below. (Both tidal effects and gravitational radiation will tend to circularize the orbits, but, even for eccentric orbits, we expect tidal disruption to be likeliest at pericenter.) In this coordinate system,

$$𝑱=\widehat{𝒚}rw_z+\widehat{𝒛}r(V+w_y)$$
(8)
$`𝑨=\widehat{𝒙}[GM+r(V+w_y)^2+rw_z^2]`$
$`+\widehat{𝒚}[rw_x(V+w_y)]+\widehat{𝒛}(rw_xw_z)`$ (9)
and

$$E=\frac{GM}{r}+\frac{w_x^2+(V+w_y)^2+w_z^2}{2}.$$
(10)
The velocity of the ejecta can thus be expressed in terms of the known explosion parameters ($`𝒘,V`$, and $`r`$) by combining equation (6) with equations (7)-(10). The result is that the total momentum carried off per unit mass of ejecta is

$$𝒖=_{E(𝒘)>0}d^3𝒘P(𝒘)𝒗_{\mathrm{}}(𝒘),$$
(11)
so that $`\eta =|𝒖|/V,`$ to zeroth order in $`m/M.`$ As $`P(𝒘)=P(w)`$ is isotropic we can separate out the integrals

$$𝒖=𝑑ww^2P(w)d^2\widehat{𝒘}𝒗_{\mathrm{}}(w\widehat{𝒘}),$$
(12)
and consider different distributions of ejection speeds separately. The integration only extends over unbound orbits. For $`w>\sqrt{2GM/r}+V`$, all ejecta escape, and for $`w<\sqrt{2GM/r}V`$, no ejecta can escape.

Because of reflection symmetry with respect to the binary orbital plane, $`u_z`$ vanishes identically. The kick is thus given in the orbital plane of the pre-explosion binary. It is easy to show that when the particles are ejected with equal velocity $`w_o`$ so that $`P(w)=\delta (ww_o)/(4\pi w_o),`$ the integrals yield $`u_x(GM/r)w_o^1`$ and $`u_yV`$ for very large values of $`w_o`$ ($`V`$) which is consistent with the requirement that at large ejection speeds, the outflow reaches infinity with an average velocity comparable to that of the exploding star. Figure 1 shows $`\eta `$ as a function of the dimensionless explosion velocity $`w_o/\sqrt{GM/r}.`$ Notice that there is no kick (i.e., $`\eta `$ vanishes) when $`E=0`$ (see eq. ). This occurs at a critical value of the expansion speed $`w_o/\sqrt{GM/r}=\sqrt{2}V/\sqrt{GM/r}`$ which is $`0.414`$ for $`V\sqrt{GM/r}`$ to zeroth order in $`m/M.`$ When $`E0`$, the mass outflow vanishes identically and the ejecta fall back to the remaining neutron star transferring all their linear momentum to it.

## 4. Final Kick Speed from the Short-Lived Binary

### 4.1. Masses and orbital separations

The final kick speed imparted to the remaining neutron star depends crucially on the evolutionary scenario leading to the formation and disruption of its low mass companion. The range of possible values of $`m_{\mathrm{in}},`$ the initial mass of the secondary (hereon), depends on the time the fragment forms and how fast its cools. Both of these complications distinguish the problem of “early” formation and decay of an evanescent proto-neutron star binary (case B) from the problem of “late” decay of a neutron star binary (case C) originally discussed by Clark & Eardley (1977), Blinnikov et al. (1984), and Imshennik & Popov (1998); the final outcome may differ from the results found by the authors for cold neutron stars.

Figure 2 shows an estimate of the binary separation $`r_{\mathrm{orb}}`$ at the time of formation (solid line) as a function of $`m_{\mathrm{in}}`$, which we regard as an independent variable here. To obtain this estimate, we assumed that the angular momentum $`J_{\beta ,\mathrm{dyn}}`$ of the bar-unstable proto-neutron star with $`T_{\mathrm{rot}}/|W|\beta _{\mathrm{dyn}}`$ goes entirely into orbital angular momentum of the binary. (Even if the lump forms in corotation, its angular momentum is smaller than the orbital angular momentum of the system, substantially so if the orbital separation is considerably larger than its radius.) As a reference value for $`J_{\beta ,\mathrm{dyn}}`$ we adopt $`4\times 10^{49}\mathrm{g}\mathrm{cm}^2\mathrm{s}^1`$, the value given by Saijo et al. (2000) for differentially rotating compact stars at $`\beta _{\mathrm{dyn}}`$ (note that this is the minimum $`J_\beta `$ since the $`m=2`$ and $`m=1`$ modes can grow at values of $`T_{\mathrm{rot}}/|W|`$ larger than $`\beta _{\mathrm{dyn}}`$). As shown in the Figure, lighter secondaries can be accommodated on wider orbits where the timescale of orbit decay for emission of gravitational waves $`\tau _{\mathrm{GW}}=5c^5r_{\mathrm{orb}}^4/(256G^3M_{\mathrm{tot}}Mm_{\mathrm{in}})`$ is longer ($`M_{\mathrm{tot}}=M+m_{\mathrm{in}}`$ is the total mass of the collapsing core).

Along the relation $`r_{\mathrm{orb}}m_{\mathrm{in}}`$ for $`J_{\mathrm{orb}}=J_{\beta ,\mathrm{dyn}}`$, the bold dots show values of $`m_{\mathrm{mmc}}`$ at different cooling times $`\tau _{\mathrm{cool}}`$ according to the cooling models of Strobel et al. (1999). The dot-dashed lines in the sample plane correspond to loci of fixed values of the gravitational radiation decay timescale for the orbit, $`\tau _{\mathrm{GW}}`$ (from the top, 1, 30, 100 and 300 seconds, as labeled in the Figure).

Notice that gravitational radiation timescales are generally longer than cooling timescales, with the only exception occurring when $`m_{\mathrm{mmc}}0.7M_{}.`$ Coalescence of the two stars is likely to occur in this case with some mass loss that can only be estimated via hydrodynamical simulations. In the opposite case when $`m_{\mathrm{in}}`$ and $`r_{\mathrm{orb}}`$ are such that $`\tau _{\mathrm{GW}}`$ is longer than the cooling time $`\tau _{\mathrm{cool}},`$ the binary survives for a while and below we study its evolution.

### 4.2. Binary evolution, mass transfer, and maximum kicks

In the case of late formation (case C), a phase of stable mass exchange may set in, driven by gravitational radiation emission, which terminates when the secondary reaches a critical mass, $`m_{\mathrm{tid}}`$ (Clark & Eardley 1977; Blinnikov et al. 1984; Jaranowski & Krolak 1992; Bildsten & Cutler 1992; Imshennik & Popov 1998). Subsequently, the low mass star (which has a core-envelope structure) loses its extended halo evolving along a sequence of hydrostatic equilibria while the orbit widens (Blinnikov et al 1984); the star reaches the minimum stable neutron star mass, whereupon it explodes (Page 1982; Colpi, Shapiro & Teukolsky 1989, 1991, 1993; Blinnikov et al. 1990; Sumiyoshi et al 1998). In the case of early formation of the binary (B) we show below that mass transfer is likely to be unstable, and we argue that because of this it may just drive the star to the point of exploding dynamically.

Consider a binary in which Roche lobe overflow is in effect; then, the radius of the low mass star $`R(m)`$ overfills the Roche lobe and this determines the typical orbital separation at the onset of mass transfer

$$r=2.2R(m)(M+m)^{1/3}/m^{1/3}.$$
(13)
The dotted lines in Figure 2 give the tidal radius as a function of $`m`$ for the two different mass-radius relations describing a cold and warm neutron star, respectively. A ring/disk of material forms around the primary star at the moment of Roche filling, and accretion begins to the primary. As the disk drains angular momentum from the orbit there may not be sufficient time to return it to the orbit through disk-donor tidal torques (see for details Bildsten & Cutler 1992; Lubow & Shu 1975). If $`1f`$ is the fraction of orbital angular momentum stored in the disk, the orbital separation would vary with time (following Bildsten & Cutler 1992) as

$$\frac{\dot{r}}{r}=\frac{2\dot{m}(fM^2m^2)}{Mm(M+m)}\frac{64G^3mM(M+m)}{5r^4}.$$
(14)
Fitting the values of $`f(m/M)`$ computed by Hut & Paczynski (1984) for the Roche problem, Bildsten & Cutler (1992) found

$$f(5/3)(m/M)^{1/3}(3/2)(m/M)^{2/3}.$$
(15)
During Roche spillover $`m/R(m)^3(M+m)/r^3`$ and the rate of change of $`r`$ is

$$\frac{\dot{r}}{r}=\left(\frac{d\mathrm{ln}R(m)}{d\mathrm{ln}m}\frac{1}{3}\right)\frac{\dot{m}}{m}$$
(16)
where $`\dot{m}`$ is the mass transfer rate. Using this result in equation (14) for the orbital evolution, we have

$$\left[\frac{1}{3}\frac{d\mathrm{ln}R(m)}{d\mathrm{ln}m}\frac{2(fM^2m^2)}{M(M+m)}\right]\frac{\dot{m}}{m}=\frac{64G^3mM(M+m)}{5r^4}.$$
(17)
Mass transfer is stable only if

$$\frac{1}{3}\frac{d\mathrm{ln}R(m)}{d\mathrm{ln}m}\frac{2(fM^2m^2)}{M(M+m)}<0$$
(18)
given the fact that $`\dot{m}/m<0.`$

In case (C), if we adopt the mass-radius relationship as in Jaranowski & Krolak (1992) valid for a cold star (and mass below $`1M_{}`$), $`R(m)=R_0m^b/(mm_0)^a`$ where $`R_0=7.5`$ km, $`m_0=0.09M_{},`$ $`b=0.79,`$ and $`a=0.83,`$ we find that the condition for stable mass transfer is

$$T\frac{1}{3}+\frac{a+(ab)(m/m_0)}{m/m_01}\frac{2(fM^2m^2)}{M(M+m)}<0.$$
(19)
In general there is only a finite interval for the mass $`m`$ for stability whose boundaries are determined by the two root of the equation $`T=0`$; $`m_{\mathrm{tid}}^{}`$ (the lower bound) and $`m_{\mathrm{tid}}^+`$ (the upper bound). When $`m_{\mathrm{tid}}^{}<m<m_{\mathrm{tid}}^+`$ mass transfer is stable. Outside this interval stability is lost. There may also be no roots, in which case the fate of the light star depends on a number of details that we will discuss below.

Adopting a constant value for $`f,`$ and setting $`\beta \gamma ,`$ to simplify analysis, we note that one regime where the necessary inequality fails is at low masses, near $`mm_0,`$

$$\frac{m_{\mathrm{tid}}^{}}{m_0}1+\frac{\beta }{2f1/3}$$
(20)
(For $`f=0.4`$ and $`\beta =\gamma =0.8`$, $`m_{\mathrm{tid}}^{}/m_0=2.7`$ while if $`f=1`$, $`m_{\mathrm{tid}}^{}/m_01.5.`$) The other regime of instability is at masses well above $`m_0,`$ in which case $`d\mathrm{ln}R(m)/d\mathrm{ln}m`$ is relatively small, and $`m=m_{\mathrm{tid}}^+`$ with

$$\frac{m_{\mathrm{tid}}^+}{M}\frac{1}{2}\sqrt{4f\frac{23}{36}}\frac{1}{12}.$$
(21)
(For $`f=0.4`$, $`m_{\mathrm{tid}}^+/M0.4`$ while for $`f=1,`$ $`m_{\mathrm{tid}}^+/M0.8.`$) Thus, if there is a disk that stores orbital angular momentum (i.e., $`f<1`$), the range of masses for which there can be stable, Roche filling transfer shrinks considerably relative to the case $`f=1`$ explored in Blinnikov et al. (1984) and Imshennik & Popov (1998).

We have computed the roots of equation (19) numerically for different $`R(m)`$. For the Jaranowski & Krolak (1992) parametrization of $`R(m)`$ for cold neutron stars (C) we find the following intervals of $`(m_{\mathrm{tid}}^{},m_{\mathrm{tid}}^+)`$ in units of $`M_{}`$: (0.14,1.08) or (0.14,1.33) for $`f=1`$ and $`M=1.4M_{}`$ or $`1.7M_{}`$; for $`f`$ given by equation (15), there are no roots for $`M=1.4M_{}`$, and the roots are (0.29,0.50) for $`M=1.7M_{}`$; the stability interval $`(m_{\mathrm{tid}}^{},m_{\mathrm{tid}}^+)`$ depends sensitively on the logarithmic derivative of $`R(m).`$ Note that the root $`m_{\mathrm{tid}}^{}`$ exists because of the special dependence of $`R(m)`$ on the parameter $`m_0.`$ Numerical hydrodynamical models for the description of these semi-detached tight binaries are necessary to address the problem of mass transfer, and preliminary results points toward the existence of stable phases (Davies, private communication). Below, we will adopt the values of $`m_{\mathrm{tid}}^{}`$ and $`m_{\mathrm{tid}}^+`$ for a 1.7$`M_{}`$ neutron star.

What is the fate of the binary then ? In case (C) if the two roots $`m_{\mathrm{tid}}^{}`$ and $`m_{\mathrm{tid}}^+`$ exist, and if the initial mass of the low mass star is $`m_{\mathrm{tid}}^{}<m_{\mathrm{in}}<m_{\mathrm{tid}}^+(M)`$, the binary may first go through a stage of evolution via gravitational radiation without mass exchange, then a phase of evolution with mass exchange at a rate determined by the rate of angular momentum loss due to gravitational radiation, and finally become unstable when $`m_{\mathrm{tid}}^{}`$ is reached (Blinnikov et al. 1984; Imshennik & Popov 1998). If, on the other hand, $`m_{\mathrm{mmc}}m_{\mathrm{in}}m_{\mathrm{tid}}^{}`$, then the binary tightens as a consequence of gravitational radiation without mass loss, until $`m`$ fills its Roche lobe, whereupon it becomes unstable. If we define $`m_{\mathrm{stab}}`$ to be the stellar mass at the onset of instability, then $`m_{\mathrm{stab}}=m_{\mathrm{tid}}^{}`$ if $`m_{\mathrm{in}}>m_{\mathrm{tid}}^{}`$, and $`m_{\mathrm{stab}}=m_{\mathrm{in}}`$ if $`m_{\mathrm{mmc}}m_{\mathrm{in}}m_{\mathrm{tid}}^{}`$.

What happens once $`m=m_{\mathrm{tid}}^{}`$ sets in has been described by Blinnikov et al. (1984). As already mentioned, they argue that although the equilibrium radius of a low mass neutron star may become larger than the Roche radius, the vast majority of its mass may still be enclosed within it (99$`\%`$ is contained inside the inner 38 km or so). The stripping of the stellar mass envelope occurs slowly enough that the companion evolves through a series of nearly equilibrium states until $`m_{\mathrm{mmc}}`$ is attained, while the orbit widens.

More accurately, equation (14), in the absence of gravitational wave losses, leads to a change in orbital radius $`r`$ as a function of mass loss equal to

$$\frac{d\mathrm{ln}r}{d\mathrm{ln}m}=2f+\frac{2fm}{M+m}+\frac{2m^2}{M(M+m)}$$
(22)
which can be integrated to yield when the mass is $`m`$

$$\frac{r}{r_i}\left(\frac{m_{\mathrm{in}}}{m}\right)^{2f}\mathrm{exp}\left[\frac{2(1f)(m_{\mathrm{in}}m)}{M+m}\right],$$
(23)
neglecting changes in $`M`$. In this approximation the orbit widens by a factor $`(m_{\mathrm{stab}}/m_{\mathrm{mmc}})^{2f}`$ as the mass decreases from the value $`m_{\mathrm{stab}}`$ toward $`m_{\mathrm{mmc}}.`$

The orbital radius at the point of explosion is somewhat larger than the radius of Roche spillover, and is accordingly

$$r_{\mathrm{mmc}}\frac{2.2R(m_{\mathrm{stab}})M^{1/3}m_{\mathrm{stab}}^{(2f1/3)}}{m_{\mathrm{mmc}}^2}.$$
(24)
The orbital speed of the low mass neutron star when it explodes is

$$V\sqrt{\frac{GM}{r_{\mathrm{mmc}}}}\frac{0.68G^{1/2}M^{1/3}m_{\mathrm{mmc}}}{m_{\mathrm{stab}}^{(f1/6)}[R(m_{\mathrm{stab}})]^{1/2}},$$
(25)
implying a maximum kick speed

$$V_{\mathrm{kick},\mathrm{max}}\frac{m_{\mathrm{mmc}}}{(m_{\mathrm{mmc}}+M)}V\frac{0.68G^{1/2}m_{\mathrm{mmc}}^2}{m_{\mathrm{stab}}^{(f1/6)}[R(m_{\mathrm{stab}})]^{1/2}}M^{2/3}$$
(26)
that scales as $`M^{2/3},`$ due to the dependence of $`V`$ on $`r_{\mathrm{mmc}}`$. This is a direct consequence of the hypothesis of Roche lobe mass transfer; the system has lost memory of the initial orbital separation. Also note that in $`V_{\mathrm{kick},\mathrm{max}}`$ we have $`m_{\mathrm{mmc}}`$ at numerator, even when $`m_{\mathrm{stab}}=m_{\mathrm{tid}}^{}.`$ When $`m_{\mathrm{stab}}=m_{\mathrm{tid}}^{}0.29M_{},`$ (or when $`m_{\mathrm{stab}}=m_{\mathrm{mmc}}=0.0925M_{}`$) the maximum kick is $`V_{\mathrm{kick},\mathrm{max}}936\mathrm{km}\mathrm{s}^1`$ ($`800\mathrm{km}\mathrm{s}^1`$), and $`f0.4`$.

Concerning case (B), warmer proto-neutron stars seem to satisfy a smoother mass-radius relation. A fit to Strobel et al. (1999) data gives $`R(m)m^{6/5}`$ which would lead to a stability interval $`(m_{\mathrm{mmc}},m_{\mathrm{tid}}^+)`$ only for $`f`$ close to unity (note that there is only a single root in this case). Since disk-donor torquing may act on a timescale longer than the lifetime of the binary, it is likely that mass transfer is always unstable under these circumstances, and as soon as the light star overfills its Roche lobe at $`m_{\mathrm{in}}`$, unstable mass transfer begins and there might be time for a phase of common envelope evolution. In this context, we wish to argue that the star is driven almost immediately toward explosion. The cooling time $`\tau _{\mathrm{cool}}50`$ sec (case B; see Figure 2) exceeds the typical time for mass transfer $`\tau _{\mathrm{mass}}m_{\mathrm{mmc}}/\overline{\rho }c_{\mathrm{sound}}R^20.02`$ sec ($`\overline{\rho }10^{11}\mathrm{g}\mathrm{cm}^3`$ denotes the density, close to neutron drip, of the mass losing star in its envelope, $`R50`$ km the size of the Roche lobe and $`c_{\mathrm{sound}}`$ is the sound speed in units of $`10^{10}\mathrm{km}\mathrm{s}^1`$). In this case, the mass losing star encounters the instability point $`m_{\mathrm{mmc}}(T)`$ before being stabilized by cooling. We speculate that the core of the light star explodes on its own dynamical time ($`0.001`$ sec) before coalescence of the two stars in completed over a orbital period of $`P_{\mathrm{orb}}=2\pi (r_{\mathrm{mmc}}^3/GM)^{3/2}0.03`$ sec. Common evolution is known to be accompanied by ejection of part of the envelope enshrouding the system, and this mass loss may provide an additional thrust to the merged object. Only hydrodynamical simulations can give credit to this scenario, and we expect a continuous transition in the physics when moving from C (cold late type mode) to A (hot hydrodynamical mode) passing through B(warm early time mode). If we were to compute $`V_{\mathrm{kick},\mathrm{max}}`$ by assuming that the warm star explodes right at the time it fills the Roche lobe, we would obtain nominal speeds close to $`10,000\mathrm{km}\mathrm{s}^1.`$ The actual value of the kick velocity is computed below combing gravitational bending with orbital phase averaging of the momentum impulse.

### 4.3. Velocity phase-averaging and kick speeds

Our estimate of $`V_{\mathrm{kick}}`$ for cases (C) and (B) using the expression for $`\eta `$ would be not complete without considering that explosion does happens non instantaneously. The star explodes on a timescale comparable to the dynamical timescale, and this can be close to the time of revolution in the binary. In case (B) in particular this effect might be severe.

The final kick speed will be diminished if the explosion is not instantaneous to a good approximation. The orbital frequency when the companion explodes is

$$\omega _{\mathrm{mmc}}=\left(\frac{GM}{r_{\mathrm{mmc}}^3}\right)^{1/2}\frac{0.314G^{1/2}m_{\mathrm{mmc}}^3}{m_{\mathrm{stab}}^{(3f1/2)}[R(m_{\mathrm{stab}})]^{3/2}},$$
(27)
which ranges from $`\omega _{\mathrm{mmc}}17`$ s-1 for $`m_{\mathrm{stab}}=m_{\mathrm{mmc}}`$ to $`\omega _{\mathrm{mmc}}360`$ s-1 for $`m_{\mathrm{stab}}=m_{\mathrm{tid}}^{}0.29M_{}`$ and $`M=1.7M_{}`$ (case C), and $`185`$ s-1 for case B with $`m_{\mathrm{mmc}}=0.3M_{}`$. At the time of explosion, the light star overfills its Roche radius $`R_{R,\mathrm{exp}}R(m_{\mathrm{stab}}).`$ Matter ejected from the unstable star at a velocity $`w_o`$ crosses it in a time $`R_{R,\mathrm{exp}}/w_o`$ and we can take $`R_{R,\mathrm{exp}}/w_o`$ as an estimate of the duration of the explosion, $`\tau _{\mathrm{exp}}.`$

The explosion is almost instantaneous as long as the dimensionless combination

$$\frac{\omega _{\mathrm{mmc}}R_{R,\mathrm{exp}}}{w_o}1.$$
(28)
For $`m_{\mathrm{stab}}0.0952M_{}`$ the dimensionless ratio is $`2,800\mathrm{km}\mathrm{s}^1/w_o0.1`$ for $`m_{\mathrm{stab}}m_{\mathrm{mmc}}`$ (case C), and $`9,000/w_o\mathrm{km}\mathrm{s}^10.3`$ for $`m_{\mathrm{stab}}=m_{\mathrm{mmc}}=0.3`$ in case B, taking $`w_o30,000\mathrm{km}\mathrm{s}^1`$ as an example.

To be more quantitative, suppose that the star loses mass at a rate $`\dot{m}(t)=m\mathrm{\Gamma }(t)`$ once it becomes unstable, with

$$_0^{\mathrm{}}𝑑t\mathrm{\Gamma }(t)=1.$$
(29)
Consider circular orbits only, and assume that although mass loss may extend over many orbital periods, the orbit of the exploding star remains unaltered during the mass loss. (The latter approximation should be valid if mass loss is rapid enough, to zeroth order in the decreasing secondary mass.) Adopt a fixed coordinate system defined at $`t=0`$, when mass loss begins, so that $`\widehat{𝒙}=𝒓(0)/r`$ and $`\widehat{𝒚}=𝑽(0)/V`$. The momentum per unit ejected mass is now

$`𝒖`$ $`={\displaystyle _0^{\mathrm{}}}𝑑t\mathrm{\Gamma }(t)\left[u_x(\widehat{𝒙}\mathrm{cos}\omega t+\widehat{𝒚}\mathrm{sin}\omega t)\right]+`$
$`{\displaystyle _0^{\mathrm{}}}𝑑t\mathrm{\Gamma }(t)\left[+u_y(\widehat{𝒙}\mathrm{sin}\omega t+\widehat{𝒚}\mathrm{cos}\omega t)\right]`$
$`u_+\widehat{𝒆}_{}\mathrm{\Gamma }_\omega +u_{}\widehat{𝒆}_+\mathrm{\Gamma }_\omega ^{},`$
where $`u_\pm =(u_x\pm iu_y)/\sqrt{2}`$, $`\widehat{𝒆}_\pm =(\widehat{𝒙}\pm i\widehat{𝒚})/\sqrt{2}`$, and

$$\mathrm{\Gamma }_\omega _0^{\mathrm{}}𝑑t\mathrm{\Gamma }(t)e^{i\omega t}.$$
(31)
Here, $`(u_x,u_y)`$ are the same as were computed in § 3 for instantaneous mass loss (at $`t=0`$). When the finite duration of the mass loss is accounted for, the efficiency of the recoil is diminished from $`\eta `$ (as computed in § 3 for instantaneous mass loss) to $`\eta \mathrm{\Gamma }_\omega `$. For example, if mass is lost at a constant rate, then $`\mathrm{\Gamma }(t)=\tau _{\mathrm{exp}}^1\mathrm{\Theta }(\tau _{\mathrm{exp}}t)`$, and we find

$$\mathrm{\Gamma }_\omega =\frac{|\mathrm{sin}(\omega _{\mathrm{mmc}}\tau _{\mathrm{exp}}/2)|}{\omega _{\mathrm{mmc}}\tau _{\mathrm{exp}}/2};$$
(32)
if instead $`\mathrm{\Gamma }(t)=\tau _{\mathrm{exp}}^1\mathrm{exp}(t/\tau _{\mathrm{exp}})`$ then

$$\mathrm{\Gamma }_\omega =\frac{1}{\sqrt{1+(\omega _{\mathrm{mmc}}\tau _{\mathrm{exp}})^2}}.$$
(33)
The correction due to phase averaging, $`\mathrm{\Gamma }_\omega `$, is only modest (and somewhat model dependent) as long as $`\omega _{\mathrm{mmc}}\tau _{\mathrm{exp}}1`$, but for long explosions, the efficiency is diminished by a factor of order $`(\omega _{\mathrm{mmc}}\tau _{\mathrm{exp}})^1`$ in general.

Figure 3 shows the kick velocity $`V_{\mathrm{kick}}`$ of the neutron star that remains (of mass 1.7$`M_{}`$; but for 1.4 $`M_{}`$ scale the velocity with $`M^{2/3}`$) as a function of the speed of the ejecta $`w_o`$ for the two cases (C) and (B). The solid curves refer to scenario (C), when $`f`$ is given by the Roche value (eq. ) and $`m_{\mathrm{stab}}`$ is equal to $`m_{\mathrm{tid}}^{}=0.29M_{}`$ (lower curve before crossover occurring at $`w_o40,000\mathrm{km}\mathrm{s}^1`$)) and $`m_{\mathrm{mmc}}`$ (upper curve before crossover). The dotted cure is for $`m_{\mathrm{stab}}=m_{\mathrm{tid}}^{}0.14M_{}.`$ The dash-dotted line, for case (B), is computed starting the explosion during common envelope, when $`m_{\mathrm{mmc}}0.3M_{}`$ and when the orbit separation is $`P_{\mathrm{orb}}/\tau _{\mathrm{cross}}=P_{\mathrm{orb}}(V/r)`$ times less than the separation at the moment of filling the Roche lobe (eq. ), to mimic orbit decay during common envelops. We have included the phase-average orbit correction as given in equation (33).

Very high kicks come from such complex hydrodynamical process that alternative models have not been able to produce.

## 5. Discussion

In this paper, we have elaborated on suggestions that substantial neutron star recoil can result from the explosion of a low mass neutron star formed in the aftermath of rotating core collapse (Imshennik & Popov 1998). A noteworthy feature of this process is that the final kick speed is determined by nuclear physics. The details about the formation mechanism and hydro-dynamical effects produce in reality a widespread range in velocities, as the one we see.

We have extended earlier work by estimating the effect of the finite speed of the material ejected in the explosion, including both gravitational deflection, phase averaging, and indirectly taking into account cooling effects. We argue that a range of neutron star recoil speeds could arise from disruption scenarios with different companion masses. The recoil speed resulting from the dissolution of an evanescent binary system can be around 1000 $`\mathrm{km}\mathrm{s}^1`$, explaining the larger values deduced observationally. It is remarkable that nuclear physics implies a value in the observational range at all, once such a binary system is presumed to form, a concordance that lends some support to the idea.

The kicks that result from this mechanism are confined to the orbital plane of the evanescent neutron star binary. We expect the spin angular momentum vector of the remnant neutron star to be nearly, if not perfectly, aligned with the orbital angular momentum of the binary, and in turn aligned with the angular momentum vector of the collapsing iron core. Near alignment would be consistent with the requirements imposed by the observation of geodetic precession of B1913+16, where the kick is constrained to lie very nearly in the plane of progenitor binary, which was most likely perpendicular to the spins of the spun-up neutron star (i.e. B1913+16) and its pre-explosion companion star (Wex, Kalogera, & Kramer 2000). In contrast, though, X-ray observations of the Vela pulsar have revealed a jet parallel to its proper motion (Pavlov et al. 2000; Helfand, Gotthelf, & Halpern 2001), and it has been argued that the proper motions of both Vela and the Crab pulsar are closely aligned with their spin axes (Lai, Chernoff, & Cordes 2001). The kick mechanism studied here would not be able to account for parallel spin and velocity. However, we note that the proper motions of both Vela and the Crab correspond to transverse speeds of $`70141\mathrm{km}\mathrm{s}^1`$ and $`171\mathrm{km}\mathrm{s}^1`$ respectively, using reasonable estimates of the distances to the pulsars (Lai, Chernoff, & Cordes 2001). For the alignment to be real, the space velocities of these two systems must lie in the plane of the sky. The inferred speeds of these two pulsars are then considerably smaller than the characteristic speed arising from explosion of a low mass, tidally disrupted companion. We propose that the formation of Vela and the Crab did not produce an evanescent binary, and therefore some other mechanism was responsible for their spin-aligned kicks (such as those explored by Lai, Chernoff, & Cordes 2001). We suggest that the larger kick component is due to formation and disruption of an evanescent binary, and is perpendicular to the spin axis, and that the smaller kick component is associated with other less vigorous kicks that tend to align with the rotation axis, perhaps because of phase averaging (e.g. Spruit & Phinney 1998; Lai, Chernoff, & Cordes 2001). A superposition of these two classes of kicks would also be consistent with the requirement of nearly but not precisely spin-perpendicular kicks to account for the observation of geodetic precession in B1913+16 (Wex, Kalogera, & Kramer 2000). If this idea is correct, then one expects that lower velocity neutron stars should have their space velocities predominantly along their spin axes, and high velocity neutron stars should have space velocities predominantly perpendicular to their spins, giving origin to two independent (low and high velocity) distributions. A contamination of low velocity stars (belonging to the low velocity tail of the high velocity distribution) would come from those explosions in the evanescent binary where light bending and rotational averaging have been more important. This could give rise to a bimodal distribution of velocities as is inferred from the observations (Arzoumanian, Chernoff, & Cordes 2001).

This scenario predicts a clear signature in the neutrino emission, in the aftermath of core collapse: Two neutrino bursts occurring several seconds or minutes after core bounce, i.e., after the main neutrino burst, should signal (i) Kelvin contraction of the condensation gathered by the instability and (ii) explosion (on the dynamical time) of the light neutron star. Neutrino emission between the two bursts should also come form material that accretes onto the neutron star that remains. Heavy $`r`$process elements, debris of the explosion of the light star, should also be ejected and found deep in the supernova expanding shells.

Finally, we note that the same phenomenon could occur if the formation of a light neutron star accompanies rotating core collapse to a black hole. In that case, we find that the recoil velocity of the black hole will be proportional to $`M_{\mathrm{BH}}^{2/3}`$, which is not the simple $`M_{\mathrm{BH}}^1`$ scaling one would expect for a kick mechanism that ejects the same amount of momentum irrespective of whether a supernova leaves behind a neutron star or a black hole 222During revision of this manuscript we learned of a paper by Davies et al. (astro-ph/0204358) in which a connection bewteen gamma-ray bursts and pulsar kicks is made, where kicks result from recoil by a short-lived binary.. The higher kick implied by this scaling could be compensated by stronger relativistic effects of gravitational bending and unstable mass transfer that can lead to a much lower escape probability of the exploding debris.

The authors thank Melvyn Davies, John Miller and Andrea Possenti for useful discussions, and the Referee for his/her critical comments. I.W. acknowledges support from NASA Grant NAG 5-8356, and M.C. from ASI and the Ministry of the University and Research (MIUR) under Grant MM02C71842-001.


format: text

text:
# On the Aizenman exponent in critical percolation

## Abstract

The probabilities of clusters spanning a hypercube of dimensions two to seven along one axis of a percolation system under criticality were investigated numerically. We used a modified Hoshen–Kopelman algorithm combined with Grassberger’s “go with the winner” strategy for the site percolation. We carried out a finite-size analysis of the data and found that the probabilities confirm the Aizenman’s proposal of the multiplicity exponent for dimensions three to five. A crossover to the mean-field behavior around the upper critical dimension is also discussed.

Percolation occurs in many natural processes, from electrical conduction in disordered matter to oil extraction from field. In the latter, the coefficient of oil extraction from oil sands (the ratio of the actually extracted to the estimated oil) can be as much as 0.7 for light oil and as low as 0.05 for viscous heavy oil. An increase in this coefficient by any new point requires appreciable investment. Additional knowledge about the percolation model could reduce the amount of additional investment.

A remarkable breakthrough in the theory of critical percolation was established in the last decade thanks to a combination of mathematical proofs, exact solutions, and large-scale numerical simulations. Recently, Aizenman has proposed a new exponent that describs the probability $`P(k,r)`$ of a critical percolation $`d`$-dimensional system with the aspect ratio $`r`$ being spanned by at least $`k`$ clusters Aiz97 ,

$$\mathrm{ln}P(k,r)\alpha _dk^\zeta r,$$
(1)
where $`\alpha _d`$ is a universal coefficient depending only on the universality class, and $`\zeta =d/(d1)`$.

In two dimensions, Aizenman’s proposal (1) was proved mathematically Aiz97 , confirmed numerically SK1 , and derived exactly Cardy98  using conformal field theory and Coulomb gas arguments. This exponent seems to be related to the exponents of two-dimensional copolymers Dup . In three dimensions, proposal (1) was checked numerically in Lev99  and, more recently and more precisely, in Grass02 .

The upper critical dimension of percolation is $`d_c=6`$, which follows from the comparison of the exponents derived on the Cayley tree with those satisfying scaling laws (see, e.g., StAh  and BH ). The fractal dimension $`D_f`$ of percolating critical clusters is equal to 4 above $`d_c`$, and the number of percolating clusters becomes infinite for $`d>d_c`$. This fact would imply that $`\zeta =0`$ at $`d=6`$ if we supposed (rather naively) that Aizenman’s formula applies at the upper critical dimension. Supposing that this is true and taking into account that the values of $`\zeta `$ for $`d=2`$ and $`d=3`$ are, respectively, $`2`$ and $`1.5`$, we can place all three points on the straight line $`\zeta =(6d)/2`$, as depicted in Fig. 1. We can then estimate the respective values of $`\zeta `$ for $`d=4`$ and $`d=5`$ to be $`\zeta =1`$ and $`\zeta =0.5`$; these values are far from those predicted by Aizenman’s formula giving $`4/3`$ and $`5/4`$, respectively. In contrast, based on simulations, Sen Sen97  claims that $`\zeta =2`$ for all dimensions from two to five.

The main purpose of our simulations is to estimate the exponents for the dimensions from two to six with an accuracy sufficient for distinguishing between the values predicted for $`d=4`$ and $`d=5`$ by the Aizenman’s formula and a naive application of cluster fractal-dimension arguments and by the straight-line fit, as discussed above.

In the rest of the paper, we briefly summarize the highlights of our study, then present some details of our research, and finally discuss the results for the Aizenman exponent and the physics of a crossover from the Aizenman picture to the mean-field picture.

Our main results can be summarized as follows:

1. Modified combination of the Hoshen–Kopelman algorithm and Grassberger’s strategy. We use the Hoshen–Kopelman (HK) algorithm HK  to generate clusters and Grassberger’s “go with the winner” strategy Grass02  to track spanning clusters. We add a new tag array in the HK algorithm, which allows the reduction of the tag memory order from $`L^d`$ to $`L^{d1}`$, where $`L`$ is the linear size of the hypercubic lattice. As a result, the amount of memory is about two orders less for large values of $`L`$, and the program is about four times faster—the complexity of the algorithm is compensated by the lower memory capacity needed for for swapping to and from the auxiliary array.

2. Efficient realization of combined shift-register random number generators. We use an exclusive-or ($``$) combination $`z_n`$ of two shift registers:

$`x_n`$ $`=`$ $`x_{n9689}x_{n5502},`$
$`y_n`$ $`=`$ $`y_{n4423}y_{n2325},z_n=x_ny_n`$ (2)
(see RNG  and the references therein). We reduce the computational time for generating random numbers by a factor 3.5 through an efficient technical modification: we use the SSE command set that is available on processors of the Intel and AMD series starting from the Intel Pentium III and AMD Athlon XP.

3. Extraction of the exponents for dimensions three to five. We first use finite-size analysis to estimate the logarithm of the probability $`P(k,r)`$ in the limit of infinite lattice size $`L`$. We then fit data as a function of the number of spanning clusters $`k`$ to obtain the Aizenman exponent $`\zeta `$.

4. Confirmation of Aizenman’s proposal. The estimates of the exponent $`\zeta `$ for the dimensions $`d=2,3,4`$, and $`5`$ coincide well with those proposed by Aizenman.

5. Qualitative interpretation of Aizenman’s conjecture. Cardy interpreted Aizenman’s result qualitatively in two dimensions on the basis of assumption that the main mechanism for reducing the number of percolation clusters is that some of them terminate. The same result can be derived for the cluster confluence (or merging) mechanism. This means that, in low dimensions, the percolation clusters consist of a number of closed paths (loops), while, in higher dimensions, clusters are more similar to trees. Indeed, it is well known that the probability of obtaining loop becomes lower for higher dimensions and goes to zero in the limit of infinite dimensions (Cayley tree) StAh ; BH .

6. Crossover to mean-field behavior. We found evidence that the probability of clusters spanning a hypercubic lattice tends to unity in the limit of high dimensions, as it follows from the well-accepted picture. We did not find any dramatic changes in the probabilities around the upper critical dimension $`d_c=6`$, but rather found evidence for a crossover. Therefore, Aizenman’s formula (1) can also apply to dimensions higher (but not too much higher) than the upper critical dimension and describe approximately the probabilities of spanning clusters in large, though finite-size systems.

We follow with the details of the critical percolation, simulations, and data analysis.

Spanning probability.

We can define the probability $`P(k,r;L)`$ that $`k`$ clusters traverse a $`d`$-dimensional hyperrectangle $`[0,L]^{d1}\times [0,Lr]`$ in the $`Lr`$ direction Aiz97 . Provided that the scaling limit exists (this was proved recently by Smirnov for the percolation in plane Smirnov ), the probability $`P(k,r)`$ can be defined as the limit of $`P(k,r;L)`$ as $`L\mathrm{}`$. Aizenman proposed that $`P(k,r)`$ should behave according to (1) in dimensions from three to five. The validity of formula (1) for the percolation in plane was well established in Aiz97 ; Cardy98 ; SK1 .

Numerical results (Lev99  and Grass02 ) for the exponent $`\zeta `$ for critical percolation on cubic lattices seems to confirm Aizenman’s proposal for the value of $`\zeta =1.5`$

Actually, we could consider the probability $`P(k,r)`$ as the probability of obtaining $`k`$ clusters at the distance $`r`$ from the left side of the hyperrectangle if clusters grow to the right. Only two processes can change the number of clusters: cluster merging and cluster terminating.

The differential $`dP`$ of the probability is

$$dPP(k,r)k^{1/(d1)}kdr,$$
(3)
where the right-hand side represents the product of the probability $`P(k,r)`$ and the differential of the total border hyperarea of $`k`$ clusters, each with the hyperarea differential $`k^{1/(d1)}dr`$. This expression follows from the fact that area unit of measure is proportional to the characteristic transverse length of “infinite” clusters. Therefore, the transverse area remains constant as $`k`$ changes, while the longitudinal length increment in these units is $`k^{1/(d1)}dr`$. Integrating (3), we recover probability (1). Thus, $`P(k,r)`$ describes the probability that $`k`$ clusters do not merge together.

The same probability could be obtained by the process of cluster termination, as given by Cardy in plane Cardy98 , which can easily be extended to dimensions $`d>2`$.

This means that the exponent $`\zeta `$ cannot be larger than the one proposed by Aizenman, and $`\zeta =d/(d1)`$ is the upper bound for the exponent.

Algorithms and realizations.

The classical realization of the HK algorithm HK  requires memory for two major structures: an array for keeping a $`(d1)`$-dimensional cluster slice and a tag array. The total memory required by the algorithm is $`L^{d1}+p_crL^d`$, where $`p_c`$ is the site percolation threshold value. Therefore, for large $`rL`$, one of the main advantages of the HK algorithm (i.e., relatively low memory consumption) is negated by the second term. Our modification of the original algorithm allows the memory for the tag array to be reduced to about $`3p_cL^{d1}`$.

Instead of keeping all tags in memory and selecting new tags with increasing tag numbers, we create two arrays, of which one keeps the tag value and the other one keeps the number $`N`$ of the slice where the corresponding tag was last used. When we build a cluster, we update this array with $`N=N_{\mathrm{current}}`$ for the tags used. If $`N<N_{\mathrm{current}}1`$, then this tag is not on the front surface of the sample, and it will never be used again so that we can, therefore, reuse it. We note that cluster size information should be taken into account before reusing the associated tag, if the size information is required.

We use the “go with the winner” strategy Grass02  as follows. If the system has $`k`$ spanning clusters for some aspect ratio $`r=n\delta r`$, it is stored in memory and is grown for $`\delta r`$. If the resulting configuration has $`k`$ spanning clusters, it is stored, and the growth process continues. Otherwise, we return to the previously saved state. Using this procedure, we calculate the probability $`P_i(\delta r)`$ that the system propagates at the distance $`\delta r`$ from the position $`r=(i1)\delta r`$. Finally, we obtain $`P(r=n\delta r)=_{i=1}^nP_i(\delta r)`$. By choosing sufficiently small values of $`\delta r`$, we can achieve rather high probabilities of $`P_i(\delta r)`$ (which can be determined from a few realizations), while the total probability may be very small (down to $`10^{100}`$ in our case).

The random number generator was optimized for the SSE instruction set as follows. Because the length of all four RNG legs is $`\{a|b\}_{\{x|y\}}>4`$, the $`n`$th step of the RNG does not intersect with the ($`n+3`$)th step. Therefore, we can pack four consecutive 32-bit values of {$`x_{n\{a_x|b_x\}}`$} and {$`y_{n\{a_y|b_y\}}`$} into 128-bit XMM registers, process them simultaneously (see Eq. (2)), and thus obtain $`z_n`$, $`z_{n+1}`$, $`z_{n+2}`$, and $`z_{n+3}`$ within one RNG cycle.

Data analysis.

The lattice size was varied from $`L_{min}`$ to $`L_{max}`$ with the step $`\delta L`$. In Table 1, particular values of the simulation parameters are presented together with the interval of the number of clusters $`k`$ depending on the dimension $`d`$. The direct result of the simulations is the probabilities $`P(k,r;L)`$ that exactly $`k`$ clusters connect two opposite surfaces (separated by the distance $`rL`$) of the rectangle with size $`L^{d1}`$ in the “perpendicular” direction in which we apply periodic boundary conditions. We use the values of the site percolation thresholds on hypercubic lattices from pc2 pc6 , as shown in Table 1.

Data analysis consists of three steps. First, we compute the slope $`s(L)`$ of $`\mathrm{ln}P(k,r;L)`$ for a given dimension $`d`$, number of clusters $`k`$, and linear lattice size $`L`$. An example of such a function is given in Fig. 2 for $`\mathrm{ln}P(5,r;16)`$ in the dimension four. We also plot the logarithm of the probability $`P_+(k,r;L)=_{k^{}k}P(k^{},r;L)`$ of the event that, at least, $`k`$ clusters span the (hyper)rectangle at the distance $`rL`$. To calculate $`s(L)`$, we use data only in the interval of the aspect ratio $`r`$ between $`1.5`$ and $`5`$. We note that the probability of five clusters spanning a rectangle with the linear size $`L=16`$ at the distance $`516=80`$ is extremely small $`10^{52}`$.

Second, we compute probabilities in the limit of an infinite system size $`L`$, fitting slopes $`s(k)`$ with the expression (see Fig. 3)

$$s(k;L)=s(k)+\frac{B}{(L+L_0)^t},$$
(4)
where $`B`$, $`t`$, and $`L_0`$ are fitting parameters Ziff92 ; Lev99 ; Ziff02 . The resulting values of the slopes $`s(k)`$ are presented in Table 2. The number of runs used to compute each particular entry in Table 2 varied from $`10^6`$ to several tens for higher dimensions.

We checked the accuracy of our simulations, as well as the validity of the approach in general for site percolation on a square lattice. Table 3 shows a comparison of our results for the slope $`s`$ with the exact values and with earlier simulations, in which the other modification of the HK algorithm, but not the Grassberger strategy, was used. We note that our results coincide well with the exact results and give a higher accuracy for larger values of $`k`$ in comparison with the previous numerical results, despite the smaller computation time used. Our data for $`k=1`$ is less accurate because of the smaller statistics ($`10^6`$ runs, compared to $`10^8`$ samples in Lev99 ). This is a direct demonstration of the effectiveness of the Grassberger strategy for large values of $`k`$.

Finally, we use values in Table 2 to determine the Aizenman exponent $`\zeta `$ by fitting data in each column to

$$s=A(k^2k_0)^{p/2}$$
(5)
in two and three dimensions, as proposed by Grassberger Grass02 , and to

$$s=A(k^pk_0)$$
(6)
in higher dimensions. Here, $`A`$, $`k_0`$, and $`p`$ are fitting parameters. We take only the leading behavior in $`k`$ into account.

Spanning, proliferation, and crossover to mean-field behavior.

The results of the final fit to (5) and (6) are shown in Table 4. The second row for each particular dimension $`d`$ is the fit with the power $`p`$ fixed to the Aizenman exponent value. This is done to check the fit stability. Indeed, the values of $`A`$ and $`k_0`$ coincide within one standard deviation for the dimensions two to five.

The larger deviations of parameters for the dimensions six and seven may be attributed to appearance of cluster proliferation—the number of clusters is known Aiz97  to grow as $`L^{d6}`$ in dimensions $`d>d_c=6`$. We plot the coefficient $`\alpha _d`$ (defined by expression (1)) in Fig. 4 as a function of the dimension $`d`$. The probability of exactly one cluster spanning at the given distance $`r`$ becomes smaller as the dimension increases from two to five and larger for larger dimensions, as can be seen from the first row ($`k=1`$) of Table 2 and from the lower curve dependence in Fig. 4. For any fixed $`d`$, the value of $`\alpha _d`$ approaches some limit for the dimensions two to five and $`k>2`$, which suggests the value of the corrections to the leading behavior in $`k`$ (see Eqs. (5) and (6)).

The fact that the value of $`\zeta `$, which we formally extracted from our data for $`d=6`$, more or less coincides with $`\zeta =d/(d1)=6/5`$, as formally calculated using the Aizenman expression, may be interpreted as an indication that the number of clusters depends logarithmically on the lattice size $`L`$. One can expect that the logarithmic behavior is visible only for somewhat larger values of $`L`$ than we have used so far (see Table 1). With the values of $`L`$ of the order we have used in simulations, we see effectively the same picture as for the lower dimensions—clusters spann according to the Aizenman formula. This means that at small (or moderate) values of $`L`$, the main mechanism is as discussed above: cluster merging and terminating. And only at sufficiently large system sizes we will see cluster proliferation. An indication of that can be seen from the values of $`\alpha _d`$ in the dimension seven in Fig. 4. The probabilities become closer, and this can be attributed to cluster proliferation and treated as a crossover to the mean-field behavior.

Discussion.

The results have shown the validity of the Aizenman’s proposal in the dimensions three to five (results on plane were already proved rigorously) and do not support Parongama Sen claims based on their simulations (Fig. 1). We have found evidence for cluster proliferation for the dimension seven. The analysis can be extended to the number of spanning clusters to distinguish exponential decay with the system size of the number of clusters for the dimension five, logarithmic growth of them for the dimension six, and linear growth for the dimension seven. The same technique can be used to establish numerically such a crossover to the mean-field picture, although a significantly longer computational time, than we used, is needed for this. In fact, the linear growth of the multiplicity of spanning clusters for seven-dimensional critical percolation was confirmed numerically in preprint ACF  posted at arXiv preprint library a few days after our cond-mat/0207605.

###### Acknowledgements.

We acknowledge useful discussions of algorithms with P. Grassberger and R. Ziff. Our special thanks to S. Korshunov and G. Volovik for the discussion of the results. This work supported by grants from Russian Foundation for Basic Research.


text:
# Mass-loss predictions and stellar masses of early-type stars

## 1. Introduction

The stellar winds of early-type stars are thought to be driven by radiation pressure on spectral lines. Predictions of wind properties are usually tested by comparing them to $`\dot{M}`$ & $`v_{\mathrm{}}`$ values derived for a set of about 30 well studied Galactic and Magellanic Cloud O and early-B stars (Puls et al. 1996). Although it has proven challenging enough to match predictions with observations, these stars provide only a limited test. To better constrain wind theory one should confront a wider range of predictions with observations. Meaningful new tests include extending the comparisons of $`\dot{M}`$ and $`v_{\mathrm{}}`$ to: i) extremely luminous Of and WN5h-6h stars; ii) stars at (both sides of) the bistability jump, an abrupt discontinuity in $`v_{\mathrm{}}`$ found to occur at spectral type B1, and iii) Luminous Blue Variables. We have applied our predictions of mass loss to all the above cases with excellent results (de Koter et al. 1997; Vink et al. 1999, 2000, 2002).

As the properties of stellar winds depend also on stellar mass, an alternative test would be to compare masses derived from line-driven wind theory with those of independent methods. This is a relevant issue in view of current problems with masses based on mass-loss rates and terminal wind velocities.

## 2. Stellar masses derived from wind theory

Although some wind properties (of O stars) are almost mass independent, such as the modified wind momentum $`\dot{M}v_{\mathrm{}}\sqrt{R}`$, the mass loss itself behaves as $`\dot{M}M_{}^{1.3}`$ (Vink et al. 2000). So, in principle, the mass may be derived when all relevant stellar properties – including $`\dot{M}`$ & $`v_{\mathrm{}}`$ – are known. So far, attempts to derive $`M`$ from wind theory yielded systematically lower masses (up to tens of percent) than those derived from evolutionairy tracks, but appeared to be in agreement with spectroscopic mass determinations (e.g. Herrero et al. 2001).

Using the latest generation of model atmospheres (see Herrero this proceedings) the spectroscopic masses, however, now appear to be in much better agreement with $`M_{\mathrm{evol}}`$. Also, direct mass determinations using spectroscopic binaries (e.g. Massey et al. 2002) seem to favor evolutionairy masses. This would leave the wind masses as the exception. However, here we present evidence that also our mass-loss predictions favor $`M_{\mathrm{evol}}`$. The figure shows a consistency check between predicted surface escape velocity $`v_{\mathrm{esc}}[M_{}(1\mathrm{\Gamma })]^{1/2}`$ and $`v_{\mathrm{esc}}`$ implied by the stellar parameters. $`\mathrm{\Gamma }`$ is the Eddington factor for electron scattering. The prediction involves the observed mass loss and terminal wind velocity, i.e. $`v_{\mathrm{esc}}(\mathrm{pred})=f(\dot{M},v_{\mathrm{}},L_{},T_{\mathrm{eff}},M_{})`$. The case for the evolutionairy masses is shown using asterisk symbols; “old” spectroscopic masses (in agreement with previous attempts to derive masses from wind theory) have triangle symbols. The mean ratio for $`M_{\mathrm{evol}}`$ (dotted line) is $`0.90\pm 0.16`$; for $`M_{\mathrm{spec}}`$ (dashed line) it is $`0.79\pm 0.12`$.

This shows that within the error our predicted wind properties are in agreement with evolutionairy masses, bringing also the masses from wind theory in accord with those derived from all other mass determination methods.

## References

de Koter A., Heap S.R., Hubeny I. 1997, ApJ477, 792

Herrero A., Puls J., Corral J., Kudritzki R.-P., Villamariz M.R. 2001, A&A 366, 623

Massey P., Penny L.R., Vukovich J. 2002, ApJ 565, 982

Puls J., Kudritzki R.-P., Herrero A. 1996, A&A305, 171

Vink J.S., de Koter A., Lamers H.J.G.L.M. 1999, A&A 350, 181

Vink J.S., de Koter A., Lamers H.J.G.L.M. 2000, A&A 362, 295

Vink J.S., de Koter A. 2002, A&A in press, astroph/0207170


text:
# Real-time determination of free energy and losses in optical absorbing media

## I Summary

We have introduced notions of free energy and losses relevant to the macroscopic behavior of passive, linear dielectric media. These notions are relevant in the same regime in which the macroscopic Maxwell equations are themselves relevant, i.e. in the regime in which that theory specifies all measurable dynamics. In particular, the macroscopic theory introduced is relevant to the regimes in which the macroscopic, energy barrow-return process describes the production of slow and fast light in passive, linear, temporally dispersive media. Further communications will describe the precise evolutions of the ideal medium-field interactions giving rise to maximum energy recovery, i.e. those ”recovery” fields suggested by the variational definition (6). The nature of the analogy of such recovery fields in dissipative media with reversible processes in the thermodynamic setting will be analyzed. Finally, the nonintuitive evolution of the medium-field interaction on such recovery fields in complex media will be exposed, e.g. the failure of monotonicity in the evolution of the internal energy to its minimum value, even on ideal recovery fields, and the identification of such as a measure of the level of macroscopic disorder created by the existence of many, competing degrees of freedom.

###### Acknowledgements.

Conversations with Kurt Oughstun, Joseph Eberly, Peter Milonni, Aephraim Steinberg, Raymond Chiao, and Michael Fleischhauer are gratefully acknowledged. This research was supported in part by the National Science Foundation under Grant No. PHY99-07949.