text:
# Ordering of Random Walks: The Leader and the Laggard

## I Introduction

Consider $`N`$ identical and independent random walkers that are initially located at $`x_1(0)<x_2(0)<\mathrm{}<x_N(0)`$ on an infinite line remark . There are many interesting questions that can be posed about the order of the particles. For example, what is the probability that all walkers maintain their relative positions up to time $`t`$, that is, $`x_1(\tau )<x_2(\tau )<\mathrm{}<x_N(\tau )`$ for all $`0\tau t`$? By mapping this vicious random walk problem onto the diffusion of a single effective particle in $`N`$ dimensions and then exploiting the image method for the diffusion equation, this ordering probability was found Fi84 ; HF84  to decay asymptotically as $`t^{\alpha _N}`$ with $`\alpha _N=N(N1)/4`$. Many additional aspects of this problem have been investigated within the rubrics of vicious random walks Fi84 ; FG88 ; For90 ; G99 ; Baik00 ; KT02 ; CK02  and “friendly” random walks GV02 ; KGV00 .

In this work, we study two related and complementary random walk ordering problems. In the “leader” problem, we ask for the probability $`_N(t)`$ that the initially leftmost particle in a group of $`N`$ particles remains to the left of all the other particles up to time $`t`$ N83 ; Br91 ; Ke92 ; Kr96 ; Re99 ; Gr02 . In the “laggard” problem, we are concerned with the probability $`_N(t)`$ that the initially rightmost particle from a group of $`N`$ particles never attains the lead (becomes leftmost). These two probabilities $`_N(t)`$ and $`_N(t)`$ decay algebraically

$$_N(t)t^{\beta _N},_N(t)t^{\gamma _N},$$
(1)
as $`t\mathrm{}`$ and our basic goal is to determine the exponents $`\beta _N`$ and $`\gamma _N`$.

These ordering problems arise in a variety of contexts. Physical applications include wetting phenomena Fi84 ; HF84  and three-dimensional Lorentzian gravity DL02 . A more probabilistic application is the ballot problem N83 , where one is interested in the probability that the vote for a single candidate remains ahead of all the other candidates throughout the counting; this is just a restatement of the leader problem. Another example is that of the lamb and the lions Kr96 ; Re99 , in which one is interested in the survival of a diffusing lamb in the presence of many diffusing lions. In one dimension, a lamb that was initially in the lead must remain the leader to survive.

For the leader problem, exact results are known for small $`N`$ only: $`\beta _2=1/2`$ and $`\beta _3=3/4`$ N83 ; Av88 ; FG88 ; Re01 , while $`\beta _N\mathrm{ln}(4N)/4`$ for large $`N`$ Ke92 ; Kr96 ; Re99 ; Gr02 . This slow increase arises because adding another particle has little effect on the survival of the leader when $`N`$ is large. For $`N4`$, no exact results are available and one focus of our work is to obtain an accurate estimate of $`\beta _N`$ for the case $`N=4`$. We accomplish this by mapping the reaction onto an equivalent electrostatic potential problem due to a point charge within an appropriately-defined three-dimensional domain. This mapping provides both an appealing way to visualize the reaction process and an accurate estimate of the survival exponent $`\beta _4`$.

For the laggard problem, we employ the same method as in the leader problem to obtain $`\gamma _3=3/8`$. We also estimate the asymptotic behavior of $`_N`$ and find $`\gamma _NN^1\mathrm{ln}N`$ as $`N\mathrm{}`$. As is expected, a laggard in a large population likely remains a laggard. Therefore the probability of remaining a laggard decays very slowly with time for large $`N`$.

In the next section, we review known results about the leader problem in a 3-particle system. In Sec. III, we outline the electrostatic formulation of the leader problem and then apply it to the 4-particle system in Sec. IV. A numerical solution of the pertinent Laplace equation gives $`\beta _4=0.91342(8)`$, a significant improvement over the previously-quoted estimate $`\beta _40.91`$ Br91 . In Sec. V, we turn to the laggard problem and give an asymptotic estimate for the exponent $`\gamma _N`$. Concluding remarks and some open questions are given in Sec. VI.

## II Conventional Approach to the Leader Problem

We begin by mapping the original problem of $`N`$ diffusing particles $`x_1,x_2,\mathrm{},x_N`$ on the line onto a single effective diffusing particle located at $`(x_1,x_2,\mathrm{},x_N)`$ in $`N`$-dimensional space. The particle order constraints on the line translate to bounding hyperplanes within which the effective particle is confined Fi84 ; Re01 . The effective particle is absorbed if it hits one of these boundaries. For the leader problem, the explicit shape of these bounding hyperplanes can be easily worked out for the cases $`N=2`$ and $`N=3`$; we will later extend this analysis to the 4-particle system.

For 2 particles, their separation undergoes simple diffusion and the process terminates when the separation equals zero. Thus the survival probability of the leader decays as $`t^{1/2}`$. To fix notation and ideas for later sections, we now study the 3-particle system. For a leader at $`x_1(t)`$ and particles at $`x_2(t)`$ and $`x_3(t)`$, we view these coordinates as equivalent to the isotropic diffusion of a single effective particle at $`(x_1(t),x_2(t),x_3(t))`$ in three dimensions. Whenever this effective particle crosses the plane $`A_{ij}:x_i=x_j`$, the original walkers at $`x_i`$ and $`x_j`$ in one dimension have reversed their order. There are three such planes $`A_{12},A_{13},A_{23}`$ that divide space into $`6`$ domains, corresponding to the $`3!`$ possible orderings of the three walkers (Fig. 1(a)). These planes all intersect along the $`(1,1,1)`$ axis.

We may simplify this description by projecting onto the plane $`x_1+x_2+x_3=0`$ that contains the origin and is perpendicular to the $`(1,1,1)`$ axis. Now the plane $`A_{12}`$ may be written parametrically as $`(a,a,b)`$ and its intersection with the plane $`x_1+x_2+x_3=0`$ is the line $`(a,a,2a)`$. Likewise, the intersections of $`A_{13}`$ and $`A_{23}`$ with the plane $`x_1+x_2+x_3=0`$ are $`(a,2a,a)`$ and $`(2a,a,a)`$, respectively (Figs. 1(b)).

The survival of the leader corresponds to the effective particle remaining within the adjacent domains 123 and 132 in Fig. 1(b). The background particles at $`x_2`$ and $`x_3`$ are allowed to cross, but the leader at $`x_1`$ always remains to the left of both $`x_2`$ and $`x_3`$. The union of these two domains defines a wedge of opening angle $`120^{}`$. Since the survival probability of a diffusing particle within a wedge of arbitrary opening angle $`\phi `$ and absorbing boundaries decays asymptotically as $`t^{\pi /2\phi }`$ CJ59 , we deduce the known result that the leader survival probability exponent is $`\beta _3=3/4`$.

## III Electrostatic Formulation

For more than 3 particles, the domain for the effective particle is geometrically more complex and the corresponding solution to the diffusion equation does not seem tractable. We therefore recast the survival probability of the effective diffusing particle in terms of the simpler problem of the electrostatic potential of a point charge within the same geometric domain Re01 . Let $`S(t)`$ be the survival probability of a diffusing particle within an infinite wedge-shaped $`d`$-dimensional domain with absorbing boundary conditions. Let $`V(r)`$ be the electrostatic potential due to a point charge within the same domain, with $`V=0`$ on the boundary. Generically, these two quantities have the asymptotic behaviors

$`S(t)t^\beta ,t\mathrm{},`$
$`V(r)r^\mu ,r\mathrm{}.`$ (2)
More relevant for our purposes, these two quantities are simply related by Kr96 ; Re99 ; Re01

$$_0^tS(t)𝑑t^\sqrt{t}V(r)r^{d1}𝑑r.$$

This equivalence arises because the integral of the diffusion equation over all time is just the Laplace equation. Thus the time integral of the survival probability has the same asymptotic behavior as the spatial integral of the electrostatic potential over the portion of the domain that is accessible by a diffusing particle up to time $`t`$. Substituting the respective asymptotic behaviors from Eqs. (III), and noting that the allowed wedge domain for an $`N`$-particle system has dimension $`N1`$, we obtain the fundamental exponent relation

$$\beta =\frac{\mu N+3}{2}.$$
(3)
Thus the large-distance behavior of the electrostatic potential within a specified domain with Dirichlet boundary conditions also gives the long-time survival probability of a diffusing particle within this domain subject to the same absorbing boundary conditions. From this survival probability, we can then determine the original ordering probability.

To illustrate this approach, let us determine the various ordering probabilities of 3 particles on the line in terms of the equivalent electrostatic problem. In fact, it is simpler to work backwards and find the equivalent ordering problem that corresponds to a specific wedge domain, For example, consider the $`60^{}`$ wedge 123 in Fig. 1. If the effective particle remains within this wedge, the initial particle ordering on the line is preserved. This corresponds to the vicious random walk problem in which no particle crossings are allowed. To obtain the asymptotic behavior of the potential of a point charge interior to this wedge, let us assume that $`V(𝐫)r^\mu f(\phi )`$ as $`r\mathrm{}`$. Substituting this ansatz into the 2-dimensional Laplace equation, we obtain the eigenvalue equation $`f^{\prime \prime }(\phi )=\mu ^2f(\phi )`$, subject to $`f(\phi )=0`$ on the wedge boundaries. For the $`60^{}`$ wedge, the solution with the smallest eigenvalue is $`f(\phi )=\mathrm{sin}(3\phi )`$. Thus $`\mu =3`$, leading to the known result $`\beta =(\mu N+3)/2=3/2`$, for 3 vicious walkers.

## IV 4-Particle System

The state of the system may be represented by an effective diffusing particle in 4 dimensions. By projecting onto the 3-dimensional subspace $`x_1+x_2+x_3+x_4=0`$ that is orthogonal to the $`(1,1,1,1)`$ axis, the order domains of the original particles can be reduced to 3 spatial dimensions. In this 3-dimensional subspace, the boundary $`A_{12}:x_1=x_2`$ becomes the plane $`(a,a,b,c)`$, with $`2a+b+c=0`$. Likewise, $`A_{13}`$ may be written parametrically as $`(a,b,a,c)`$, where again $`2a+b+c=0`$. Thus the locus $`L_{123}A_{12}A_{13}:x_1=x_2=x_3`$ is the line $`(a,a,a,3a)`$. This body diagonal joins the nodes 4 and $`\overline{4}`$ in Fig. 2. Along this axis, the original particle coordinates on the line obey the constraint $`x_1=x_2=x_3`$, with $`x_4<x_3`$ on the half-axis closer to node 4 and $`x_4>x_3`$ on the half-axis closer to node $`\overline{4}`$. A similar description applies to the axes $`L_{124}=(a,a,3a,a)`$ between 3 and $`\overline{3}`$, $`L_{134}=(a,3a,a,a)`$ between 2 and $`\overline{2}`$, and $`L_{234}=(3a,a,a,a)`$ between 1 and $`\overline{1}`$.

The locus where $`x_1=x_2`$ and $`x_3=x_4`$ simultaneously, is the line $`L_{12,34}A_{12}A_{34}=(a,a,a,a)`$. Likewise, $`L_{13,24}=(a,a,a,a)`$, and $`L_{14,23}=(a,a,a,a)`$. Viewed in the orthogonal 3-subspace, the $`6`$ planes $`A_{ij}`$ intersect in $`7`$ lines ($`4`$ $`L_{ijk}`$ and $`3`$ $`L_{ij,kl}`$), and divide the subspace into $`24`$ semi-infinite wedges, as shown in Fig. cube(b). Each of these wedges corresponds to one of the $`4!`$ orderings of the walkers in one dimension.

We first illustrate the electrostatic formulation of this system for the vicious random walk problem. Since the initial particle order is preserved, the effective diffusing particle remains within a single wedge $`ijkl`$ in Fig. 2. As outlined in the previous section, the survival probability of this effective particle corresponds to the electrostatic potential of a point charge within this one wedge, with the boundary surfaces held at zero potential. To solve this electrostatic problem, it is convenient to place the point charge at the symmetric location $`(u,v,w)=(0,1,1/2)`$ within the wedge 1234, where the $`(u,v,w)`$ axes are defined in Fig. 2. From the image method, the potential due to a point charge within one wedge is equivalent to the potential of an array of 24 symmetrically-placed point charges consisting of the initial charge and 23 image charges, with positive images at $`(0,1,1/2)`$, $`\pm (0,1,1/2)`$, $`\pm (\pm 1/2,0,1)`$, $`\pm (1,\pm 1/2,0)`$, and negative images at $`\pm (\pm 1/2,1,0)`$, $`\pm (0,\pm 1/2,1)`$, and $`\pm (0,1,\pm 1/2)`$. Using Mathematica the asymptotic behavior of the potential in wedge 1234 (where the original charge is placed) due to this charge array is

$$V(r)=a_1r^7+a_2r^{11}+a_3r^{13}+a_4r^{15}+\mathrm{}.$$

(The coefficients $`a_i`$ depend on the location of the charge and on the orientation of r.) Using the exponent relation (3), the asymptotic survival probability of $`4`$ vicious walkers is given by

$$S(t)=b_1t^3+b_2t^5+b_3t^6+b_4t^7+\mathrm{}.$$

The leading behavior confirms the known result Fi84 ; HF84 , and we obtain the form of the corrections as well. As an amusing aside, notice that an $`r^7`$ dependence for the potential is normally achieved by a 64-pole charge configuration. The high symmetry of the 24 charges in the ordering problem leads to a multipole field normally associated with at least 64 point charges.

We next turn to the leader problem. This system corresponds to the electrostatic problem within the combined domain of the 6 wedges marked 1234, 1243, 1423, 1432, 1342, and 1324, in Fig. 2. The resulting domain is a tetrahedral corner, with its apex at the center of the cube, that is flanked by the rays $`0\overline{2}`$, $`0\overline{3}`$, and $`0\overline{4}`$. Despite the simplicity of this domain, we are unable to solve this electrostatic problem analytically and we have instead studied the problem numerically.

In the allowed region of the cube in Fig. 2, we discretize space and solve the electrostatic potential of a point charge by using successive over-relaxation, with the domain boundaries at zero potential. For simplicity, the charge is chosen to be at the symmetric point $`(1,1,1)`$. For the outer faces of the cube we use two different boundary conditions: (a) absorbing ($`V=0`$), and (b) reflecting, ($`dV/dn=0`$, where $`n`$ is the direction normal to the surface). The true potential – that of the infinite wedge – lies between these two extremes. We also exploit the symmetry about the $`(1,1,1)`$ diagonal and use only the domain marked 1234 in Fig. 2, with absorbing boundary conditions on the plane $`34\overline{3}\overline{4}`$, and reflecting boundary conditions on the planes $`12\overline{1}\overline{2}`$ and $`14\overline{1}\overline{4}`$. As already discussed, the outer cube face, $`1\overline{3}2\overline{4}`$, is taken to be reflecting or absorbing. This space savings allows us to carry out computations for a cube of $`500`$ lattice spacings on a side.

While the power-law decay of the potential sets in quickly, the finite-size effect is pronounced and it is perceptible already at 25 lattice spacings away from the charge. This is the primary limitation on the accuracy of our exponent estimate. Nevertheless, the local exponent $`\mu (r)=d\mathrm{ln}V(r)/d\mathrm{ln}r`$ varies only at the fourth digit (Fig. 3). The approach of the local exponent to the asymptotic limit also suggests that $`V(r)`$ has the form $`V(r)r^\mu +Ar^4`$. Assuming that this is the case, extrapolation of the data in Fig. 3 gives $`\mu =2.82684\pm 0.00016`$, where the error bar is the difference in the extrapolated value of $`\mu (r)`$ from the two different boundary conditions. From the exponent correspondence given in Eq. (3), we thereby obtain, for the lead probability,

$$_4(t)t^{\beta _4}+At^{3/2},\beta _4=0.91342(8),$$
(4)
It is hard to match this numerical accuracy with that from direct simulations of the survival of the leader. We simulated $`10^9`$ realizations of the system; this gives extremely linear data for the time dependence of the leader survival probability on a double logarithmic scale. To estimate the exponent $`\beta _4`$, we computed the local slopes of the survival probability versus time in contiguous time ranges between $`t`$ and $`1.5t`$ when plotted on a double logarithmic scale. These local exponents are plotted against $`1/\mathrm{ln}_{1.5}t`$ (Fig. 4). The results are compatible but much less accurate than Eq. (4).

## V The Laggard Problem

In the laggard problem, we study the probability that the initially rightmost particle at $`x_N`$ has never been the leader during the time interval $`(0,t)`$. The laggard problem can also be recast into the diffusion of a single effective particle within an $`N`$–dimensional wedge-shaped region, with absorbing domain boundaries. This mapping leads to the basic conclusion that every particle that is initially not in the lead exhibits the same asymptotic behavior as the last particle. Indeed, for any particle $`i`$ initially at $`x_i`$, the regions $`x_ix_1,\mathrm{},x_{i1},x_{i+1},\mathrm{},x_N`$ are isomorphic. The initial condition merely fixes the location of the effective particle in this allowed region. Another fundamental observation is that the allowed regions of the effective particle for the leader and the laggard problems are complementary for all $`N`$.

For two particles, the probability that the laggard does not become the leader obviously decays as $`t^{1/2}`$, i.e., $`\gamma _2=1/2`$. The case $`N=3`$ is more interesting but also solvable. The condition that a particle initially at $`x_3`$ never attaining the lead ($`x_3`$ $`x_1,x_2`$) is equivalent to the effective particle remaining with the lighter shaded region in Fig. 1(b). Since the opening angle of the resultant wedge is $`\phi =4\pi /3`$, the corresponding survival probability decays as $`t^{3/8}`$, implying that $`\gamma _3=3/8`$.

We have performed direct numerical simulations of the process to estimate the exponent $`\gamma _N`$ for $`N=3,4,5`$, and 6. Each simulation is based on $`10^6`$ realizations in which $`N1`$ particles are initially at the origin, while the laggard is at $`x=1`$. Each simulation is run until the laggard achieves the lead or $`10^5`$ time steps, whichever comes first. From the survival probability, we estimate $`\gamma _N0.35,0.30,0.26`$, and $`0.23`$ for $`N=3`$ – 6, respectively. Since we know that $`\gamma _3=0.375`$, the discrepancy of 0.025 between the simulation result and the theory is indicative of the magnitude of systematic errors in this straightforward numerical approach.

While it appears difficult to determine the exponent $`\gamma _N`$ analytically for general $`N>3`$, the situation simplifies in the large $`N`$ limit because the position of the leader becomes progressively more deterministic. Indeed, the probability density $`P_N(y,t)`$ that the leader is located at distance $`y`$ from the origin (assuming that all particles are initially at the origin) is Re99

$`P_N(y,t)={\displaystyle \frac{Ne^{y^2/4Dt}}{\sqrt{4\pi Dt}}}\left[1{\displaystyle \frac{1}{2}}\mathrm{erfc}\left({\displaystyle \frac{y}{\sqrt{4Dt}}}\right)\right]^{N1},`$
where $`\mathrm{erfc}(x)=1\mathrm{erf}(x)`$ is the complementary error function. Performing a large $`N`$ analysis, we find that the probability density $`P_N(y,t)`$ approaches a Gaussian

$$P_N(y,t)\frac{1}{\sqrt{2\pi \sigma ^2}}\mathrm{exp}\left\{\frac{(y+y_{})^2}{2\sigma ^2}\right\},$$
(5)
with the mean and the variance of this distribution given by

$$y_{}=z\sqrt{4Dt},\sigma ^2=\frac{Dt}{z^2},$$
(6)
where $`z`$ is determined from the transcendental relation

$$ze^{z^2}=\frac{N}{\sqrt{4\pi }}.$$
(7)
Consequently, the parameter $`z`$ diverges as $`z\sqrt{\mathrm{ln}N}`$ when $`N\mathrm{}`$. The ratio of the dispersion to the mean displacement thus decreases as $`\sigma /y_{}(2\mathrm{ln}N)^1`$, so that the position of the leader indeed becomes more deterministic as $`N\mathrm{}`$.

Therefore in the large $`N`$ limit we can assume that the leader is moving deterministically and its position is given by $`y_{}(t)`$. Then the probability that the laggard never achieves the lead is equivalent to the probability that a diffusing particle initially at the origin does not overtake a receding particle whose position is varying as $`y_{}(t)`$. This corresponds to the solution to the diffusion equation in the expanding region $`x(y_{}(t),\mathrm{})`$ with an absorbing boundary condition at the receding boundary $`x=y_{}`$. When $`y_{}t^{1/2}`$, this diffusion equation can be solved exactly by reducing it to a parabolic cylinder equation in a fixed region. However, in the limit $`t\mathrm{}`$, we can obtain asymptotically correct results much more simply. Because the absorbing boundary recedes from the laggard particle relatively quickly, we solve the problem by assuming that the density $`P(x,t)`$ of the laggard particle approaches a Gaussian with yet unknown amplitude $`_N(t)`$ Kr96 ; Re99 :

$$P(x,t)=\frac{_N(t)}{\sqrt{4\pi Dt}}\mathrm{exp}\left\{\frac{x^2}{4Dt}\right\}.$$
(8)
Although this distribution does not satisfy the absorbing boundary condition, the inconsistency is negligible since the exponential term in Eq. (8) is of order $`N^1`$ at the boundary $`y=y_{}`$.

The probability $`_N(t)`$ is now found self-consistently by equating the “mass” loss to the flux:

$$\frac{d_N}{dt}=D\frac{P}{x}|_{x=y_{}}$$
(9)
Using Eq. (8) to compute the flux we convert Eq. (9) into

$$\frac{d_N}{dt}=\frac{z^2}{N}\frac{_N}{t},$$
(10)
from which the exponent $`\gamma _N`$ is $`z^2/N`$. This is of course valid only in the large $`N`$ limit. Taking this limit in Eq. (7) we obtain

$$\gamma _N=\frac{\mathrm{ln}N}{N}\frac{\mathrm{ln}\mathrm{ln}N}{2N}+\mathrm{}.$$
(11)
Thus as $`N`$ gets large, the probability that the laggard never attains the lead decays extremely slowly with time. This fits with the naive intuition that if the number of particles is large a laggard initially is very likely to remain a laggard. Each additional particle makes it even less likely that the laggard could achieve the lead. Amusingly, this asymptotic exponent predictions is numerically close to the previously-quoted results from direct numerical simulations of the laggard problem.

## VI Concluding Remarks

We investigated two dual random walk ordering problems in one dimension: (i) what is the probability that a particle, that is initially in the lead, remains in the lead and (ii) what is the probability that a particle, that is initially not in the lead, never achieves the lead? These problems are most interesting in one spatial dimension because of the effective correlations between the interacting particles. These correlations are absent in two dimensions and greater, so that an $`N`$-particle system reduces to a 2-particle system Kr96 ; Re99 ; CK02 .

We determined the respective exponents $`\beta _N`$ and $`\gamma _N`$ associated with the lead and laggard probabilities for general $`N`$. Both exponents can be determined by elementary geometric methods for $`N=2`$ and 3 and by asymptotic arguments for $`N\mathrm{}`$. Our new results are the following: (i) a precise estimate for $`\beta _4`$ and (ii) the large-$`N`$ behavior of $`\gamma _N`$.

A simple generalization is to allow each particle $`i`$ to have a distinct diffusion coefficient $`D_i`$. The exponents $`\beta _N`$ and $`\gamma _N`$ will now depend on the diffusion coefficients, except for $`N=2`$, where $`\beta _2`$ and $`\gamma _2`$ always equal 1/2. The case $`N=3`$ is still solvable by introducing rescaled coordinates $`y_i=x_i/\sqrt{D_i}`$ to render the diffusion of the effective particle isotropic, after which the mapping to the wedge can be performed straightforwardly. We thus find

$`\beta _3`$ $`=`$ $`\left\{2{\displaystyle \frac{2}{\pi }}\mathrm{cos}^1{\displaystyle \frac{D_1}{\sqrt{(D_1+D_2)(D_1+D_3)}}}\right\}^1,`$
$`\gamma _3`$ $`=`$ $`\left\{2+{\displaystyle \frac{2}{\pi }}\mathrm{cos}^1{\displaystyle \frac{D_3}{\sqrt{(D_1+D_3)(D_2+D_3)}}}\right\}^1.`$
An amusing special case is the case of a stationary laggard, for which $`\gamma _3=1/3`$.

Finally, it is also worth mentioning a promising development to solve the diffusion equation in the domains defined by the ordering of one-dimensional random walks. This is the recent discovery of deep connections between vicious walkers and random matrix theory G99 ; Baik00 ; KT02 . These allow one to not only re-derive the exponent $`\alpha _N`$ of the original vicious random walk problem, but also lead to many new results. It would be extremely useful if these techniques could be extended to the leader and laggard problems.

###### Acknowledgements.

We are grateful to NSF grants PHY-0140094 (DbA) and DMR9978902 (PLK and SR) for partial financial support of this research.


format: text

text:
# Peculiar, Low Luminosity Type II Supernovae: Low Energy Explosions in Massive Progenitors? Based on observations collected at the European Southern Observatory and Cerro Totolo Inter-American Observatory, Chile (Proposals ESO 63.H-0141 and 66.D-0683).

## 1 Introduction

Recently, a number of supernovae have been discovered that, according to their spectral properties, are classified as Type II but that, at the same time, show remarkably peculiar properties. The first clearly identified example was the exceptionally faint SN 1997D in NGC 1536 (Turatto et al. ; Benetti et al. ). SN 1997D showed both a very faint radioactive tail in the light curve, indicating an ejected 56Ni mass of only a few $`\times 10^3M_{}`$, and a low expansion velocity of $`1000`$ km s-1. Modeling the spectra and light curve Turatto et al.  concluded that a low energy explosion in a 26 $`M_{}`$ progenitor star could successfully fit the early observations. The scenario of the low energy explosion of a high mass progenitor is also consistent with the late time ($`400`$ days) spectral and photometric data (Benetti et al. ). Chugai & Utrobin  have presented an alternative analysis in which the progenitor was a star at the low end of the mass range of core-collapse supernovae (8–12 $`M_{}`$).

The very low luminosity at discovery (10 times less than the peak luminosity reached by SN 1987A) and very small expansion velocity of the ejecta (3-4 times less than that of a normal Type II) of SN 1997D suggest that this explosion event was under-energetic. The mechanism that causes the energy of the explosion to be so low is a challenging and important problem in the physics of core-collapse supernovae. In fact, the exact evolution of the star after neutrino reheating depends on the rate of early infall of stellar material on the collapsed core and on the binding energy of the envelope. If both are large, as in high mass stars ($`M>20M_{}`$), the energy available to accelerate and heat up the ejecta can be greatly reduced and, eventually, may become insufficient to cause a successful explosion. Therefore, after the passage of the shock wave (and possibly the reverse shock formed at the H–He interface) a variable amount of matter may remain gravitationally bound to the collapsed remnant and fall back onto it (Woosley & Weaver ). For this reason it was early suggested that SN 1997D could host a black hole remnant formed during the explosion (Turatto et al. ) and that the luminosity powered by fallback of envelope material onto the central black hole could emerge at about 1000-1200 days after the explosion (Zampieri, Shapiro & Colpi ). Unfortunately, to date it has not been possible to confirm or disprove observationally this prediction.

After SN 1997D, a number of supernovae with similar observational properties have been identified. These objects appear to define a fairly homogeneous group of explosion events (Pastorello et al. ). As for SN 1997D, these supernovae provide a unique opportunity to probe the physics of the explosion and reach a better understanding of both the explosion mechanism and the conditions for the formation of black hole remnants. In this Paper we present the results of a joint analysis of the observations of two peculiar Type II supernovae with very low luminosity and expansion velocity, SN 1997D and SN 1999br, which are representative of the properties of the whole group. In particular, SN 1999br represents the most extreme case of a low-luminosity event to date and appears to follow a recently reported correlation between expansion velocities of the ejecta and bolometric luminosities during the plateau phase (Hamuy & Pinto ).

In Section 2 we describe the basic spectral and photometric data of these supernovae. Section 3 summarizes the semi-analytic method employed to model the light curve and the evolution of the line velocity and continuum temperature of Type II supernovae. Finally, in section 4 the main results are presented and their consequences regarding the nature of the progenitors and the energy of the explosion are briefly discussed.

## 2 Peculiar Type II supernovae

Supernova 1997D, serendipitously discovered on 14 January 1997 during an observation of the parent galaxy NGC 1536 (De Mello & Benetti ), is the first, clearly identified example of a peculiar Type II supernova with a very low luminosity and expansion velocity (Turatto et al. ). It was detected when it was already decaying from the plateau and was at least 2 mag fainter than a typical Type II supernova (Patat et al. ). The decline rate of the last segment of the bolometric light curve is consistent with complete thermalization of the gamma rays from the radioactive decay of 56Co into 56Fe. Assuming for SN 1997D the same deposition as in SN 1987A, the ejected 56Ni mass is $`10^310^2M_{}`$ (Benetti et al. ). The spectra are dominated by a red continuum and P-Cygni profiles of H I, Ba II, Ca II, Na I and Sc II (see Figure 1). The most striking property of these spectra is the very low expansion velocity inferred from the spectral lines. The minima of the absorption lines give an expansion velocity of about 1100-1200 km s-1 for H (see Benetti et al.  for details).

Recently, other supernovae with properties similar to those of SN 1997D have been identified. A comprehensive analysis of the observations of these peculiar supernovae, including the recently discovered SN 2001dc, will be presented elsewhere (Pastorello et al. ). Here we focus on one particularly representative object, SN 1999br in NGC 4900, discovered on 12 April 1999 (King ). Filippenko et al.  pointed out the exceptionally low luminosity of this event. Later, Patat et al.  recognized that SN 1999br shows other properties similar to those of SN 1997D, such as spectra with narrow P-Cygni lines, prominent Ba II features and a low continuum temperature. This supernova had a long plateau lasting at least 110 days with a mean luminosity of only $`\text{ }<10^{41}`$ erg s-1. It was extensively monitored during the first $`110`$ days with the CTIO and ESO telescopes, yielding extremely good photometric and spectroscopic coverages (Hamuy et al. ). The spectrum of SN 1999br taken at a phase of $`100`$ days shows narrow metal lines (Figure 1). As in the case of SN 1997D, the red continuum and the strength of Ba II lines are caused by the low temperature of the ejecta. The U-B bands are strongly affected by line blanketing as the temperature decreases.

As shown in Figure 1, the spectrum of SN 1997D obtained soon after discovery and that of SN 1999br at $`100`$ days are strikingly similar in both the continuum and the line components. Taking the ratio of the two spectra results in a rather flat curve where the main deviations are caused by slightly different line width (i.e. different expansion velocities of the ejecta) and spectral resolution. The expansion velocity of SN 1997D is only $`5\%`$ larger than that of SN 1999br. The close resemblance of these two spectra, the low inferred expansion velocity ($``$1000 km s-1) and the fact that the two supernovae have comparable luminosities strongly suggest that SN 1997D and SN 1999br may be similar events. In the following we will investigate in detail the consequences of this assumption.

SN 1997D was discovered at the end of the plateau stage when the light curve was plummeting. The duration and luminosity of the plateau were inferred by comparison between models and observations and the consequent estimate of the explosion epoch ($`60`$ days before discovery) was uncertain (Turatto et al. ). On the other hand, the date of the explosion of SN 1999br has an uncertainty of only a few days (thanks to a stringent pre-discovery limit). Then, from the similarity of the observational properties (luminosities and spectra) of SN 1997D and SN 1999br, we tentatively assume that the phase of SN 1997D at discovery is $`90100`$ days. This assumption will be checked by means of a detailed comparison of spectral-synthesis models with observations, presently under way.

## 3 The Model

In order to determine the physical properties of a supernova, the observed light curve and spectra are compared to numerical simulations. Usually, it is important to carry out a preliminary investigation in order to obtain an approximate but reliable estimate of the physical conditions of the ejected gas and to establish a framework in which more detailed follow-up calculations can be performed. For these reasons, we have implemented a semi-analytic model in spherical symmetry that can provide a robust estimate of the parameters of the ejecta of Type II supernovae. The novelty of the present approach lies mainly in the fact that the physical properties of the envelope are derived by performing a simultaneous comparison of the observed and simulated light curve, evolution of the line velocity and continuum temperature. The model follows the approach originally introduced by Arnett  and later developed by Arnett & Fu  and Arnett . The treatment of the motion of the recombination front follows in part the work of Popov . The present analysis includes all the relevant energy sources powering the supernova and evolves the envelope in 3 distinct phases that cover the whole evolution from the photospheric up to the late nebular stages. We include also the energy input from recombination that may be particularly important for low energy events like SN 1997D. The model has been tested against numerical radiation-hydrodynamic computations in spherical symmetry under the same assumptions, giving good agreement. The details of the method and its validation will be presented in a companion paper (Zampieri et al. ). Here we briefly summarize the basic ideas.

Full hydrodynamical calculations using realistic pre-supernova models show that the dynamical evolution of the envelope during shock passage is quite complex. The propagation of the shock determines how the explosion energy is distributed in the envelope of the progenitor star. The innermost layers comprised of helium and heavier elements transfer almost all of their kinetic energy and momentum to the Hydrogen envelope. The star mixes but does not homogenize. The actual velocity, density and heavy element distributions of the post-shock material affect the light curve and estimation of the envelope mass in a major way. In the present analysis we do not consider the evolution of the star during this complex phase but rather assume idealized initial conditions that provide an approximate description of the ejected material after shock (and possible reverse shock) passage, as derived from hydrodynamical calculations. The evolution starts at time $`t_{in}`$ (since core collapse) when the envelope is essentially free-coasting and in homologous expansion. In realistic explosion calculations, it takes several hours in order for the envelope to relax to this state. We adopt $`t_{in}5`$ hrs. At this stage, the velocity $`V`$ of each gas shell is approximately constant and proportional to its position $`r`$, $`V(r)=V_0(r/R_0)`$, where $`R_0`$ and $`V_0`$ are the initial radius and velocity of the outermost shell of the envelope. Therefore, the radius of each gas shell increases linearly with time: for the outermost shell, $`R=R_0+V_0(tt_{in})`$. The post-shock, ejected envelope is assumed to have spatially constant density $`\rho `$ and total mass $`M_{env}=4\pi \rho R^3/3`$. Mass conservation gives $`\rho =\rho _0(R_0/R)^3`$, where $`\rho _0`$ is the initial density. In reality, the outer part of the star develops a steep power-law density structure that affects the light curve during the first few days after shock breakout (which is not included). However, this is only $``$ 1% by mass of the star, while most of the stellar material resides in the inner part with roughly constant density. The latter region dominates the light curve after 10–20 days. The expansion velocity of the envelope as a function of interior (Lagrangian) mass $`m(r)=4\pi \rho r^3/3`$ is then $`V(m)=V_0(m/M_{env})^{1/3}`$ (see Figure 2). The initial thermal+kinetic energy of the envelope is $`E=(3/10)M_{env}V_0^2/f_0`$, where $`f_0`$ is the initial fraction of kinetic energy. Elements are assumed to be completely mixed throughout the envelope. In our simple spherically symmetric model, their distribution depends only on $`r`$ (or $`m`$). In particular, Hydrogen, Helium and Oxygen are assumed to be uniformly distributed, whereas 56Ni is more centrally peaked (see Figure 2). Our assumptions about the velocity and elemental distributions provide only an approximate description of the actual post-shock structure of the envelope and should be regarded as a potential source of systematic uncertainty in the present model.

The evolution of the supernova envelope starts at time $`t_{in}`$ and is schematically divided into 3 phases. During the first stage (a few tens of days) the envelope is hot and completely ionized because of the energy deposited by the shock wave. In the following phase (from 30–40 days up to $`100`$ days), because of the decrease in temperature caused by expansion and radiative diffusion, the envelope recombines and can be schematically divided into two regions, below and above the position of the recombination wavefront. In the third stage (after $`100`$ days), the ejecta are completely recombined and transparent to optical photons. The evolution is computed solving the energy balance equation for the envelope gas. The thermal balance is governed by the competition among the energy input from trapped gamma-rays, the $`PdV`$ work and the energy losses through radiative diffusion. Assuming that radiation is in LTE with the gas throughout the envelope, the energy balance equation becomes a second order partial differential equation for the temperature $`T`$. Arnett  has shown that, under specific assumptions on the spatial distribution of 56Ni, the energy equation (with appropriate boundary conditions) can be solved by separation of variables. The function that describes the radial dependence $`\psi (r)`$ is solution of an eigenvalue equation (see e.g. equation \[D.48\] of Arnett ). The final solution has the form

$$T^4(t,r)=T_0^4\left(\frac{R_0}{R}\right)^4\psi (x)\varphi (t),$$
(1)
where $`x=r/R`$, $`T_0`$ is a constant reference temperature, $`\varphi (t)`$ is the function that describes the time dependence \[$`\varphi (0)=1`$\] and, in the limit of zero mean free path, $`\psi (x)=\mathrm{sin}(\pi x)/(\pi x)`$ (Arnett’s “radiative zero” solution). The term $`(R_0/R)^4`$ accounts for adiabatic expansion. At $`t=t_{in}`$, the radial temperature profile is $`T/T_0=\psi (x)^{1/4}`$ (see Figure 2). Once a solution for $`\varphi (t)`$ is calculated, the emitted luminosity is approximately given by the expression for the diffusion luminosity $`L=(4\pi r^2c)/(3\kappa \rho )[aT^4/r]_{r=R}`$, where $`\kappa `$ is the gas opacity.

When the gas starts to recombine at $`t=t_{rec,0}`$, the part of the envelope below the position $`r_i`$ of the wavefront is assumed to be in LTE with radiation. We adopt the approximation that in this region radiative diffusion can effectively readjust the radial temperature distribution to the changing position of the outer boundary at $`r_i=r_i(t)`$ (“slow approximation”, Arnett & Fu ). In this assumption, the spatial dependence of $`T`$ is essentially given by Arnett’s “radiative zero” solution. Therefore, we search for approximate solutions of the energy equation of the form given in equation (1), with $`x`$ replaced by $`x/x_i`$, $`x_i=r_i/R`$ and $`\psi (x/x_i)=\mathrm{sin}(\pi x/x_i)/(\pi x/x_i)`$ (Arnett ). The function that gives the time dependence of the temperature profile, $`\varphi (t)`$, can be approximately computed setting the diffusion luminosity at recombination equal to the luminosity emitted by a blackbody at the effective temperature, $`4\pi r_i^2\sigma T_{eff}^4`$ (Popov ). This gives

$$\varphi (t)=\frac{3}{4}\kappa \rho _0R_0\left(\frac{T_{eff}}{T_0}\right)^4\left(y^2\frac{d\psi }{dy}\right)_{y=1}^1\left(\frac{R}{R_0}\right)^2x_i,$$
(2)
where $`y=x/x_i`$ and $`T_{eff}`$ is the effective temperature, approximately constant and nearly equal to the gas ionization temperature. Using equations (1) (with $`x`$ replaced by $`x/x_i`$) and (2), and integrating the energy equation over the ionized region, we obtain an equation for the motion of the recombination wavefront $`x_i(t)`$. Once a solution for $`x_i`$ is computed, the bolometric luminosity from the inner envelope is approximately given by (see Zampieri et al.  for details)

$$L_{r_i}=L+4\pi r_i^2v_i(aT_{eff}^4/2+\rho Q_{ion}),$$
(3)
where $`L=(4\pi caT_0^4)/(3\kappa \rho _0)[y^2d\psi /dy]_{y=1}R_0x_i\varphi (t)`$ is the diffusion luminosity, the second term on the right hand side is the total advection luminosity released because of the wavefront motion, $`v_i=\dot{x}_iR`$ is the wavefront velocity relative to the envelope gas and $`Q_{ion}`$ is the recombination energy per unit mass.

In the recombined region ($`r_irR`$), the deposition of gamma ray photon energy through Comptonization and photoelectric absorption of heavy elements is the dominant thermal and radiative process. The optical luminosity emitted in this region originates by the reprocessing of the gamma rays and becomes important only at the end of the recombination stage (that coincides with the end of the plateau), when the radioactive decay time of 56Co becomes comparable to the expansion timescale. An approximate expression for this luminosity can be obtained neglecting the internal energy and the $`PdV`$ work in the energy equation and integrating it over the interval $`r_irR`$. The total bolometric luminosity of the envelope during this phase is then

$$L_{tot}=L_{r_i}+4\pi \rho _0R_0^3X_{Ni,0}f(t)_{x_i}^1x^2\xi _{Ni}(x)𝑑x.$$
(4)
where $`X_{Ni,0}`$ and $`\xi _{Ni}(x)`$ are the central mass fraction and radial distribution of 56Ni, $`f(t)=[3.9\times 10^{10}e^{t/\tau _{Ni}}+7.2\times 10^9(e^{t/\tau _{Co}}e^{t/\tau _{Ni}})]`$ erg g-1 s-1, $`\tau _{Ni}=8.8`$ days and $`\tau _{Co}=111`$ days are the Nickel and Cobalt decay times. Following Arnett , we take $`\xi _{Ni}(x)=\psi (x)=\mathrm{sin}(\pi x)/(\pi x)`$ (see Figure 2, where the 56Ni abundance is shown as a function of $`m`$).

In the third stage, when the envelope has completely recombined, the 56Co radioactive decay time is shorter than the expansion timescale and the luminosity is given by the second term in equation (4) with $`x_i=0`$.

## 4 Results and Discussion

Following the discussion in § 2, we assume that the phase of SN 1997D at discovery is $`90100`$ days and that the early light curve and the evolution of the line velocity resemble closely those of SN 1999br. Figure 3 shows the evolution of the UBVRI luminosity, the velocity of metal (Sc II) lines and the continuum temperature of SN 1997D and SN 1999br along with the results of our semi-analytic model calculations. In the top panel, the observed and calculated light curves of SN 1987A are also shown for comparison. The continuum temperatures have been estimated fitting the spectra with a Planckian and are affected by a significant uncertainty ($`\pm 500`$ K) because of the severe line blanketing at short wavelengths ($`<4500`$ Å, see Figure 1). The line velocity and continuum temperature inferred from the model refer to the photospheric epoch. Hence the last two plots in Figure 3 are truncated at the end of the recombination phase ($`110`$ days).

Considering the approximations adopted in the model, the general agreement with observations is satisfactory. We emphasize that obtaining a simultaneous “fit” of the light curve and the evolution of the line velocity and continuum temperature is an essential requirement to obtain a meaningful estimate of the parameters of the ejected envelope. In particular, the long plateau, the apparent break in the line velocity profile at $`30`$ days and its fast decline impose rather severe constraints on the model. The break appears to be related to the onset of recombination, that terminates when the light curves plummets at $`110`$ days.

The parameters of the post-shock, ejected envelope required to reproduce the observations of SN 1997D and SN 1999br are listed in Table 1. Only the values of the radius $`R_0`$, mass $`M_{env}`$, 56Ni mass $`M_{Ni}`$ and outer velocity $`V_0`$ (and hence initial energy) of the envelope have been significantly varied. The other model parameters were maintained essentially equal to the values required to reproduce the observations of SN 1987A. In particular, a colour correction factor $`f_c=T_c/T_{eff}=1.11.3`$ has been adopted to reproduce the observed continuum temperature for all the supernovae. Small colour corrections are rather common in stellar and supernova atmospheres and are induced by distortions of the continuum caused by radiative transfer processes (e.g. scattering, wavelength dependent opacities; see e.g. Eastman, Schmidt & Kirshner ).

Interpolating simultaneously the light curve and the evolution of the line velocity and continuum temperature results in a rather robust “fit”. The estimated value of the initial thermal+kinetic energy of the ejecta $`E`$ indicates that both events were rather under-energetic compared to a typical Type II supernova. The inferred 56Ni mass of SN 1999br is extremely small, testifying that the energy available to produce and eject nucleosynthetic elements was very low. Furthermore, because of the low luminosity in the plateau stage, the post-explosion envelope is rather compact. The ejected envelope masses are quite large, comparable to those required to reproduce the plateau of typical Type II supernovae. In particular, the ejected envelope mass of SN 1997D is almost three times larger than that estimated by Chugai & Utrobin  and only 30% smaller than that inferred by Turatto et al. . It is worth noting that the gross properties of the light curve and line velocity of SN 1997D and SN 1999br can be roughly accounted for by SN 1987A-like parameters, but simply decreasing the expansion velocity $`V_0`$ (and hence the energy $`E`$) and 56Ni mass.

Although it is not straightforward to determine the error in the parameters, we estimate that the intrinsic uncertainty of the “fit” is not larger than $``$ 30% (see also Zampieri et al. ). Additional sources of systematic errors are related to the approximations introduced in the present analysis, in particular to the choice of the initial conditions. Although it is known that significant mixing occurred in SN 1987A (see e.g. Woosley ), this effect may be less pronounced in other supernovae. The light curve is very sensitive to the prescription for mixing and to the actual velocity distribution as a function of mass $`V(m)`$. In particular, the luminosity and duration of the plateau depend on the energy and mass of the high velocity, Hydrogen rich part of the envelope. If the innermost Helium and heavier elements layers did not mix appreciably with the Hydrogen envelope, they would not produce any observable effect (having very low energy). In this case, the estimated value of $`M_{env}`$ would simply refer to the Hydrogen envelope mass. Then, in general, $`M_{env}`$ represents a lower limit to the total mass of the ejecta. Furthermore, the estimate of the envelope mass is sensitive to uncertainties in the actual value of the gas opacity and the details of the recombination physics. Here we adopted $`\kappa =0.2`$ cm2 g-1, $`t_{rec,0}30`$ days and $`T_{eff}=40004400`$ K, values similar to those used for the “fit” of SN 1987A. Only a larger opacity and/or a delayed onset of recombination could in principle allow for smaller envelope masses, but this is not in agreement with what found for SN 1987A and other Type II supernovae (Zampieri et al. ). The fact that the line velocities at $`100`$ days are slightly larger than what predicted by the model (see Figure 3) seems to indicate that the assumption of spatially constant density is only approximately correct and that more mass is concentrated in the innermost, low velocity part of the envelope. This effect may reduce the estimated energy $`E`$ because the bulk of the kinetic energy is carried by the outer, high velocity layers.

We emphasize that the present analysis for SN 1997D is based on the assumption that it is an event similar to SN 1999br. However, although the inferred parameters of the ejecta of SN 1997D rely on this hypothesis, the estimates for SN 1999br are certainly valid because they are not affected by uncertainties on the early light curve, the duration of the plateau and the evolution of the line velocity.

Determining the mass of the progenitor is quite a difficult task and, unless a pre-discovery identification is available (see, e.g., Smartt et al. ), the inferred value is usually rather uncertain. Adding to the mass of the envelope reported in Table 1 the mass of the collapsed core ($`2M_{}`$), the progenitors of SN 1997D and SN 1999br have at least 19 and 16 $`M_{}`$, respectively. We stress that these are most probably lower bounds and the actual values of the progenitor masses are likely to be larger. Therefore, the present estimates situate SN 1997D and SN 1999br in the intermediate mass progenitor range and rule then out an origin from the low end of the mass range of core-collapse supernovae.

It should be noted that the results of the present analysis are essentially insensitive to the inclusion of Rayleigh scattering from neutral Hydrogen. In fact, towards the end of the plateau, $`\tau _{ra}\kappa _{ra}\rho r0.1`$, where $`r10^{15}`$ cm, $`\rho 10^{12}`$ g cm-3 and the Rayleigh scattering opacity $`\kappa _{ra}10^4`$ cm2 g-1. Therefore, the outer recombined region is essentially transparent to Rayleigh scattering.

Understanding why the energy of these peculiar Type II supernovae is low is of paramount importance in connection with the hydrodynamics of the explosion and the formation of the central compact object. Because the mass of the iron core before the explosion does not vary significantly with progenitor mass $`M_{}`$ (Woosley & Weaver ), the gravitational potential energy liberated during the collapse of the core is roughly independent of $`M_{}`$. Therefore, the energy of the ejecta depends mostly on neutrino reheating and the hydrodynamic interaction of the supernova shock with the surrounding material. To unbind a spherically symmetric envelope, the shock must overcome the ram pressure of the gas that started to accrete after the collapse of the core (see, e.g., Fryer ). With increasing progenitor mass both the ram pressure and the binding energy of the envelope become larger. Therefore, comparatively less energy remains available to heat up and accelerate the ejecta. The final luminosity and expansion velocity are then small and the resulting supernova is under-energetic (Zampieri ). In this respect, it is tempting to note that the difference between the inferred total (kinetic+thermal+binding) energy of SN 1987A ($`2.5\times 10^{51}`$ erg) and the ejected kinetic+thermal energy of SN 1997D is $`1.5\times 10^{51}`$ erg, comparable to the binding energy of the envelope in a $`2530M_{}`$ progenitor star (Woosley & Weaver ). Therefore, in SN 1997D, a large fraction of the gravitational potential energy liberated by the collapse of the core may have been spent in trying to unbind the massive envelope.

An important consequence of this scenario is that, because of the low kinetic energy acquired by the envelope, a large fraction of the stellar material (in particular the innermost layers that have smaller velocities) is likely to remain gravitationally bound to the core after shock passage and fall back onto it. This is confirmed also by the small amount of 56Co present in the ejecta. The fallback of stellar material may also turn the newly formed neutron star into a black hole. If this happens, the late time fallback onto the central black hole may give rise to detectable emission of radiation with a characteristic power-law decay (Zampieri et al. , Balberg et al. ). The detection of this emission in the late time light curve would provide the first direct evidence for the presence of a black hole in the site of its formation.

We note that the overall picture may be different for a supernova from a massive progenitor whose post-shock envelope retains a significant fraction of its initial angular momentum. In this case, the hydrodynamics of the explosion becomes more complex and may give rise to a high energy, asymmetric explosion (MacFadyen & Woosley ).

The present analysis indicates that the parameters of the ejecta of SN 1997D and SN 1999br are consistent with them being intermediate mass core-collapse supernovae. If the mass of their progenitors is sufficiently large, they could form a black hole remnant, giving rise to significant fallback and late-time accretion onto the central compact object. In order to confirm these findings, more detailed investigations using radiation hydrodynamic simulations and spectral synthesis calculations are being planned. Clearly, monitoring this type of supernovae from discovery to late phases is of the utmost importance.

ACKNOWLEDGMENTS

We thank the referee Stan Woosley for valuable comments. We acknowledge support from the Italian Ministry for Instruction, University and Research (MIUR) through grant Cofin MM02905817 and the Italian Space Agency (ASI) under grant ASI I/R/70/00. M. H. acknowledges support for this work by NASA through Hubble Fellowship grant HST-HF-01139.01-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555.


text:
# Critical conditions for a stable molecular structure

## Abstract

Here, we show how the molecular structure appears and becomes stable for supercritical physical conditions. In particular we consider, for the ammonia molecule in a gas, a simple model based on a standard non-linear double-well Schrödinger equation with a dissipative term and a term representing weak collisions.

The existence of a well defined molecular structure for a symmetric molecule is an old and ongoing problem in chemistry W , V . Various attempts have been made to make quantum mechanics agree with apparently stationary asymmetric states. It was clear from the beginning that the action of the environment would be the basic reason of this effect, since from a quantum mechanics point of view an isolated symmetric molecule has no structure. The point is to find the simplest model simulating the influence of the environment, and to exploit it.

In 1927 Hund H  understood that in the case of certain molecules, such as ammonia $`NH_3`$, the basic model is a double well one, this simple model is able to give the splitting and the inversion frequency of an isolated molecule, but not the structure related to localization. From the experimental point of view, experimental data BL  showing the decrease of the splitting, as a decrease of the inversion line (red shift), for increasing values of the gas pressure appeared in 1948. Theoretical models for inter-molecular interactions were considered, and they show an interesting fit of the data at small pressure M , BR .

Further researchers proposed some simple models in order to explain the localization phenomenon of a symmetric molecule D , CJ . Later, some works GMS1 , GS  have exploited the strong effect of metastability, both static and dynamic, on red shift and localization.

In 1995 it was shown GM  that a non-linear Schrödinger equation, previously suggested by Pratt P  for the hydrogen case, is able to give spontaneous symmetry breaking and the bifurcation of the ground state at the value

$`\mu _c=1`$ (1)
of the dimensionless non-linearity parameter $`\mu `$ defined below (which is expected to be monotically dependent on the gas pressure). Later, it became clear that this non-linear model is able to give the red shift and localization too. In particular, in V  and GMS  it is suggested that localization appears for $`\mu =2\mu _c`$; while in JPT  it is suggested that the critical value (1) found in GM  coincides with the vanishing of the inversion line, and corresponds to the pressure of about 1,7 atmospheres at room temperature nota1 .

In this letter we consider the non-linear double-well model as above GM  in order to obtain the stability of the localization for finite values of the parameters. Let us recall Wi  that actually we consider a substituted ammonia molecule $`NHDT`$, so that a localized state has one of the two possible chiral configurations. In the approximation of a two-level system, an invariant quantity (energy) appears, and we observe periodic motions of two kinds: vibrational periodic motions around just one well (inside the grey regions of Fig. 1 and 2) or beating periodic motions between the two wells. In order to have a more physical model, we add a dissipative term, related to the photonic emission in the radio frequency range, giving the decreasing of the energy. We also take into account the effect of weak collisions. As a result we find that a chiral configuration is unstable if the non-linearity parameter is less than $`3\mu _c`$. For larger values of this parameter, the chirality states become stable provided that the interval between two collisions is larger than the relaxation time, in order that when a collision occurs, the state be near an asymmetric stationary state.

In order to better understand the situation, consider Figs. 1 and 2 where the space of states is shown, as a sphere projected on a square by a Mercator map. The two coordinates $`z`$ and $`\theta `$ respectively represent the imbalance variable (which measures the localization) and the relative phase of the components of the state with respect to the two localized states. A collision actually changes $`\theta `$ but not $`z`$. Thus, for $`\mu `$ greater than $`3\mu _c`$ (see Fig. 2), the line defined by a value of $`z`$ equal to the one of an asymmetric stationary state is fully contained in the vibrational region. In contrast, for $`\mu <3\mu _c`$ (see Fig. 1) we have that a change of $`\theta `$, due to a collision, could shift an asymmetrical stationary state into the beating region.

The Hamiltonian for a single ammonia molecule takes the form $`H_0=\frac{\mathrm{}^2}{2m}\mathrm{\Delta }+V`$ where $`V`$ is a double-well potential invariant under a coordinate reflection

$`V=\text{P}V\text{P},`$
where P is the unitary symmetry operator representing the inversion of the $`n`$–th coordinate

$`(x^{},x_n)(x^{},x_n),x^{}=(x_1,\mathrm{},x_{n1}),n1.`$
The tunneling time through the inter-well potential barrier is inversely proportional to the energy splitting $`\mathrm{\Delta }E`$ between the odd– and even–parity eigenstates $`|`$ and $`|+`$ with energies $`\lambda _{}`$ and $`\lambda _+`$, and the solution of the unperturbed equation

$`i\mathrm{}{\displaystyle \frac{}{t}}|\psi (t)=H_0|\psi (t)`$
shows a beating motion between the two wells with period

$`\tau ={\displaystyle \frac{2\pi \mathrm{}}{\mathrm{\Delta }E}},\mathrm{\Delta }E=\lambda _{}\lambda _+.`$
The actual semi-classical parameter is the energy splitting $`\mathrm{\Delta }E`$, and we choose the units such that $`\mathrm{}=1`$ and $`\mathrm{\Delta }E1`$.

Let us consider the interaction of a single molecule with the other molecules of the gas. In the mean field approximation, the new Hamiltonian takes the form

$`H=H_0+W`$
where the term $`W`$ is given by the polarization of the external environment due to the presence of the single ammonia molecule itself. Because of the dissipative terms, the perturbation $`W`$ could be written by means of a non-linear and non-Hermitian term

$`W=\nu g\left(ϵ\text{1}+i\eta \text{P}\right),\nu =\psi (t)|g|\psi (t),ϵ,\eta <0,`$
where $`g(x)`$ is a given odd function, $`\text{P}g\text{P}=g`$, and where the parameters $`ϵ`$ and $`\eta `$ respectively measure the strength of the dipole-dipole and the dissipation interactions. As in the case of complex Ginzburg-Landau equations AK  and of Gross-Pitaevskii equations with weakly dissipative effects CV , we expect to observe vortex solutions here.

We underline that in the non-dissipative case, where $`\eta =0`$, then we have the conservation of the charge:

$`\psi (t)|\psi (t)=\psi (t)^2=1`$ (2)
and the conservation of the energy functional

$`(\psi )=\psi ^2+\psi |V|\psi +{\displaystyle \frac{1}{2}}ϵ\psi |g|\psi ^2.`$
In the dissipative case, where $`\eta <0`$, since the term $`W`$ is such that $`[\text{T}\text{P},W]=0`$, where $`\text{T}|\psi (t)=|\psi (t)^{}`$ is the time-reversed operator, we still have the conservation of the charge (2), and we expect that the energy functional decreases. The system then relaxes towards local minima of the functional.

Here, we make use of the two-level model where we approximate the wave-function $`|\psi (t)`$ by means of its projection on the two lowest states. In particular, in such an approximation the total wavefunction of the system may be expanded as

$`|\psi (t)=a_R(t)|R+a_L(t)|L,`$ (3)
where

$`|R={\displaystyle \frac{1}{\sqrt{2}}}\left(|++|\right),|L={\displaystyle \frac{1}{\sqrt{2}}}\left(|+|\right)`$
are the right and left hand-side states; they are such that $`\text{P}|R=|L`$. The normalization condition (2) on the wavefunction $`|\psi (t)`$ implies that $`|a_R|^2+|a_L|^2=1`$.

By substituting $`|\psi (t)`$ by (3) in the time-dependent Schrödinger equation

$`i\mathrm{}{\displaystyle \frac{}{t}}|\psi (t)=H|\psi (t)`$
it follows that the expansion coefficients $`a_R`$ and $`a_L`$ have to satisfy to the following system of ordinary differential equations

$`\{\begin{array}{ccc}i\dot{a}_R\hfill & =& \mathrm{\Omega }a_R\omega a_L+ϵ\nu ca_R+i\eta \nu ca_L\hfill \\ i\dot{a}_L\hfill & =& \mathrm{\Omega }a_L\omega a_Rϵ\nu ca_Li\eta \nu ca_R\hfill \\ \nu \hfill & =& c(|a_R|^2|a_L|^2)\hfill \end{array}`$ (7)
where

$`\omega ={\displaystyle \frac{1}{2}}(\lambda _{}\lambda _+),\mathrm{\Omega }={\displaystyle \frac{1}{2}}(\lambda _{}+\lambda _+)`$
and where we set

$`c=R|g|R=L|g|L.`$
In order to study the beating motion, it is convenient to introduce the relative phase

$`\theta =\mathrm{arg}(a_R)\mathrm{arg}(a_L),`$
which is a torus variable, and the imbalance variable

$`z=|a_R|^2|a_L|^2,`$
which takes values in the interval $`[1,+1]`$. They have to satisfy to the system of ordinary differential equation

$`\{\begin{array}{ccc}\dot{z}\hfill & =& \omega Z(z,\theta )\hfill \\ \dot{\theta }\hfill & =& \omega \mathrm{\Theta }(z,\theta )\hfill \end{array}`$ (10)
where

$`Z(z,\theta )=2\sqrt{1z^2}\left[\mathrm{sin}\theta \zeta z\mathrm{cos}\theta \right]`$
and

$`\mathrm{\Theta }(z,\theta )=2{\displaystyle \frac{z}{\sqrt{1z^2}}}\left[\mathrm{cos}\theta +\zeta z\mathrm{sin}\theta \right]+2\mu z`$
and

$`\mu =c^2ϵ/\omega ,\zeta =c^2\eta /\omega .`$
Here, $`\mu `$ and $`\zeta `$ represent positive dimensionless parameters that measure the effective non-linearity and the dissipation, respectively.

We observe the existence of the critical value (1) for the non-linearity parameter $`\mu `$. We consider, at first, the weak non-linearity case such that $`\mu <\mu _c`$. In such a case, equation (10) admits just two stationary solutions corresponding to the unperturbed even– and odd–parity eigenstates: $`z_1=0`$ and $`\theta _1=0`$, corresponding to the even–parity eigenstate, is a stable stationary solution, and $`z_2=0`$ and $`\theta _2=\pi `$, corresponding to the odd-parity eigenstate, is an unstable stationary solution.

Then, we consider the strong non-linearity case such that $`\mu >\mu _c`$. We observe that when $`\mu `$ takes the value $`\mu _c`$ then the stable stationary solution makes experience of a bifurcation phenomenon GM . More precisely, for $`\mu >\mu _c`$ we have 4 stationary solution; two of them still correspond to the unperturbed even– and odd-parity eigenstate, the other two correspond to asymmetrical states that, in the limit of large non-linearity, are fully localized on just one of the two wells. The stationary solution $`z_1=0`$ and $`\theta _1=0`$, corresponding to the even–parity eigenstate, is a saddle point for $`\mu >\mu _c`$; the stationary solution $`z_2=0`$ and $`\theta _2=\pi `$, corresponding to the odd–parity eigenstate, is still an unstable solution; $`z_3=\sqrt{(\mu ^21)/(\mu ^2+\zeta ^2)}`$ and $`\theta _3=\text{arctan}(\zeta z_3)`$, and $`z_4=z_3`$ and $`\theta _4=\theta _3`$ are stable asymmetrical stationary solutions. Thus any state generically goes near to one of these two asymmetrical stationary states.

Therefore, in the strong non-linearity case we have that any initial state, except the two even– and odd-parity unperturbed eigenstates, finally goes to one of the two asymmetric stationary eigenstates giving a chiral configuration for the ammonia molecule.

Now, we show that this chiral configuration is stable with respect to collisions when the pressure is large enough. We underline that in the ammonia case the thermal energy at room temperature is smaller than the distance between the doublet $`\{\lambda _\pm \}`$ and the other energy levels, so that the validity of the two-level approximation holds, and it is much larger than the splitting energy, so that a collision could produce a strong variation of the energy $`(\psi )`$.

To this end we introduce a simplified model for molecular collision. When the single molecule undergoes a collision we add to the Hamiltonian $`H`$ a perturbative term of the type $`f(x)v(t)`$ where $`f(x)`$ is a function with compact support and $`v(t)`$ is a given time dependent function. For instance, let $`v(t)=\chi _{t_1,t_2}(t)`$ where $`\chi `$ is the characteristic function on the interval $`[t_1,t_2]`$, and where we assume that the perturbation acts for a time much shorter than the beating period, that is

$`t_2t_1\tau ,`$ (11)
Equation (7) then takes the form

$`\{\begin{array}{ccc}i\dot{a}_R\hfill & =& \mathrm{\Omega }a_R\omega a_L+ϵ\nu ca_R+i\eta \nu ca_L+v(t)c_Ra_R\hfill \\ i\dot{a}_L\hfill & =& \mathrm{\Omega }a_L\omega a_Rϵ\nu ca_Li\eta \nu ca_R+v(t)c_La_L\hfill \\ & & c_R=R|f|R,c_L=L|f|L\hfill \end{array}`$
since

$`R|f|L=\overline{L|f|R}0`$
Hence, the system (10) takes the form

$`\{\begin{array}{ccc}\dot{z}\hfill & =& \omega Z(z,\theta )\hfill \\ \dot{\theta }\hfill & =& \omega \mathrm{\Theta }(z,\theta )+(c_Rc_L)v(t)\hfill \end{array}`$
from which it follows that $`|z(t_2)z(t_1)|1`$, since (11), and where it is not possible obtain a similar bound for $`\theta `$ when $`c_Rc_L`$.

This fact, that is the relative phase is strongly modified after a generic collision (such that $`c_Rc_L`$), does not actually destroy the chiral configuration of a localized ammonia molecule if the pressure is large enough, i.e. such that $`\mu >3\mu _c`$. Indeed, in such a case we have that the stable solution $`(z_3,\theta _3)`$ (respectively $`(z_4,\theta _4)`$) has a basin of attraction containing the strip $`zz^{}`$ (respectively $`zz^{}`$) if $`\zeta `$ is small enough and where $`z^{}=2\sqrt{2}/3<z_3`$. We explain this fact by means of a continuity argument in the limit of $`\zeta =0`$. Indeed, in such a limit we have the existence of two separatrix lines GMS , AV  starting from the stationary solution $`(z_1,\theta _1)`$ and satisfying the equation

$`\sqrt{1z^2}\mathrm{cos}\theta +{\displaystyle \frac{1}{2}}\mu z^2=1.`$
It is not hard to see (Fig. 2) that these paths are contained in the strip $`z^{}zz^{}`$ if $`\mu >3\mu _c`$. As a result, it follows that a perturbation due to a collision acting in an interval of the order (11) shifts a state initially near to one of the asymmetric stationary stable eigenstates, to a state belonging to the basin of attraction of the stationary eigenstate itself. In particular, the state will be always far enough from the unstable stationary state and thus, after a finite time depending on $`\zeta `$, it returns near to the initial stationary state without visiting the other well, provided that in this period another collision does not occur.

Finally, let us notice that the molecular structure is not completely stable for strong collisions; indeed we have the phenomenon of racemization which makes the statistical mean of the chirality vanish at large time, and that an external electromagnetic field in the radio frequency range could be able to destroy the molecular structure.

In conclusion, in this letter we have shown that a dissipative non-linear model is able to explain the molecular structure of symmetric molecules. Such chiral configurations are stable for weak collisions provided that the non-linearity parameter is larger than the critical value $`3\mu _c`$, and that the frequency collisions is small enough. The kind of trap we propose here is simple, but not trivial, and in any case is able to make stable the spontaneous symmetry breaking given by the non-linearity. It could be relevant in the theory of dechoerence, related to the appearance of classical mechanics, and in the study of many irreversible precesses. As it clearly appears in this letter, we have only considered the spontaneous symmetry breaking effect. In fact, in the case of organic molecules, one enantiomer may be dominant. This effect could be explained by means of a small initial enantiomeric excess, due to the parity violation for weak interactions Q , largely amplified during a very long time by the dynamics of the system K .

###### Acknowledgements.

This work is partially supported by the Italian MURST, INFN and INDAM-GNFM. We would thank Prof. André Martinez for helpful discussions and remarks.


text:
# References

physics/0210117

SNBNCBS-2002

Jacobi Identity for Poisson Brackets: A Concise Proof

R.P.Malik *** E-mail address: malik@boson.bose.res.in

S. N. Bose National Centre for Basic Sciences,

Block-JD, Sector-III, Salt Lake, Calcutta-700 098, India

Abstract: In view of the recent interest in a short proof of the Jacobi identity for the Poisson-brackets, we provide an alternative simple proof for the same. Our derivation is based on the validity of the Leibnitz rule in the context of dynamical evolution.

PACS number(s): 45.20.-d

Keywords: Poisson brackets; Jacobi identity; Leibnitz rule; dynamical evolution

In the context of a precise classical description of the particle dynamics, the Poisson brackets (PBs) play a very prominent role as far as the Hamiltonian formulation of the particle mechanics is concerned (see, e.g., ). The bilinearity, antisymmetry property and celebrated Jacobi identity, etc., are some of the key properties that are respected by PBs. Many of the text books on classical mechanics provide the proof of Jacobi identity by exploiting the tedious (but straightforward) algebra connected with the basic definition of the PB in the (momentum) phase space . “There seems to be no simple way of proving Jacobi’s identity for the Poisson bracket without lengthy algebra” says Goldstein in his famous book . To mitigate the complexity of manipulations, a certain specific set of differential operators are defined in the well-known text books (see, e.g., ). In a recent article , a short proof for this identity has been given by using the concept of infinitesimal canonical transformations for the dynamical variables. The key role in this derivation is played by the generator formalism  which owes its origin to the innate group property associated with the above infinitesimal transformations. The purpose of the present note is to demonstrate that there is yet another simple proof for the Jacobi identity where (i) mainly the Leibnitz rule is exploited for the derivation, and (ii) the emphasis is laid on the dynamical evolution of the system. The Leibnitz rule plays a very important and decisive role in the Hamiltonian description of the particle dynamics. In fact, a dynamical system is said to be Hamiltonian  if and only if the time derivative acts on a PB as if the latter were a product of two dynamical variables (see, e.g., (2) below).

Let $`f(q,p)`$ and $`g(q,p)`$ be a couple of dynamical variables, defined in the (momentum) phase space. It is evident that, right from the beginning, these variables have no explicit time dependence. Thus, their time evolution w.r.t. the Hamiltonian function $`H(q,p)`$ is

$$\begin{array}{ccc}\frac{df}{dt}=\frac{f}{q_i}\dot{q}_i+\frac{f}{p_i}\dot{p}_i=\frac{f}{q_i}\frac{H}{p_i}\frac{f}{p_i}\frac{H}{q_i}\{f,H\}_{PB},\frac{dg}{dt}=\{g,H\}_{PB},\hfill & & \end{array}$$
$`(1)`$
where $`i=1,2,3\mathrm{}.s`$ corresponds to the $`s`$ number of degrees of freedom associated with the mechanical system. Here the summation convention is adopted for the definition of the basic PB in the phase space characterized by $`(q_i,p_i)`$. The application of the Leibnitz rule leads to

$$\begin{array}{ccc}\frac{d}{dt}\{f,g\}_{PB}=\{\frac{df}{dt},g\}_{PB}+\{f,\frac{dg}{dt}\}_{PB}.\hfill & & \end{array}$$
$`(2)`$
The above equation is valid even for the case when $`f`$ and $`g`$ are explicitly dependent on time (see, e.g., ). Now exploiting the basic definition (1) of the time evolution for a dynamical variable, it can be seen that each term of (2) can be separately expressed as

$$\begin{array}{ccc}\frac{d}{dt}\{f,g\}_{PB}\hfill & =& \{\{f,g\}_{PB},H\}_{PB},\{\frac{df}{dt},g\}_{PB}=\{\{f,H\}_{PB},g\}_{PB},\hfill \\ \{f,\frac{dg}{dt}\}_{PB}\hfill & =& \{f,\{g,H\}_{PB}\}_{PB}.\hfill \end{array}$$
$`(3)`$
Substitutions of the above expressions into (2) and rearrangements (corresponding mainly to the antisymmetry property of the PB) yield the Jacobi identity for all the above three dynamical variables as

$$\begin{array}{ccc}\{\{f,g\}_{PB},H\}_{PB}+\{\{g,H\}_{PB},f\}_{PB}+\{\{H,f\}_{PB},g\}_{PB}=0.\hfill & & \end{array}$$
$`(4)`$
As far as the emphasis on the evolution of the dynamical variable is concerned, the derivation for the Jacobi identity in this note highlights the key role played by the Leibnitz rule in the context of Hamiltonian description of the particle mechanics. Of course, here the Hamiltonian $`H`$ plays similar kind of exceptional role as the generator $`h`$ plays for the infinitesimal canonical transformation $`\delta _c`$. For the derivation of Jacobi identity in , one starts with the application of $`\delta _c`$ on a PB as

$$\begin{array}{ccc}\delta _c\{f,g\}_{PB}=\{\{f,g\}_{PB},h\}_{PB}.\hfill & & \end{array}$$
$`(5)`$
The l.h.s. of the above equation can be further expanded as $`\{\delta _cf,g\}_{PB}+\{f,\delta _cg\}_{PB}`$. Now, exploiting the basic definition of the generator $`h`$ (i.e. $`\delta _cf=\{f,h\}_{PB}`$) and making some rearrangement in the above equation, one derives the Jacobi identity as

$$\begin{array}{ccc}\{\{f,g\}_{PB},h\}_{PB}+\{\{g,h\}_{PB},f\}_{PB}+\{\{h,f\}_{PB},g\}_{PB}=0.\hfill & & \end{array}$$
$`(6)`$
It should be noted that the generator $`H`$ for the dynamical evolution and the generator $`h`$ for the infinitesimal canonical transformation are on equal footing as far as the derivation of the Jacobi identity in (4) and (6) is concerned. Whereas the emphasis in the former case is on the “dynamical evolution”, the latter case revolves around the concept of “infinitesimal symmetry” and its innate group property.