text:
# Driven low density granular mixtures.

## I Introduction

Granular materials present a rich and intriguing phenomenology, which has attracted the interest of the scientific community since the nineteenth century general . However, in spite of its recent progress the theoretical study of granular gases, i.e. of fluidized granular particles, is certainly less advanced than that concerning ordinary molecular fluids. The reason for this state of affairs is the presence of dissipation due to inelastic collisions and of friction with the surroundings, which prevents these system to reach thermodynamic equilibrium. In fact, in order to render stationary a granular system one needs to inject continuously energy into the system. This can be done, for instance by shaking or vibrating the grains.

In the present paper we illustrate the results of a numerical investigation concerning the properties of a two component granular mixture, modeled, following an established tradition, by inelastic hard spheres (IHS) with different masses, restitution coefficients, radii, and subject to different forms of external drive. The physical motivation for our study stems from the fact that in nature most granular materials are polydisperse from the point of view of their sizes and/or of their physical and mechanical properties. The theoretical study of granular mixtures has attracted so far the attention of several researchersdufty ; ioepuglio ; Barrat ; duftynew . These studies comprehend both freely cooling and uniformly heated granular mixtures and have been performed almost contemporaneously with laboratory experiments feitosa ; Parker . The most striking outcome is the lack of energy equipartition, i.e. the presence of two different kinetic temperatures, one for each species.

In the present paper our interest will be concentrated on the stationary state obtained by applying an energy feeding mechanism, represented by a uniform stochastic driving force or by a vibrating wall. In the first case the work of Barrat and Trizac Barrat  and Garzó et al. duftynew  provides a simple theoretical account, based on the Boltzmann equation approach and on the Direct Simulation Monte Carlo method (DSMC) bird1 ; bird2 , of the lack of energy equipartition.

The present work differs from the recently appeared papers dealing with granular mixtures in several aspects:

* the heat bath, which drives the system, is at finite temperature, since we have introduced a finite friction of the particles with the surroundings and we have considered the possibility of having spatially inhomogeneous states;
* in addition, we considered a second situation in which the particles subject to a vertical gravitational field receive energy anisotropically from the bottom vibrating wall of the container;
* in both cases we provide information about the presence of inhomogeneities in the system, e.g. density clustering and non-uniform density and temperature profiles.

In Sec. II we present the model fluid and two different mechanisms of energy supply. In sect III we discuss the first sub-model, the one with the heat bath, and obtain by means of an approximate analytic method an estimate of the partial temperature of each component. Subsequently we study the same sub-model with a DSMC algorithm. In sect. IV we study numerically the second sub-model (the one with gravity and vibrating wall) by means of DSMC simulations. Finally in the last section, V, we present our conclusions.

## II Definition of the models

We shall consider a dilute inelastic gas constituted of $`N_1`$ particles of mass $`m_1`$ and $`N_2`$ particles of mass $`m_2`$ subject to some kind of external driving (this will be specified in the following). We suppose that the interactions between the grains can be described by the smooth inelastic hard sphere model (IHS) campbell , thus we specify only the radius of the spheres, their masses and the fraction of the kinetic energy dissipated at each collision. This can be done by defining three different restitution coefficients $`\alpha _{ij}`$, i.e. $`\alpha _{11},\alpha _{22}`$, and $`\alpha _{12}=\alpha _{21}`$, which account for normal dissipation in collisions among particles of type $`i`$ and $`j`$. No internal degrees of freedom (e.g. rotations) are included.

One can describe the velocity changes induced by the instantaneous inelastic collisions of smooth disks labeled $`1`$ and $`2`$ of diameter $`\sigma _1`$ and $`\sigma _2`$ by the following equations:

$`𝐯_1^{}`$ $`=𝐯_1{\displaystyle \frac{1+\alpha _{\kappa _1\kappa _2}}{2}}{\displaystyle \frac{m_{\kappa _2}}{m_{\kappa _1}+m_{\kappa _2}}}((𝐯_1𝐯_2)\widehat{𝐧})\widehat{𝐧}`$ (1a)
$`𝐯_2^{}`$ $`=𝐯_2+{\displaystyle \frac{1+\alpha _{\kappa _1\kappa _2}}{2}}{\displaystyle \frac{m_{\kappa _1}}{m_{\kappa _1}+m_{\kappa _2}}}((𝐯_1𝐯_2)\widehat{𝐧})\widehat{𝐧}`$ (1b)
where $`\widehat{𝐧}=2(𝐱_\mathrm{𝟏}𝐱_\mathrm{𝟐})/(\sigma _{\kappa _1}+\sigma _{\kappa _2})`$ is the unit vector along the line of centers $`𝐱_1`$ and $`𝐱_2`$ of the colliding disks at contact and $`\kappa _1,\kappa _2`$ are the species ($`1`$ or $`2`$) to whom particles $`1`$ and $`2`$ belong. An elementary collision conserves the total momentum and reduces the relative kinetic energy by an amount proportional to $`(1\alpha _{\kappa _1\kappa _2}^2)/4`$. The collision rule we have adopted excludes the presence of tangential forces, and hence the rotational degrees of freedom do not contribute to the description of the dynamics.

Since the particles suffer mutual collisions and loose kinetic energy, in order to achieve a steady state, one needs to supply from the exterior some energy. The energy source has been modeled in two different fashions.

In the first sub-model we have assumed that the particles, besides suffering mutual collisions experience a uniform stochastic force and a viscous dampingpuglisi1 ; puglisi2 . The presence of the frictional, velocity-dependent term in addition to the random forcing, not only is motivated by the idea of preventing the energy of a driven elastic system ($`\alpha _{\kappa _1\kappa _2}1`$), to increase indefinitely, but also mimics the presence of friction of the particles with the container. Moreover a fluctuation dissipation relation is assumed between the viscous force and the intensity of the noise. Even in extended systems with small inelasticity the absence of friction may cause some problems of stability tobepubl . We shall be mainly interested in the stationary situation determined by the balance between the energy feeding mechanism and the dissipative forces due both to the inelasticity of the collisions and to the velocity-dependent frictional force.

In the second sub-model the grains are constrained to move on a frictionless inclined plane and the bottom boundary vibrates (as a thermal japan  or deterministic kudrolli2  oscillating wall) giving energy only to the particles bouncing on it. Periodic boundary conditions are assumed laterally.

Since we consider throughout only sufficiently low density systems successive binary collisions are effectively uncorrelated and Boltzmann equation can be used to describe the non equilibrium dynamics. In fact, by assuming the validity of the Boltzmann molecular chaos hypothesis, introduced to treat the collision term, it is straightforward to derive the governing equations for the probability density distribution function of each species. Moreover in this situation we are allowed to perform Direct Simulation Monte Carlo to integrate numerically the inhomogeneous Boltzmann equation of the sub-models.

## III Uniformly heated system

In order to see the effect of the heat bath let us consider the system in the absence of collisions. In this case, the evolution of the velocity of each particle is described by an Ornstein-Uhlenbeck process. If we require that the two components must reach the same granular temperature in the limit of vanishing inelasticity we have two different possibilities to fix the heat bath parameters:

$$_t𝐱_i(t)=𝐯_i(t)$$
(2)
$`m_i_t𝐯_i(t)`$ $`=\gamma 𝐯_i(t)+\sqrt{2\gamma T_b}\xi _i(t)`$ (3a)
$`m_i_t𝐯_i(t)`$ $`=m_i\beta 𝐯_i(t)+\sqrt{2m_i\beta T_b}\xi _i(t)`$ (3b)
where $`i=1,2`$ and $`T_b`$ is the heat bath temperature and $`\xi (t)`$ is a Gaussian noise with the following properties:

$`<\xi _i(t)>`$ $`=0`$ (4a)
$`<\xi _i(t_1)\xi _j(t_2)>`$ $`=\delta (t_1t_2)\delta _{ij}`$ (4b)
The associated Fokker-Planck equations for the two cases are respectively:

$`_tf_i(𝐫,𝐯,t)`$ $`={\displaystyle \frac{\gamma }{m_i}}_v(𝐯f_i(𝐫,𝐯,t))+{\displaystyle \frac{\gamma T_b}{m_i^2}}_v^2f_i(𝐫,𝐯,t)+𝐯_rf_i(𝐫,𝐯,t)`$ (5a)
$`_tf_i(𝐫,𝐯,t)`$ $`=\beta _v(𝐯f_i(𝐫,𝐯,t))+{\displaystyle \frac{\beta T_b}{m_i}}_v^2f_i(𝐫,𝐯,t)+𝐯_rf_i(𝐫,𝐯,t)`$ (5b)
### III.1 Spatially Uniform solutions

When we take into account collisions among particles equations (5) become two coupled Boltzmann equations modified by the presence of a diffusion term due to the thermal noise. We shall first consider a spatially homogeneous case in order to derive the temperature of each species in the homogeneous stationary state. This represents a sort of mean field approximation. In fact, in the original Boltzmann equations the collisions are considered to occur only between particles at the same spatial location, whereas here these can occur between arbitrary pairs of particles regardless their spatial separation. In the present subsection we compute the second moments, $`<v_i^2>`$, of the distribution functions in order to determine the partial temperatures and their ratio. Although the method of derivation of the equations for the partial temperatures is not original, we present it in order to render the paper self-contained and because it shows the differences between the particular heat bath we employed and those chosen by other authors Barrat . First, indicating by $`n_i=N_i/V`$ the partial density of species $`i`$, we notice that both eqs. (5) possess the same equilibrium solution:

$$f_i(𝐯)=n_i\left(\frac{m_i}{2\pi T_b}\right)^{\frac{d}{2}}e^{\frac{m_iv^2}{2T_b}}$$
(6)
although the relaxation properties are different. Only upon adding the inelastic collision term the two species display different temperatures. The resulting Boltzmann equation for a granular mixture is dufty ; Barrat ; ioepuglio :

$$_tf_i(𝐯_1;t)=\underset{j}{}J_{ij}[𝐯_1|f_i,f_j]+\frac{\xi _{0i}^2}{2}_v^2f_i+\beta _i_v(𝐯_1f_i)$$
(7)
where we have used a compact notation to represent the two different choices of heat bath:

* Case 1

$$\begin{array}{c}\hfill \xi _{0i}^2\frac{2\gamma T_b}{m_i^2}\\ \hfill \beta _i\frac{\gamma }{m_i}\end{array}$$
(8)
* Case 2

$$\begin{array}{c}\xi _{0i}^2\frac{2\beta T_b}{m_i}\\ \beta _i\beta \end{array}$$
(9)

and $`J_{ij}[v_1|f_i,f_j]`$ is the collision integral:

$$J_{ij}[v_1|f_i,f_j]\sigma _{ij}^2𝑑𝐯_2𝑑\widehat{\sigma }\mathrm{\Theta }(\widehat{\sigma }𝐠_{12})(\widehat{\sigma }𝐠_{12})[\alpha _{ij}^2f_i(𝐯_\mathrm{𝟏}^{})f_j(𝐯_\mathrm{𝟐}^{})f_i(𝐯_\mathrm{𝟏})f_j(𝐯_\mathrm{𝟐})]$$
(10)
The primed velocities are pre-collisional states, which can be obtained by inverting eqs. (1)

Due to the presence of the heat bath terms the system reaches asymptotically a steady state, characterized by time independent pdf’s. By imposing the vanishing of the time derivatives we obtain:

$$\underset{j}{}J_{ij}[𝐯_1|f_i,f_j]+\frac{\xi _{0i}^2}{2}_v^2f_i+\beta _i_v(𝐯_1f_i)=0$$
(11)
or after integrating over $`𝐯_\mathrm{𝟏}`$:

$$\underset{j}{}𝑑𝐯_1v_1^2J_{ij}[𝐯_1|f_i,f_j]+\frac{\xi _{0i}^2}{2}𝑑𝐯_1v_1^2_v^2f_i+\beta _i𝑑𝐯_1v_1^2_v(𝐯_1f_i)=0$$
(12)
After simplifying the second and the third integral by integration by parts and using the normalization property $`f_i𝑑v_i=n_i`$ we find:

$$\underset{j}{}𝑑𝐯_\mathrm{𝟏}v_1^2J_{ij}[𝐯_1|f_i,f_j]+n_id\xi _{0i}^22\beta _i𝑑𝐯_1v_1^2f_i(𝐯_1)=0$$
(13)
By defining the partial temperature:

$$n_iT_i\frac{1}{d}𝑑𝐯_\mathrm{𝟏}m_iv_1^2f_i$$
(14)
we can rewrite eq. (13) in terms of the temperatures :

$$T_i=\frac{m_i}{2d\beta _i}\left(\frac{1}{n_i}\underset{j}{}𝑑𝐯_\mathrm{𝟏}v_1^2J_{ij}[v_1|f_i,f_j]+d\xi _{0i}^2\right)$$
(15)
Eq. (15) determines the partial temperatures once the $`f_i`$ are known. In practice one can obtain an estimate of $`T_i`$ by substituting two Maxwell distributions:

$$f_i(v)=n_i\left(\frac{m_i}{2\pi T_i}\right)^{\frac{d}{2}}e^{\frac{m_iv_1^2}{2T_i}}$$

and performing the remaining integrals (see dufty ; Barrat ) one gets:

$$\begin{array}{cc}& \frac{d\mathrm{\Gamma }(d/2)}{m_i\pi ^{(d1)/2}}2\beta _i(T_bT_i)=\sigma _{ii}^{d1}n_i\frac{2(1\alpha _{ii}^2)}{m_i^{3/2}}T_i^{3/2}\hfill \\ & +\sigma _{ij}^{d1}n_j\mu _{ji}\left[\mu _{ji}(1\alpha _{ij}^2)\left(\frac{2T_i}{m_i}+\frac{2T_j}{m_j}\right)+4(1+\alpha _{ij})\frac{T_iT_j}{m_1+m_2}\right]\left(\frac{2T_i}{m_i}+\frac{2T_j}{m_j}\right)^{1/2}\hfill \end{array}$$
(16)
where $`\mu _{ij}=m_i/(m_i+m_j)`$.

By solving numerically the nonlinear system of eqs. (16) one obtains the steady values of the partial temperatures in the spatially homogeneous situation.

In figs. 1 and 2 we report the temperature ratio $`\frac{T_1}{T_2}`$ as a function of a common restitution coefficient $`\alpha `$, having chosen equal coefficients $`\alpha _{11}=\alpha _{22}=\alpha _{12}=\alpha `$. Assuming identical concentrations, and a varying the mass ratio $`\frac{m_1}{m_2}`$, we studied the cases 1 and 2.

We notice that with the second recipe (case 2) the temperature ratio is an increasing function of the mass ratio $`m_1/m_2`$: this is exactly the opposite of what happens for the case 1. The experimental observation feitosa  suggest that the trend of case 2 is the right one. In other words, in the case 2 both the friction term and the power supplied are proportional to the mass of the two species and such a property is probably true also in the experiments. This is why we shall use the recipe of case 2 in the following DSMC simulations.

### III.2 Non uniform solutions

In the present section we illustrate the results obtained simulating the system with the heat bath (with recipe 2, i.e. equation (3b)) by the so called Direct Simulation Monte Carlo according to the implementation described in Baldassa .

At every time step of length $`\mathrm{\Delta }t`$ each particle is selected to collide with a probability $`p_c=\mathrm{\Delta }t/\tau _c`$ (where $`\tau _c`$ is an a priori fixed mean free time established consistently with the mean free path and mean squared velocity) and seeks its collision partner among the other particles in a neighborhood of radius $`r_B`$, choosing it randomly with a probability proportional to their relative velocities. Moreover in this approximation the diameter $`\sigma `$ is no more explicitly relevant, but it is directly related to the choices of $`p_c`$ and $`r_B`$ in a non trivial way: in fact the Bird algorithm allows the particles to pass through each other, so that a rigorous diameter cannot be defined or simply estimated as a function of $`p_c`$ and $`r_B`$. In this section, to indicate the degree of damping, we give the time $`\tau _b=1/\beta `$ instead of $`\beta `$. This is useful to appreciate the ratio between the mean collision time $`\tau _c`$ and the mean relaxation time due to the bath which indeed is $`\tau _b`$.

In the present section we chose $`N_1=N_2=500`$ and $`T_b=1`$, and equal restitution coefficients for all collisions and $`\tau _c=0.16`$. As illustrated in figures 3 and 4 the two components display different granular temperatures in agreement with the analytical predictions of the homogeneous Boltzmann equations. We also comment that the ratio $`T_1/T_2`$ is very sensitive to the mass ratio and much less to differences in the restitution coefficients and concentrations. This is in agreement with experimental observations feitosa .

At a finer level of description we consider the rescaled velocity pdf for different values of the inelasticity parameter and different mass ratios $`m_2/m_1=2`$. One sees that not only the deviations from the Gaussian shape become more and more pronounced as we increase the inelasticity parameter $`(1\alpha )`$, but also the shape of the two distributions differ appreciably in the tails even after velocity rescaling to make the two pdf’s have the same variance. Such a property is similar to the one already reported in ioepuglio . One can also observe that the rescaled pdf of the lighter species has slightly broader tails. The mass ratio also controls the deviations from the Gaussianity of the velocity pdf’s. Whereas it is known that the departure from the Maxwell-Boltzmann statistics is triggered by the inelasticity of the collisions, the larger the inelasticity the stronger the deviation, to the best of our knowledge this is the first numerical evidence of the phenomenon, predicted within a Maxwell model in ref. ioepuglio . Comparing figs. 5 and 6 one sees that the mass asymmetry enhances the non gaussianity of the pdf.

We have also studied the limits of low and high $`\tau _b`$, in figure 7, to show how the velocity distributions change. For values of the characteristic time of the heat bath, $`\tau _b`$, comparable with the collision time, $`\tau _c`$, the dynamics is essentially controlled by the stochastic acceleration term. This fact renders the two partial temperatures very close and makes the velocity distributions nearly Maxwellian. As $`\tau _b`$ increases we have observed that the energy dissipation due to the inelasticity makes the temperatures of the two species different. Moreover, the temperature ratio displays power law decreasing trend as a function of $`\tau _b`$ whose strength depends on the mass ratio (see figure 8).

In order to obtain information about the spatial structure of the mixture we have performed an analysis of the following correlation function, already introduced in the context of granular media by puglisi1 ; puglisi2 ; procaccia :

$$C_{\alpha \beta }(r)=\frac{1}{N(N1)}\underset{ij}{}\mathrm{\Theta }(r|𝐱_i^\alpha 𝐱_j^\beta |)$$
(17)
For a spatially homogeneous system we expect that $`C_{\alpha \beta }(r)r^{d_2}`$ with $`d_2=d`$ the dimension of the embedding space. However, in the cases we have studied, we found that $`d_2<d`$ (see fig 9), a signature that the system tends to clusterize. The larger the value of $`\tau _b`$ the stronger the deviation.

## IV System with gravity and heated from below

We turn, now, to illustrate the results relative to an inelastic mixture subject to gravity and confined to a vertical plane of dimensions $`L_x\times L_y`$. In the horizontal direction $`x`$ we assumed periodic boundary conditions. Vertically the particles are confined by walls. Energy is supplied by the bottom wall vibrating stochastically or periodically according to the method employed in Baldassa  for a one component system.

The vibration can have either a periodic character (as in ref. kudrolli ) or a stochastic behavior with thermal properties (as in japan ). In the periodic case, the wall oscillates vertically with the law $`Y_w=A_w\mathrm{sin}(\omega _wt)`$ and the particles collide with it as with a body of infinite mass with restitution coefficient $`\alpha _w`$, so that the vertical component of their velocity after the collision is $`v_y^{}=\alpha _wv_y+(1+\alpha _w)V_w`$ where $`V_w=A_w\omega _w\mathrm{cos}(\omega _wt)`$ is the velocity of the vibrating wall. In the stochastic case we assume that the vibration amplitude is negligible and that the particle colliding with the wall have, after the collision, new random velocity components $`v_x(\mathrm{},+\mathrm{})`$ and $`v_y(0,+\mathrm{})`$ with the following probability distributions:

$$P(v_y)=\frac{v_y}{T_w}exp(\frac{v_y^2}{2T_w})$$
(18)
$$P(v_x)=\frac{1}{\sqrt{2\pi T_w}}exp(\frac{v_x^2}{2T_w})$$
(19)
In this model we assume that the particle do not feel any external (environmental) friction. In all simulation we have chosen $`N_1=N_2=100`$, $`g=1`$, $`m_1=m_2/2=1/4`$, $`L_x=L_y=L=\sqrt{N}`$.

The Boltzmann equations for the partial distribution functions $`f_\kappa (𝐫,𝐯;t)`$ with $`\kappa =1,2`$ read

$$\left(\frac{}{t}+𝐯_r+g_i\frac{}{v_i}\right)f_\kappa (𝐫,𝐯,t)=\underset{\beta }{}J_{\kappa \beta }(f_\kappa ,f_\beta ).$$
(20)
The simulation results show that the tendency of the grains to form clusters is enhanced by such a choice of driving mechanism with respect to the homogeneous heat bath of the previous section. In the latter the noise acting uniformly was more effective in breaking the clusters. In addition the gravitational force tends to group the particles in the lower portion of the container for not too large driving frequencies Baldassa .

Fig 10 illustrates the partial density profiles and granular temperature profiles in the presence of a thermal wall of intensity $`T_p=20`$. The density profiles differ slightly near the bottom wall where both present a maximum (see inset). The temperature profiles are different, putting again in evidence a strong lack of equipartition in the system. Interestingly, the temperature ratio is almost constant along the vertical direction, notwithstanding the partial profiles are non constant. We have also performed numerical simulations with a harmonically vibrating wall: as the maximum velocity of the vibrating wall (which, for $`A=1`$, is equal to $`\omega `$) increases, the position of the density maximum raises indicating that gravity becomes less and less relevant. As far as the partial temperatures are concerned figure (11) illustrates the corresponding situation.

One sees that the temperatures (see inset) next to the wall attain their largest value, then drop to increase again. This indicates that the region far from the bottom is hotter because of the lower density of the gas and of the lower collision rate. As the wall temperature increases the temperature becomes more homogeneous far from the bottom due to the major homogeneity in density. We also notice a small segregation effect.

Finally, we present the velocity distributions measured in the system with a thermal wall. Since the temperature depends on the vertical coordinate $`y`$ we have computed the velocity pdf at different heights. The various pdf are plotted in figure (12) after a suitable rescaling. We observe that the distributions deviate appreciably from a Gaussian shape and display overpopulated high energy tails.

Interestingly the value obtained by our simulation for the temperature ratio, having chosen $`\frac{m_1}{m_2}=0.3`$ and a single restitution coefficient $`\alpha =0.93`$, $`\frac{T_1}{T_2}0.75`$ is not too far from the value obtained experimentally by Feitosa and Menon feitosa  which is $`0.66\pm 0.06`$. One should recall that in our “setup” only the lower wall supplies energy to the system, whereas in the experiment both walls vibrate.

## V Conclusions

To summarize we have studied the steady state properties of a granular mixture subject to two different classes of external drive.

In the case of an heat bath acting homogeneously on the grains, we have first obtained the temperatures of each species by employing an approximate analytic method. We have then solved numerically the equations by allowing the density and the temperatures to be spatially varying and analyzed the spatial correlations and the velocity distributions.

Finally, we turned to investigate the properties of the mixture in the presence of gravity and inhomogeneous drive.

The present results not only confirm the predictions of different temperatures for the two species for all types of drive, but also show the existence of different shape functions for the velocity distributions and of overpopulated high energy tails in agreement with the findings based on Maxwell models ioepuglio .

In order to obtain a direct comparison with experiments we employed parameters comparable to those utilized in the experimental work of ref. feitosa  and found temperature ratios not too far from the ones there reported. In spite of this success, the agreement seems to be limited to the average properties, while it is clear from the observation of the density profiles computed numerically that a better treatment of the excluded volume effect is necessary in order to obtain a more realistic description of the system. Finally, in view of the approximations inherent to the Boltzmann approach, we have not tried to include neither the effect of rotations of the grains nor the friction with the lateral walls which might be both relevant Paolotti .

###### Acknowledgements.

This work was supported by Ministero dell’Istruzione, dell’Università e della Ricerca, Cofin 2001 Prot. 200102384. A. P. acknowledges support from the INFM Center for Statistical Mechanics and Complexity (SMC).


format: text

text:
# Shear-Induced Stress Relaxation in a Two-Dimensional Wet Foam

## Abstract

We report on experimental measurements of the flow behavior of a wet, two-dimensional foam under conditions of slow, steady shear. The initial response of the foam is elastic. Above the yield strain, the foam begins to flow. The flow consists of irregular intervals of elastic stretch followed by sudden reductions of the stress, i.e. stress drops. We report on the distribution of the stress drops as a function of the applied shear rate. We also comment on our results in the context of various two-dimensional models of foams.

Foams are ubiquitous in nature  and part of a larger class of materials which exhibit a type of behavior that has come to be called “jamming” . Loosely speaking, a “jammed” material is one that is unable to flow, typically due to the packing of the constituent particles. For the case of foams, which are composed of gas bubbles separated by fluid walls, the jamming is a consequence of the topological constraints that develop as the bubbles press against each other. The existence of a jammed state contributes to the complex flow behavior of foams. This flow behavior is one reason that foams are so interesting from both a fundamental and applied point of view . For sufficiently small strains or stresses, a foam acts as an elastic solid (the jammed state). However, when the strain (or the stress) exceeds a critical value, known as the yield strain, a foam begins to flow. Of particular interest is the flow behavior of foams for small rates of strain. This flow occurs through “avalanches”, or sudden nonlinear rearrangements of the bubbles in the foam that correspond to a decrease in the average energy, or stress, of the foam. In this paper, we will refer to these events as “stress drops”. This type of behavior is common to many jammed systems where flow occurs in an irregular or stick-slip type manner when the system is near the transition to jamming . One of the outstanding questions for flowing foams is the nature of the distribution of the sizes of the stress drops during this irregular flow.

One reason that this remains an open issue is that predictions for the distribution of stress drop sizes are somewhat model dependent. Various models of two-dimensional foams have been used to study the connection between topological rearrangements and flow properties as a function of applied rate of strain . All of these simulations predict the same qualitative behavior for the macroscopic flow for sufficiently small rates of strain that was described above: elastic behavior below a yield strain and nonlinear bubble rearrangements above the yield strain. However, the quantitative predictions for the size distribution and frequency of these events are model dependent. The dependence is predominately due to different assumptions concerning the source of dissipation in foams and treating foams with different degrees of “dryness”. Distinguishing between these different models, and their underlying assumptions, experimentally is an important step in our understanding of foam rheology.

Two-dimensional foams may be characterized by the area fraction of gas, $`\varphi `$. For $`\varphi <0.84`$, foams “melt” into a froth of exclusively circular bubbles. Near this transition, where bubbles are predominately circular, a foam is said to be “wet”. In the limit $`\varphi `$ approaches 1, the bubbles become polygonal with infinitely thin walls, and a foam is said to be “dry”. The simulations studied by Kawasaki, et al.  are based on the vertex models of foams. This model assumes infinitely thin walls and only treats the dynamics of the vertices. As such, it is applicable to dry foams and unlikely to apply to the wet-foam system studied in this work. The simulations of Weaire, et al.  focus on the behavior of foams under quasi-static, extensional flow. These are applicable to wet foams and, because they are quasistatic, contain no dissipation. These simulations measure the number of T1 events. A T1 event is a topological change that involves the switching of neighbors between bubbles. Reference  reports events that include large numbers of T1 events, implying a distribution of stress drop with a power-law like behavior. In contrast, the bubble model of Durian, which is also applicable to a wet foam, predicts power-law distributions with an exponential cutoff . This model includes viscous dissipation between bubbles. Further work on this model demonstrated that the exponential cutoff increases as one approaches the melting point . However, the dependence on area fraction is weak, and true power-law behavior was never observed, even extremely close to the melting point . Finally, simulations of a q-potts model of a sheared foam  did not assume a particular form for the dissipation. This work suggests that the distribution of topological rearrangements is not power-law like; however, the distribution of energy drops may be consistent with power-law behavior .

In order to distinguish between these various theoretical models, there has been limited experimental work on two-dimensional model foams  and three-dimensional foams . The work in Ref.  made indirect measurements of the bubble rearrangements and concluded that only local rearrangements occurred during flow. The work in Ref.  used Langmuir monolayers to form two-dimensional foam. This work directly measured the distribution of T1 events and concluded that the size distribution was not consistent with a power-law. In contrast, the work in Ref.  used a single layer of bubbles trapped between a glass plate and a water surface. This work focused on the bubbles whose topological class had changed. These measurements suggest that very large events can occur.

The apparent discrepancies in these measurements may not be as severe as it first appears. Detailed simulations of various models show that the size of T1 events and the number of bubbles involved in a rearrangement are not necessarily correlated . Also, it is possible that the number of T1 events is not correlated with the size of a stress drop . This is reasonable when one considers that not all T1 events will relieve an equivalent amount of stress. Therefore, to address the question of the scaling of the size distribution, it is critical to have direct measurements of the stress or energy drops. We have accomplished this for the case of another model, two-dimensional foam: a bubble raft .

A bubble raft consists of a single layer of bubbles placed on the surface of water. Bubble rafts have been used to model the flow behavior of amorphous solids . The bubble rafts are an ideal system for the study of two-dimensional foams for a number of reasons. They allow for both direct measurements of the macroscopic properties of the foam and the bubble dynamics. There are no confining glass plates, which can add complications as the shear rate is increased. Finally, one can control both the degree of order in the foam (by varying the distribution of bubble sizes) and the density of the foam with great precision. In this paper, we will report on the flow behavior of a disordered bubble raft with $`\varphi 0.9`$.

We generate flow of the bubble raft using a two-dimensional Couette viscometer that is described in detail in Ref. . The apparatus consists of two concentric cylinders oriented vertically. Water is placed between the cylinders, and the upper surface of the water is free. The outer cylinder consists of 12 individual pieces, so it can be expanded and compressed to adjust the density of the bubble raft. The working radius of the outer cylinder was 7.43 cm. The inner cylinder is two pieces: a solid cylinder placed in the water and a second piece that is in contact with the bubble raft and fits over the solid cylinder. The second piece is hung by a torsion wire and is used to measure the stress, $`\sigma `$, on the inner rotor due to the bubble raft. It is important to note that this is a “two-dimensional stress” given by the force on the inner rotor in the tangential direction divided by the circumference of the rotor. Therefore, in terms of the torque, $`\tau `$ on the rotor and the radius of the rotor $`r`$, the stress is given by $`\sigma =\tau /(2\pi r^2)`$. The radius of the inner cylinder was 3.84 cm. A constant rate of strain is applied to the system by rotating the outer cylinder at a constant angular speed in the range 0.0005 rad/s to 0.01 rad/s

The stress on the inner cylinder is determined from the angular displacement of the torsion wire supporting the inner cylinder. The angular displacement was measured using magnetic flux. A coil was attached to the torsion wire and suspended within a high-frequency magnetic field. The induced voltage was used to determine the angle of the coil. Typical values of the angle during shear ranged from $`22^{}`$ to $`40^{}`$, corresponding to stresses in the range of 4 - 23 dyne/cm. The resolution in stress was set by our resolution in angular measurement and the torsion constant, $`\kappa =570\mathrm{dyne}\mathrm{cm}/\mathrm{rad}`$. The voltage signal was digitized using a 12-bit A to D converter in the computer. As a test of the noise level in the signal, the stress was monitored as a function of time without shear. In the absence of shear, the noise in the measured stress signal was at the level of the lowest bit in digitized signal, corresponding to changes in the stress of $`\pm 0.026\mathrm{dyne}/\mathrm{cm}`$. Therefore, when measuring the changes in the stress, changes of $`\pm 0.026\mathrm{dyne}/\mathrm{cm}`$ were ignored, providing a lower limit on the size of the stress drops.

The bubble raft is generated by flowing nitrogen through a solution of 44% by weight glycerine, 28% by weight water and 28% by weight Miracle Bubbles (Imperial Toy Corp.). The bubble size is fixed by the pressure and needle diameter. For the experiments reported on here, three different size bubbles (2 mm, 3 mm, 5 mm) were used, with approximately 500 bubbles of each size. The experiments were carried out over two hours. The coarsening of the foam without shear was monitored during this time, and no significant coarsening or stress drops due to coarsening were observed. However, after two-hours, the stability of the foam deteriorated rapidly, and shear resulted in severe rupture of bubbles.

Figure 1 is an image of a section of the bubble raft between the two cylinders. The bubble raft was monitored for slippage at both the outer and inner cylinder. For all of the experiments reported here, there was no slip between the first row of bubbles and the corresponding cylinder. Also, it was clear from the images that the stress drop events corresponded to rearrangements of the bubbles in the bulk of the foam. The average azimuthal velocity distribution was consistent with what one would expect for Couette flow. This is illustrated in Fig. 2 for an angular rotation of $`0.01\mathrm{s}^1`$. The dashed line is the theoretical velocity profile for a Newtonian fluid, and the solid curve is the theoretical velocity profile for a shear-thinning fluid with a viscosity given by $`\eta =m\dot{\gamma }^{n1}`$, where $`m`$ is a constant, $`\dot{\gamma }`$ is the shear rate, and $`n=1/3`$ . It should be noted that the bubble model predicts shear-thinning with an exponent of approximately $`1/3`$ . The fact that the velocity profile is consistent with the expected shear-thinning profile, and not the velocity profile of the water substrate, is strong evidence that the foam is acting independent of the water. Because the rate of strain varies as a function of radius, all reported values of for rates of strain are taken at the inner cylinder. Finally, it should be noted that this is very different behavior from that observed for two-dimensional foams confined between glass plates. In this case, exponential decay of the velocity profile was observed .

Figure 3 shows a typical response of the bubble raft to shear. The stress on the inner cylinder is shown as a function of the applied strain for a rate of strain of $`3.1\times 10^3\mathrm{s}^1`$. The two key features to notice are the initial elastic region where the stress increases linearly with strain and the subsequent region in which flow occurs. The intermittent nature of the flow is obvious in this plot. The yield strain was always of order 1 for all rates of strain.

Figure 4 shows the distribution of the stress drops for three different rates of strain. Here the size of the drop is normalized by the average stress per bubble. Also, as discussed above, any change in stress of $`\pm 0.026\mathrm{dyne}/\mathrm{cm}`$ was ignored and does not enter into the plotted distribution. The distribution of stress drops is not consistent with a power-law for all length scales; however, for small stress drops, the distribution is consistent with a power law with an exponent of -0.8. For comparison, the exponent predicted for the bubble model is $`0.70\pm 0.05`$ . Also, as predicted by the bubble model , the distribution is essentially independent of the rate of strain for small stress drops, and the large stress drop cutoff is weakly dependent on the rate of strain. This is difficult to see from Fig. 4, but is made clear by considering the average stress drop size. For example, for a rate of strain of $`0.031\mathrm{s}^1`$, the average stress drop (normalized by the average stress per bubble) is 27.9; whereas, it is 47.0 for a rate of strain of $`0.31\mathrm{s}^1`$.

We have conclusively demonstrated that the general features of the distribution of stress drops observed in the bubble raft are in agreement with the simulations reported in Ref.  and . Specifically, the cutoff in the distribution for large stress drops is clear. Also, the viscosity of the bubble raft is consistent with a shear-thinning fluid with an exponent of 1/3. Finally, the distribution of stress drops, particularly for the small stress-drops, is essentially independent of the rate of strain for the values studied here. This strongly suggest the existence of a quasi-static limit. Such a limit is particularly important for our geometry where the rate of strain varies across the system.

A number of important questions remain. First, what is the correspondence between the size of the stress drops and the number of bubbles that exchange neighbors? It is still possible that the larger stress drops are due to rearrangements involving a wide distribution of the number of bubbles, as reported in . Therefore, “system-wide” events are still a possibility. The current imaging system was not appropriate for careful tracking of bubbles throughout the system. Future improvements in the imaging of this system will enable this question to be addressed. Second, a detailed study of the flow behavior of bubble rafts as a function of both disorder (polydispertiy) and area fraction will provide further tests of the various two-dimensional models of foams. Of particular importance is the behavior as a function of area fraction to distinguish further between the predictions of the quasi-static simulations  and the bubble model .

The authors acknowledge funding from NSF grant CTS-0085751. M. Dennin acknowledges additional funding from the Research Corporation and Alfred P. Sloan Foundation. J. Lauridsen was supported by NSF Research Experience for Undergraduates grant PHY-9988066 and a grant from UCI-UROP program. The authors thank Mingming Wu, Andrea Liu, and Doug Durian for useful conversations.


text:
# Precise measurement of the decay 𝐾_𝐿→𝜋⁰𝛾𝛾

## 1 Introduction

The measurement of the $`K_L\pi ^0\gamma \gamma `$ decay is useful to constrain the CP conserving amplitude of the decay $`K_L\pi ^0e^+e^{}`$ via two photon exchange. Previous measurements of the decay  have been compared with calculations performed in the framework of the Chiral Perturbation Theory ($`\chi PT`$), the effective theory of the Standard Model at low energy in the hadronic sector. These predictions are best described in terms of the two amplitudes A and B (referring to angular momentum states $`J=0`$ and $`J=2`$ of the two photons respectively) in the Lorentz invariant expression for the double differential decay rate :

$$\frac{^2\mathrm{\Gamma }}{yz}=\frac{m_K}{2^9\pi ^3}\left[z^2|A+B|^2+\left(y^2y_{\mathrm{max}}^2\right)^2|B|^2\right]$$
(1)
$$z=\frac{(k_3+k_4)^2}{m_K^2}=\frac{m_{3,4}^2}{m_K^2}y=\frac{|p_K(k_3k_4)|}{m_K^2}$$

with $`m_K`$ and $`p_K`$ the kaon mass and momentum, $`k_3`$ and $`k_4`$ the two photon momenta and $`y_{\mathrm{max}}`$ the kinematic bound for the $`y`$ variable, given by:

$$y_{\mathrm{max}}=\frac{1}{2}\sqrt{\left(1\frac{m_{\pi }^{}{}_{}{}^{2}}{m_{K}^{}{}_{}{}^{2}}\right)^2+\left(1+\frac{m_{\pi }^{}{}_{}{}^{2}}{m_{K}^{}{}_{}{}^{2}}\right)^2\times z+z^2}$$

The leading $`O(p^4)`$ $`\chi PT`$ calculation predicts B=0, in qualitative agreement with the experimental observation of a $`z`$ spectrum peaked at high values, but underestimates the $`K_L\pi ^0\gamma \gamma `$ branching ratio by about a factor of three . At $`O(p^6)`$ the rate and the $`m_{\gamma \gamma }`$ spectrum can be reproduced by adding a contribution from the VMD mechanism , via the coupling constant $`a_v`$ . The VMD mechanism could enhance the state $`J=2`$ for the two photons, hence allowing a sizeable CP conserving contribution which is not helicity suppressed for the $`K_L\pi ^0e^+e^{}`$ decay .
In this paper we present a new measurement of the branching ratio and decay spectrum of the decay $`K_L\pi ^0\gamma \gamma `$, based on data collected in the years 1998 and 1999 with the apparatus of the NA48 experiment. Quantitative information on the VMD contribution is obtained from the spectrum at low z values.

## 2 Experimental set-up

The $`K_L\pi ^0\gamma \gamma `$ events are collected using the NA48 detector located at the CERN-SPS and primarily devoted to measure the parameter $`ϵ^{}/ϵ`$ which characterises direct CP violation in neutral Kaon decays into two pions . Neutral Kaons are produced in interactions of $`1.4\times 10^{12}`$ 450 GeV/$`c`$ protons on a beryllium target during 2.4 seconds every 14.4 seconds. A system of sweeping magnets and collimators defines a neutral beam of $`2\times 10^7K_L`$ per burst and divergence of $`\pm 0.15`$ mrad. A fraction of the non interacting protons are redirected by channelling in a bent crystal, to a second target, 120 metres downstream of the first, to generate a $`K_S`$ beam . In order to tag $`K_S`$ decays, the protons before hitting this target are detected by an array of scintillation counters (tagger). The two beams converge with an angle of $`0.6`$ mrad at the calorimeter, $``$ 120 metres downstream of the final collimators. The decay region is contained in an evacuated cylindrical vessel 89 metres long separated from the NA48 detector volume by a thin Kevlar window $`0.3\%`$ of a radiation length ($`X_0`$). The neutral beam traverses the detector inside a 16 cm diameter vacuum pipe. A scintillation counter (AKS) placed on the $`K_S`$ beam defines the upstream edge of the fiducial region for decays from that beam.
The most important detector element for this analysis is the quasi-homogeneous liquid Krypton electro-magnetic calorimeter (LKr) structured in 13212 readout tower cells $`27X_0`$ deep . The ionization signal from each of the cells is integrated, amplified, shaped and digitised by 10-bit FADCs at 40 MHz sampling frequency. The energy resolution can be parameterised as:

$$\sigma (E)/E(0.09\pm 0.01)/E(0.032\pm 0.002)/\sqrt{E}(0.0042\pm 0.0005),$$

with E in GeV. The spatial and time resolutions are better than $`1.3`$ mm and $`300`$ ps, respectively, for photons with energy greater than 20 GeV.
Seven ring shaped counters (AKL), consisting of plastic scintillator and iron converters, surround the decay and spectrometer region in order to veto events with photons outside the calorimeter acceptance. Their efficiency has been estimated to be around $`95\%`$. The charged particles are reconstructed by a magnetic spectrometer consisting of four drift chambers (DCH) and a magnetic dipole. The space resolution for each projection is $`90\mu `$m and the average efficiency is better than $`99.5\%`$ per plane. For this analysis the spectrometer is used to veto events with reconstructed tracks.

## 3 Data taking

Events of the $`K_L\pi ^0\gamma \gamma `$ decay channel and of the channel $`K_L\pi ^0\pi ^0`$, used for normalisation, are collected by the same trigger, since their final states appear identical in all detector elements. The trigger decision is based on quantities which are derived from the sums of the energy deposited in the LKr calorimeter in groups of cells corresponding to 64 horizontal and 64 vertical slices . The trigger requirements are:

* at maximum 5 peaks in each projection;
* total energy larger than 50 GeV;
* centre of gravity of the event, computed from the first moments of the energy peaks in the projections to be within 15 cm from the beam axis;
* proper time of the decay, computed from the second moments, less than $`5\tau _{K_S}`$ from the beginning of the decay region.

The trigger efficiency has been checked to be better than $`99.9\%`$  using a minimum bias sample triggered by a scintillating fibre hodoscope placed in the LKr volume at a depth of $`9.5X_0`$.
During the experimental runs roughly 180 Terabytes of raw data were recorded and pre-selected by an on-line software filter, whose criteria are subsequently tightened in the off-line analysis.

## 4 Analysis

The first step is to reconstruct the showers in the liquid krypton calorimeter summing the energies of the cells in a circle of $`11`$ cm radius from a local maximum. The time and the position of the shower are derived from the most energetic cells and from the centre of gravity of the energy. A partial overlap between two showers is resolved using the shower profile obtained from simulation. The shower energy is corrected for:

* the residual 0.4% non-uniformity in the calibration electronics chain, comparing energy and momentum of electrons from $`K_{e_3}`$ events as a function of the hit cell. The uniformity is improved to the 0.15% level;
* the 1% peak-to-peak modulation of the response as a function of the impact point in the LKr cell due to the effects of a $`2`$ mm gap between electrodes in the vertical direction and to the finite integration time in the horizontal one;
* the energy-leakage in the beam pipe at small radius, the material before the active volume and the presence of not properly working cells;
* the small (0.1%) non-linearity.

Since the $`K_L\pi ^0\gamma \gamma `$ event is characterised by four clusters with the only constraint of two showers coming from a $`\pi ^0`$ decay, all the events with at least four clusters are considered. Each cluster must be located well inside the calorimeter volume, i.e. within an octagon of 113 cm apothem, with 15 cm inner radius. It must be at a distance larger than 2 cm from any non working cell and separated from another cluster by at least $`10`$ cm. The cluster energy must be between 3 GeV and 100 GeV and the time within 3 ns of the event time, defined as the average of the cluster times. The total energy of the clusters must be in the 70-170 GeV interval.
The two photon invariant mass is:

$$m_{\gamma _i\gamma _j}=\frac{\sqrt{E_iE_jd_{i,j}^2}}{z_{LKr}z_K}$$
(2)
where $`E_i,E_j`$ are the cluster energies, $`d_{i,j}^2=[(x_ix_j)^2+(y_iy_j)^2]`$ with $`x_i`$ and $`y_i`$ the horizontal and vertical positions of the i-th cluster at the LKr front face. The longitudinal coordinate of the decay vertex, $`z_K`$, is computed from the electro-magnetic calorimeter information assuming the nominal kaon mass for the parent particle:

$$z_K=z_{LKr}\frac{\sqrt{_{i,j,i>j}E_iE_jd_{i,j}^2}}{m_K},$$
(3)
where $`z_{LKr}`$ is the longitudinal coordinate of the LKr front face. The overall energy scale is fixed by the fit of the position of AKS counter, that vetoes all $`K_S`$ decays occurring upstream .
We only consider events for which there is no proton detected in the $`K_S`$ tagger in time (within $`\pm 2`$ ns) with the event; this rejects $`10.65\%`$  of the $`K_L`$ decays, due to accidental coincidences with a proton in the tagger.
Because of the small branching ratio, the $`K_L\pi ^0\gamma \gamma `$ signal is subject to large backgrounds such as $`K_L\pi ^0\pi ^0\pi ^0`$ decays, with missing and/or overlapping photons, mis-reconstructed $`K_L\pi ^0\pi ^0`$ decays with $`\gamma `$ conversion in the spectrometer and events generated by accidental pile up of particles from two different events. We distinguish two classes of cuts in the analysis: $`2\pi ^0`$ rejection and $`3\pi ^0`$ rejection.

### 4.1 $`2\pi ^0`$ rejection.

Three different combinations of $`\gamma \gamma `$ pairs can be formed from four showers. A $`\chi ^2`$ variable is defined:

$$\chi ^2=\left[\frac{(m_{1,2}+m_{3,4})/2m_{\pi ^0}}{\sigma _+}\right]^2+\left[\frac{(m_{1,2}m_{3,4})/2}{\sigma _{}}\right]^2$$
(4)
where $`\sigma _\pm `$ are the resolutions of $`(m_{1,2}\pm m_{3,4})/2`$ measured from the data and parametrised as a function of the lowest photon energy, typically $`\sigma _+0.42`$ MeV/$`c^2`$ and $`\sigma _{}0.83`$ MeV/$`c^2`$. In order to select a $`\pi ^0\pi ^0`$ candidate the best combination is required to have $`\chi ^213.5`$, while the $`\pi ^0\gamma \gamma `$ events are selected with $`\chi ^2>300`$. We retain for further analysis all combinations that give a photon pair (that we call $`m_{1,2}`$) within 3 MeV/$`c^2`$ from $`m_{\pi ^0}`$ and the other pair (corresponding to $`m_{3,4}`$) outside the window $`110160`$ MeV/$`c^2`$.
The $`2\pi ^0`$ background concentrates in the most interesting region for the determination of $`a_v`$, namely at low $`m_{3,4}`$ value ($`m_{3,4}<240`$ MeV/$`c^2`$) where the $`K_L\pi ^0\gamma \gamma `$ signal is smaller. Due to mis-reconstructions one of the $`\pi ^0`$ may be reconstructed with a wrong mass. This can happen in case of a $`\gamma `$-conversion or $`\pi ^0`$ Dalitz decay with one electron outside the calorimeter acceptance that biases the photon energy reconstruction. We use the first three drift chambers as a veto against conversions and Dalitz decays if there are hits in coincidence, $`\pm 20`$ ns, with the event time. In addition we reject events with an overflow condition in the readout of more than one plane of the drift chambers. This condition is defined by the occurrence of more than seven hits in a DCH plane in any 100 ns time interval. In this case the front-end readout buffers are reset and the DCH information cannot be used reliably.
The remaining $`2\pi ^0`$ background is reduced requiring $`p_{T,4}>\mathrm{\hspace{0.17em}40}`$ MeV/$`c`$, where $`p_{T,4}`$ is the transverse momentum to the kaon direction of $`\gamma _4`$, which is the lowest-energy photon of the non-$`\pi ^0`$ pair of photons. The variable $`p_{T,4}`$ is on average smaller for $`2\pi ^0`$ events than for $`\pi ^0\gamma \gamma `$ events and even lower when the photon energy is underestimated. Although this requirement reduces the acceptance for $`\pi ^0\gamma \gamma `$ at high $`y`$ and low $`m_{3,4}`$, it does not degrade the sensitivity to the $`B`$ amplitude as the $`\left(y^2y_{\mathrm{max}}^2\right)`$ factor in equation 1 is already small.

### 4.2 $`3\pi ^0`$ rejection.

In order to reduce the background from $`K_L\pi ^0\pi ^0\pi ^0`$ decays, we require the event to have only four clusters in-time, within 3 ns, and no activity in the AKL veto counters in coincidence with the event. As most of the $`3\pi ^0`$ background events have one or two photons outside the LKr acceptance, the reconstructed $`z_K`$ according to formula (3) moves downstream the beam axis by at least 6 metres. Therefore, the fiducial volume is restricted to the first $`30`$ metres of the $`K_L`$ decay region. Events with missing energy are also reduced by requiring the center of gravity, computed from the first momenta of the clusters ($`R_{COG}`$), to be less than $`4`$ cm from the beam axis.
The simulation of $`3\times 10^9K_L\pi ^0\pi ^0\pi ^0`$ events shows that after applying these cuts the background is composed of events with missed photons: $`36\%`$ of the sample has two missing photons, while $`61\%`$ of the sample results from one photon missing the detector and two showers overlapping in the LKr. The remaining events have all the six photons inside the calorimeter but only four reconstructed because of two overlaps. The corresponding showers are broader than for a single isolated photon.
Events with missed photons are rejected on the basis of a variable, $`Z_{max}`$, designed to estimate, for these events, the true kaon decay vertex position. The same estimate, when applied to signal events, produces an unphysical result, with a distribution showing many events in the region upstream of the final collimator. The distribution of this variable is shown in figure 1. A condition $`Z_{max}5`$ m rejects almost $`99\%`$ of the $`3\pi ^0`$ background while retaining $`46\%`$ of the signal.
The procedure to define $`Z_{max}`$ is based on the analysis of the possible topologies of events with missed photons, as shown in the following figure:

* $`\underset{¯}{\mathrm{𝐜𝐚𝐬𝐞𝐬}\mathrm{𝟏},\mathrm{𝟐}\mathrm{𝐚𝐧𝐝}\mathrm{𝟑}:}`$
a $`\pi ^0`$ can be reconstructed and the true vertex coordinate $`z_{\gamma \gamma }`$ is obtained using equation (2) with $`m_{\gamma _i\gamma _j}=m_{\pi ^0}`$. There are five possible $`\gamma \gamma `$ combinations, other than the one of the reconstructed $`\pi ^0`$, but we consider only those that have $`z_{\gamma \gamma }`$ smaller than ($`z_K6`$ m). We call $`Z_\pi `$ the vertex estimate which is closest to this limit.
* $`\underset{¯}{\mathrm{𝐜𝐚𝐬𝐞}\mathrm{𝟒}:}`$
one $`\pi ^0`$ has a missing photon and two photons from the other $`\pi ^0`$’s overlap. Starting from the hypothesys of two $`\pi ^0`$s which share a photon, we compute, for each group of three photons, the vertex coordinate $`z_{\gamma \gamma \gamma }`$, assuming each photon in turn to be the superposition of two photons with the same position in the calorimeter. For example if photons 1, 2 and 3 come from two $`\pi ^0`$s and photon 2 is overlapped, then $`z_{\gamma \gamma \gamma }`$ is computed as:

$$z_{\gamma \gamma \gamma }=z_{LKr}\frac{1}{m_{\pi ^0}}\left(\frac{1}{E_1E_2d_{1,2}^2}+\frac{1}{E_2E_3d_{2,3}^2}\right)^{\frac{1}{2}}$$
(5)
Out of the twelve possible combinations we consider only the triplets with an invariant mass below $`(m_Km_{\pi ^0})`$ and the vertex estimate closest to, but smaller than, $`(z_K6m)`$.

For each event we define the variable $`Z_{max}=max(Z_{\gamma \gamma },Z_{\gamma \gamma \gamma })`$.

The residual $`3\pi ^0`$ background originates from two overlaps of photons such that the $`z_K`$ is reconstructed unbiased, but the shape of one or two showers is different from that of a single isolated photon. As estimator of the cluster width we use the RMS of its energy profile in $`x`$ and $`y`$:

$$w_{x(y)}=\sqrt{\frac{_iE_i\mathrm{\Delta }x(y)_i^2}{_iE_i}}i5\times 5\mathrm{cells}\mathrm{around}\mathrm{the}\mathrm{impact}\mathrm{cell}$$
(6)
where $`\mathrm{\Delta }x(y)_i`$ is the distance in the $`x(y)`$ direction of the $`i`$-th cell from the center of gravity of the cluster energy. An energy-dependent cut is applied to the width of each cluster. The value of the cut is defined in terms of the average $`\overline{w}`$ and the standard deviation $`\sigma `$ of the distribution of $`w`$ measured in photons from good $`2\pi ^0`$. We require $`w<(\overline{w}+1.8\sigma )`$ on the photons of $`m_{3,4}`$ and $`w<(\overline{w}+3\sigma )`$ on the photons of $`m_{1,2}`$. The signal reduction is at the level of $`15\%`$.

## 5 Background evaluation

The $`K_L\pi ^0\gamma \gamma `$ candidates that survive all the above selections are 2558 in the signal region defined as $`132<m_{1,2}<138`$ MeV/$`c^2`$. For 345 of these events more than one combination satisfies the $`m_{1,2}`$ cut and therefore the mass assignment is ambiguous.
We estimate the contamination from $`2\pi ^0`$ decays and pile-up processes directly from the data. To quantify the $`K_L\pi ^0\pi ^0\pi ^0`$ background we rely on Monte Carlo.
We study the $`K_L\pi ^0\pi ^0`$ background in the signal region using $`K_S\pi ^0\pi ^0`$ decays from different $`K_S`$ data samples: the proton tagged $`K_S`$ sample, collected concurrently with the $`K_L`$ data, and a special $`K_S`$ high intensity run ($`K_SHI`$) taken in 1999. In the case of the $`K_S`$ tagged sample we ask for a tight coincidence within $`\pm 0.5`$ ns between the tagger and the LKr signals. In this way we enrich the number of $`2\pi ^0`$ events contributing to the background. Conversely the only contribution to $`\pi ^0\gamma \gamma `$ in this sample comes from the $`K_L`$ contamination due to accidental coincidences in the tagger, which is estimated to be $`3.5\%`$ with the cut used.

After applying the $`\pi ^0\gamma \gamma `$ selection criteria, the $`m_{3,4}`$ invariant mass distribution is shown in figure 2 for events in the tagged sample compared to the original $`K_L\pi ^0\gamma \gamma `$ mass distribution scaled by the $`K_L`$ contamination factor. The $`2\pi ^0`$ background is estimated from the tiny excess of tagged events measured where the $`\pi ^0\gamma \gamma `$ signal is smaller, namely for $`m_{3,4}<0.240`$ MeV/$`c^2`$ while for $`m_{3,4}>0.240`$ MeV/$`c^2`$ it is linearly estrapolated from control regions on the sidebands of the $`\pi ^0`$ peak ( $`127<m_{1,2}<130`$ MeV/$`c^2`$ or $`140<m_{1,2}<143`$ MeV/$`c^2`$). This procedure leads to a background evaluation of $`(4.1\pm 1.6_{(stat)}\pm 1.4_{(syst)})`$ events. Here the systematic error is related to the subtraction of the $`K_L`$ contamination and to the fact that $`K_S`$ events have a distribution of the decay vertex populated only in the region between 6 and 18 metres.
The high intensity $`K_S`$ data sample provides about $`60`$ % of the tagged $`K_S\pi ^0\pi ^0`$ statistic with negligible $`K_L`$ component. The number of background events found is in agreement with the previous estimate.
The cuts applied to reduce the $`K_L\pi ^0\pi ^0\pi ^0`$ background are mainly of geometrical nature and can safely be simulated. The most critical cut is the one on the cluster width, for which an accurate description of the LKr response for overlapping showers is needed. For this last variable we compare data and Monte Carlo distributions in figure 3.

The simulation describes the data within $`3\%`$ both in shape and quantity being normalised to the Kaon flux. After the cut the estimated residual background is $`(70.2\pm 8.2_{(stat)}\pm 7.0_{(syst)})`$ events, where the systematic uncertainty comes from the non-Gaussian tail of the calorimeter response mainly due to photo-hadron reactions.
Finally, the pile-up events, which could originate for example from two overlapped in-time kaon decays, can be quantified comparing the normalised distributions of the center of gravity at the LKr for $`K_L\pi ^0\gamma \gamma `$ candidates and for a good $`2\pi ^0`$ data sample. An excess at $`R_{COG}>4`$ cm is found in the $`K_L\pi ^0\gamma \gamma `$ distribution, after having subtracted the $`3\pi ^0`$ background (figure 4). When the excess is linearly extrapolated under the $`K_L\pi ^0\gamma \gamma `$ signal we obtain a pile-up background level of $`(8.1\pm 3.5_{(stat)}\pm 4.0_{(syst)})`$ events. The choice of extrapolation function is dominating the systematic error; however the linear one is expected from a uniform center of gravity distribution and is confirmed by the $`R_{COG}`$ distribution for events with out-of-time clusters.

In Table 1 is reported the summary of the background evaluation in the signal region and, as cross-check, in the control regions defined above.

The overall background is $`3.2\%`$. The invariant mass distributions $`m_{1,2}`$ and $`m_{3,4}`$ are shown in figure 5.

The background estimate is in agreement with that found in the two control regions where 44 background events are observed in the data with an expectation of $`38\pm 7`$.

## 6 Results and discussion

As confirmed by the Monte Carlo, most of the events in the low $`m_{3,4}`$ mass region shown in figure 5 are genuine $`K_L\pi ^0\gamma \gamma `$ events which suffer from the wrong combination of photons. Events with an ambiguous assignment for $`m_{1,2}`$ have a flat distribution in the $`m_{3,4}`$ mass region which will inflate the amount of events in the low-mass $`m_{3,4}<0.240`$ Mev/$`c^2`$ region. To extract a value for the vector coupling constant $`a_v`$ we exclude events with ambiguous mass assignment. The events in the low-mass region are thus reduced by $`50\%`$. In Table 2 we report the numbers of $`\pi ^0\gamma \gamma `$ candidate, background and acceptance in the 20 bins of the $`[0400]`$ MeV$`/c^2`$ $`m_{3,4}`$ interval.

We perform a fit of the bidimensional distribution of the two relevant kinematic variables $`m_{3,4}`$ and $`y=|E_3E_4|/m_K`$, defining a likelihood function as follows:

$$lnL=\underset{i,N_{bins}}{}[N_i(y,m_{\gamma \gamma })lnE_i(y,m_{\gamma \gamma })E_i(y,m_{\gamma \gamma })]$$
(7)
where $`N_i(y,m_{\gamma \gamma })`$ is the measured number of events in the $`i`$-th bin while $`E_i(y,m_{\gamma \gamma })`$ is the corresponding expectation value given by the simulated $`K_L\pi ^0\gamma \gamma `$ plus the estimated background. The model adopted for the simulation of the $`K_L\pi ^0\gamma \gamma `$ process is described in  and it is normalised to the total number of $`K_L`$ decays in the fiducial volume. The uncertainty of the parameterisation of the $`K_L3\pi `$ vertex  is included in the systematic error. The best fit, with a $`\chi ^2`$ of $`31.1`$ for $`30`$ degrees of freedom, gives:

$$a_v=0.46\pm 0.03_{(stat)}\pm 0.04_{(syst)}$$

The systematic error includes the following contributions: 0.02 from the uncertainty in the parameters of the theoretical model, 0.03 from the uncertainty of the acceptance evaluation and 0.02 from background evaluation and subtraction.
In figure 6 we show the $`m_{1,2}`$ distribution together with the evaluated background and the $`O(p^6)`$ $`\chi PT`$ expectation for the fitted value of the vector meson exchange parameter $`a_v`$ in the three low-mass regions a) $`m_{3,4}[30110]`$MeV/$`c^2`$, b) $`m_{3,4}[160240]`$MeV/$`c^2`$ and c) $`m_{3,4}[240260]`$MeV/$`c^2`$. Our data clearly disfavour $`a_v=0`$.

The branching ratio is obtained from the number of $`\pi ^0\gamma \gamma `$ candidates, including those with ambiguous mass assignment, normalising to the number of $`\pi ^0\pi ^0`$ events observed in the same sample of $`K_L`$ decays and selected by the same cuts, in order to minimise the uncertainty in the acceptance corrections and is found to be:

$$BR(K_L\pi ^0\gamma \gamma )=(1.36\pm 0.03_{(stat)}\pm 0.03_{(syst)}\pm 0.03_{(norm)})\times 10^6$$

Both the correction due to the Dalitz decay difference between the $`2\pi ^0`$ normalisation mode and the $`\pi ^0\gamma \gamma `$ decay and the one due to the $`K_LK_S`$ interference in $`\pi ^0\pi ^0`$ decays  are taken into account. The uncertainty related to the experimental knowledge of the $`K_L\pi ^0\pi ^0`$ branching ratio  is quoted separately. The computed average value of the acceptance for the $`\pi ^0\gamma \gamma `$ process is $`7.62\%`$ and for the $`\pi ^0\pi ^0`$ is $`8.72\%`$. The uncertainty related to this computation is the main source of the $`2.5\%`$ systematic error in the branching ratio measurement since, as shown in Table 2, it has sharp variations. The residual background estimation, the uncertainty in $`a_v`$ and the exclusion of the ambiguous events each gives less than 1% variation on the branching ratio value.
As shown in figure 6 there is no evidence for a signal in the invariant mass region $`m_{3,4}<m_{\pi ^0}`$ where we count 9 $`\pi ^0\gamma \gamma `$ candidates and ($`4.5\pm 2.7`$) background events. Using events with $`y<0.2`$ in 5 bins of $`m_{3,4}`$ between $`30`$ and $`110`$ MeV/$`c^2`$, we compute a model-independent upper limit :

$$BR(K_L\pi ^0\gamma \gamma )|_{m_{3,4}[30110]Mev/c^2,y[00.2]}<0.6\times 10^890\%C.L.$$

For this calculation we use the acceptance computed from events generated according to phase space, in a restricted region where it is almost uniform. Given the negligible contribution of the amplitude $`A(J=0)`$, our measurement of the vector coupling constant $`a_v`$ allows the CP-conserving component of the $`K_L\pi ^0e^+e^{}`$ decay  to be predicted as:

$$BR(K_L\pi ^0e^+e^{})|_{CPconserving}=(4.7_{1.8}^{+2.2})\times 10^{13}$$

This suggests that the CP-violation effects dominate in the $`K_L\pi ^0e^+e^{}`$ mechanism.

## ACKNOWLEDGEMENT

It is a pleasure to thank the technical staff of the participating laboratories, universities and affiliated computing centres for their efforts in the construction of the NA48 apparatus, in the operation of the experiment, and in the processing of data. We also grateful to G.D’Ambrosio and G.Valencia for useful discussions.


text:
# GMRT HI imaging survey of ERIDANUS group of galaxies

## 1 Introduction

Eridanus is a nearby group of galaxies at a distance of $``$20 Mpc. The group members are arranged in a filament like structure joining Fornax and Dorado clusters of galaxies. Total number of optically identified galaxies in this group is $`100`$. This group has almost equal number of early type and late type systems. The redshift is known for all galaxies through optical spectroscopy. The velocity dispersion of this group is $`350`$ km s<sup>-1</sup>. The dynamical study of this group shows that the group is a bound system comprising of several sub-groups . These sub-groups are merging together and will eventually evolve to a massive cluster.

## 2 Motivations for the present study

The GMRT data of HI emission and radio continuum from the Eridanus galaxies, combined with the existing optical and infrared data will be used to study –

1. Rotation curves of the disk galaxies
2. Tully-Fisher relation using rotation curve and 2MASS magnitudes
3. M/L ratios and dark matter haloes of the disk galaxies
4. Star formation rates in the disk galaxies
5. Radio continuum morphology and spectral indices of galaxies

## 3 Results

We have observed 26 pointings, each with a field of view of $``$24’, using the GMRT and most of the data have been analysed. Further observations are proposed with the GMRT. The HI image and the rotation curve of NGC 1385, a member of this group, are shown in figure 1.


text:
# A simple model for the relative stabilities of DNA hairpin structure

## Abstract

A model of self-avoiding walk with suitable constraints on self-attraction is developed to describe the conformational behaviour of a single stranded short DNA molecule that form hairpin structure. Using exact enumeration method we calculate the properties associated with coil-hairpin transition. Our results are in qualitative agreement with the experiment.

The structural and dynamic behaviour of a single stranded DNA molecule that form a stem-and-loop (hairpin) structures in solutions is a subject of current interest 1 ; 2 ; 3 ; 4 . This is because DNA hairpin structure is known to participate in many biological functions, such as the regulation of gene expression 5 , DNA recombination 6  and facilitation of mutogenic events 7 . The DNA hairpin-containing domains may provide potential binding sites for exogenous drugs and endogenous proteins 8 . It can also be used as DNA biosenser (e.g. molecular beacons) 9 .

The single stranded DNA that form a stem-and-loop structures under a given set of solution conditions are formed from synthetic oligonucleotides. For example, the oligonucleotide $`5^{}`$CCCAA \- $`(X)_m`$ \- TTGGG$`3^{}`$, where $`X`$ is any one of the nucleotides and $`m`$ its number, is such a molecule which can easily be synthesized. Due to thermal fluctuations such a molecule may acquire different conformations. In a simplified description all of the configurations can be divided into two main states; the open and the closed one as shown in Fig. 1. The closed hairpin structure is stabilized due to pairing of complimentary bases at the two arms of the molecule when hydrogen bonds are formed between them. The open state has high entropy due to the large number of configurations achieved by a single stranded DNA chain. The closed-to-open transition requires sufficiently large energy to unzip all the base pairs, whereas the closing transition requires the two arms of hairpin to come close to each other in space so that hydrogen bonding between the complimentary nucleotides can take place. Experimentally one finds the rate of closing depends strongly on the properties of the hairpin loop, such as the length and rigidity, whereas the rate of opening is relatively unaffected by these properties 1 .

In this note we describe a model and show that it explains some of the conformational behaviour associated with the formation of hairpin structure.

We represent a linear polymer chain of $`N`$ monomers by a self-avoiding walk of $`(N1)`$ steps on a two dimensional square lattice. Visited lattice sites represent the monomers of the chain. The first $`n`$ sites of the walk represent nucleotides, say, A or C and the last $`n`$ sites the complimentary nucleotides, T or G. The remaining $`(N2n)`$ = $`m`$ sites represent the nucleotides of one type (either A or T or C or G). The repulsion at short distance between monomers (i.e. excluded volume) interaction is taken into account by the condition of self avoidance. When the complimentary monomers in a walk occupy nearest sites in a particular direction (say, clock wise), they form a pair (known as base pair). The energy associated with this pair is $`ϵ`$ ($`ϵ<0`$). The base paring cannot take place if the nucleotides approach to each other in any other way. This condition has been imposed to take into account the fact that in a single strand DNA chain the pairing between complementary nucleotides takes place only when they approach each other directly without the sugar phosphate strand coming in between as shown in Fig. 2a. Fig. 2b shows the situation in which the sugar phosphate strand is in between and therefore pairing can not take place.

All possible conformations of walks of N sites mapped by a self- avoiding walk are generated using exact enumeration technique 10 ; 11 . Since the time involved in enumerating these conformations increases as $`\mu ^N`$ where $`\mu `$ is the connectivity constant of the lattice, we have restricted ourselves to $`12N22`$.

The partition function of our interest is

$$Z_N=\underset{p=0}{\overset{n}{}}C_N(p)(e^{ϵ/k_BT})^p$$
(1)
where $`p`$ is the number of base pairs and $`C_N(p)`$ the total number of configurations corresponding to walk of $`(N1)`$ steps with $`p`$ number of pairs (base pairs). The Helmholtz free energy of the system is given as

$$F=k_BT\mathrm{ln}Z_N(T)$$
(2)
where $`k_B`$ is the Boltzmann constant and $`T`$ the temperature of the system. From Eq. (2) one can calculate the specific heat and entropy using the following relations:

$`C_v`$ $`=`$ $`T\left({\displaystyle \frac{^2F}{T^2}}\right)`$ (3a)
$`S`$ $`=`$ $`\left({\displaystyle \frac{F}{T}}\right)`$ (3b)
In our calculation we have taken $`n=5`$ and $`ϵ=0.08`$ eV. Dependence of the free energy ($`F`$) on temperature for $`N=22`$ is shown in Fig. 3a.

From this figure it is clear that there is a “discontinuous” change in the free energy at $`T=338^o`$K. In Fig. 3b we plot the specific heat as a function of temperature for several values of $`N`$. This curve has a peak suggesting a transition from the hairpin structure to the coil (or open) state. For $`N=22`$, the peak is at $`337\pm 2^o`$K. Since we are describing a finite system, the transition is not sharp; it takes place over an interval of temperature. This interval is obvious from the width of the peak in $`C_v`$ and also from Fig. 3c in which we plot entropy ($`S`$) as a function of $`T`$ for $`N=22`$. The entropy is low at low temperature as the molecule is in a hairpin structural form whereas at high temperature the entropy is higher as the molecule in coil form assesses a large number of conformations.

The partition function defined by Eq. (1) has six terms corresponding to six values of $`p`$ from 0 to 5 ($`p=0`$ corresponds to open state of the chain whereas $`1p4`$ corresponds to partially closed state and $`p=5`$ the final state of hairpin structure. We can find the probability of each of these states in the following way.

$$P_p(T)=\frac{Z_p}{Z}$$
(4)
where $`Z=_{p=0}^5Z_p`$

We compare $`P_p(T)`$ for different values of $`p`$ in Fig. 4a. It shows that the dominating structure are either all closed or all open.

The other quantity of interest is the fraction of open to closed structures at a given temperature. The fraction of all open structure can be derived from the relation

$$\chi _{}=\frac{Z_0}{Z_0+Z_c}$$
(5)
where $`Z_0`$ is for $`p=0`$ and $`Z_c=_{p=1}^5Z_p`$ and the fraction of closed structures

$$\chi _+=\frac{Z_c}{Z_0+Z_c}$$
(6)
From Eq. (5) and (6) we get

$$\frac{\chi _{}}{\chi _+}=\mathrm{exp}\left[\frac{1}{k_BT}(F_0F_c)\right]$$
(7)
where $`F_0=k_BT\mathrm{ln}Z_0`$ and $`F_c=k_BT\mathrm{ln}Z_c`$

The fraction of all closed structure at a given temperature is found from the relation

$$\chi _+^{}=\frac{Z_5}{Z_0+Z_c}$$
(8)
From this figure \[Fig. 4b\] one may note that while $`\chi _{}`$ and $`\chi _+^{}`$ intersects at $`T=337\pm 2^o`$K in agreement with the coil-hairpin transition temperature, $`\chi _{}`$ and $`\chi _+`$ intersects at $`T=366\pm 2^o`$K. As shown in Fig. 4a the probability of partially connected conformations are small and confined to narrow temperature region. One may therefore conclude that the conformational properties of this kind of single DNA strand is primarily controlled by either all open or all closed states.

By varying the value of $`ϵ`$, one can study the effect of base pair interaction on coil-hairpin transition. The change in $`ϵ`$ value from -0.08 eV to -0.05 eV leads to the decrease in transition temperature by $`100^o`$K.

In conclusion, we have studied a simple model for DNA hairpin loop structure and their relative stabilities. We use lattice models of polymers along with exact enumeration technique to study the stability of different conformations of DNA hairpin loop and the corresponding hairpin-coil transition. Our results are in qualitatively agreement with the reported experimental work by Bonnet et al 1 .

###### Acknowledgements.

We thank Navin Singh for many helpful discussions. Financial assistance from $`INSA`$, New Delhi and DST, New Delhi are acknowledged.


text:
# The Axially Symmetric Ejecta of Supernova 1987A

## 1 Introduction

SN 1987A is a Type II supernova in which the central iron core of a massive star collapsed to form a dense neutron star or black hole. Most of the conventional models assume that the explosion process is roughly spherical. The model assumption of sphericity can be tested by observations. It has been shown that core-collapse supernovae are normally linearly polarized at a level around 1% (Trammell et al. 1993; Wang et al. 1996; Tran et al. 1997; Wang et al. 2001; Leonard et al. 2000, 2001; Leonard & Filippenko 2001). For a given supernova, the degree of polarization evolves with time in a way that the closer we see into the center the higher the degree of polarization (Wang et al. 1996; Wang et al. 2001; Leonard et al. 2001). The degree of polarization seems also to be correlated with the mass of the stellar envelope; more massive envelopes show a lower degree of polarization (Wang et al. 2001; Leonard & Filippenko 2001). The linear polarization measures the asphericity of the supernova photosphere, and the polarization data suggest models in which the central engine that drives a core-collapse supernova explosion is highly aspherical (Höflich, Wheeler, & Wang 1999; Khokhlov & Höflich 2000; Wheeler, et al. 2000; Höflich, Khohklov & Wang 2001).

As the brightest and the best observed supernova of modern astronomy, SN 1987A provides a unique test case of the geometry of a core-collapse supernova. Recent Hubble Space Telecsope images and spectroscopy are able to resolve the ejecta and it is now possible to re-visit SN 1987A and see if the early and late observations can be integrated into a single self-consistent picture. In the following, we will discuss the observed structures around SN 1987A from the outermost circumstellar material inward to the central rapidly expanding ejecta and provide a scenerio for many of the features in terms of a jet-induced explosion.

## 2 The Circumstellar Material of SN 1987A

The circumstellar material (CSM) of SN 1987A is highly structured (Wampler et al. 1990; Crotts & Heathcote 1991). The nebulosity consists of a diffuse bow-shock shaped structure called “Napolean’s Hat” and several nebular rings (Wampler et al. 1990; Wang & Wampler 1992). The nebular rings were shown beautifully by images taken by the HST (Plait et al. 1995; Burrows et al. 1995; Panagia et al. 1996). The bright inner nebular rings were found to be enriched in helium and nitrogen and are believed to be material lost by the progenitor star before the explosion (Wang 1991; Wang & Mazzali 1992; Lundqvist & Fransson 1996). Napolean’s Hat is also likely to be associated with the mass loss of the progenitor star (Wang & Wampler 1992). The innermost bright (CSM) ring observed by the HST reveals a symmetry axis of the nebulosity at a position angle (PA) of 89o that is roughly consistent with the symmetry axis of the Napolean’s Hat nebula (Burrows et al. 1995). Napolean’s Hat and the circumstellar rings were formed long before the explosion of SN 1987A, but they define a geometry that extends deep into the core of the supernova.

## 3 The Supernova Ejecta

Interior to the smallest circumstellar ring are the ejecta of the supernova, the material of which is expanding at velocities from about 40,000 km s-1 in the outermost hydrogen envelope (Pun et al. 1995) down to about 1,300 km s-1 in the carbon/oxygen-rich core. The ejected material was too compact at early times for normal direct imaging, but the asymmetry of the ejecta was revealed by speckle observations (Papaliolios et al. 1989) beginning in June 1987, after the second optical maximum. Extensive linear spectropolarimetry measurements were obtained during the first year after the supernova explosion (e.g. Bailey 1988; Mendez et al. 1988; Schwartz 1987; Jeffery 1991). The linear polarization is produced by electron scattering (Shapiro & Sutherland 1982; Höflich 1991) and probes the geometry of the photosphere (see, however, Wang & Wheeler 1996 for an alternative model involving an hypothesized circumstellar dust clump). A general increase of the polarization from about 0.1% 2 days after explosion to 1.1% at $``$200 days was observed. The dominant polarization position angle (PA) defined by the loci of the data points on the Q-U diagram (Wang et al. 2001) was nearly constant at 110o from early to late observations (cf. Q-U plots of Bailey et al. 1988: Cropper et al. 1988). Since the degree of linear polarization increases with the degree of asymmetry of the ejecta, these polarimetry observations show that the scattering surface in the ejecta becomes increasingly asymmetric as we probe deeper into the core. The fact that the polarization PA is constant implies that the photosphere of the supernova has a single well-defined axis of symmetry.

More evidence of deviations from spherical symmetry is provided by the so-called “Bochum event.” During days 20-100, the H$`\alpha `$ line exhibited two additional emission-like features at a velocity of about 4500 km s-1 nearly symmetric on both the blue and red side of the H$`\alpha `$ line (Hanuschik et al. 1988; Phillips & Heathecote 1989). The emission components are most likely to be caused by 56Ni clumps ejected during the explosion (Lucy 1988; Utrobin, Chugai, & Andronova 1995). The existence of high-velocity 56Ni was also suggested by the earlier-than-expected detection of X-rays (Dotani et al. 1987; Sunyaev et al. 1987) and $`\gamma `$-rays (Matz et al. 1988).

The line of sight velocity of the Bochum clumps can be measured from the Doppler shift of the spectral profiles and compared to the photospheric velocity using arguments based on geometric obscuration to deduce the transverse velocity and space velocity (Utrobin et al. 1995). By assuming photospheric velocities of 2800 and 2300 km s-1 on days 29 and 41, respectively, Utrobin, et al. (1995) estimated the transverse velocity of the 56Ni clump to be about 2400 km s-1 and the angle between the line of sight and the vector of the absolute velocity to be 31o. These results depend on the location of the photosphere since the velocities must be such that at a given epoch the photosphere does not obscure the clump. The observed blueshift of the minima of the P-Cygni profiles of various lines show a large range. For example, on day 27.76, the velocities are -7603, -7262, -5968, -3136, -2891, -5107, -4337, -3673, -2869 km s-1 for H$`\alpha `$, H$`\beta `$, H$`\gamma `$, Fe II $`\lambda `$5018, Fe II $`\lambda `$5169, Na I $`\lambda `$5893, Ba II $`\lambda `$4554, Ba II $`\lambda `$6142, and Sc II $`\lambda `$5527, respectively (Phillips et al. 1988). The Bochum event appeared on day 20 when the photosphere was at a velocity around 3800 km s-1 as indicated by the measured velocity of the weak Fe II $`\lambda `$5018 line. Detailed NLTE spectral analysis and the Doppler shifts of weak lines gave the photospheric velocity to be around 4300 km s-1 on day 29, and 3200 km s-1 on day 41 (Lucy 1988; Höflich 1988).

These higher photospheric velocities mean that at a given epoch the photosphere was at a larger radius than assumed by Utrobin et al. and hence that the clump must be at a larger angle to the line of sight to avoid obscuration by the photosphere. We find that the motion of the 56Ni clump is broadly consistent with being tilted by roughly 45o from the line of sight. Chugai (1992) had earlier estimated the transverse velocity of the 56Ni clump to be from 4000 km s-1 to 5000 km s-1 with a representive mean of $``$ 4500 km s-1. Such a transverse velocity is in agreement with the conclusion that the 56Ni clump is directed at 45o from the line of sight. The 45o orientation we deduce for the 56Ni clump is close to the inclination of the symmetry axis of the circumstellar rings. Such a location of the 56Ni clump is also consistent with the profile of the Ni II 6.6 $`\mu `$m line observed at day 640 (Colgan et al. 1994). While we cannot determine the PA of these 56Ni clumps on the plane of the sky from the Doppler shift data alone, the strong implication is that they share a common axis with the circumstellar rings.

Recent HST images are able to spatially resolve the ejecta. Figure 1 shows an image taken in the F439W filter on 2000 June 11 (see also Jansen & Jakobsen 2001). The ejecta are obviously elongated at a PA of about 14$`\pm `$5o. At such a late stage, the images in this filter show a mixture of emission from Fe I and Fe II, with some minor contributions from hydrogen Balmer lines (Wang et al. 1996; Chugai et al. 1997). The images display the locus of the late-time excitation by 44Ti and hence show that the late-time deposition is asymmetric and substantially bipolar. The PA of the ejecta is remarkably consistent with the alignment derived from early speckle observations of the ejecta (20$`\pm `$5o; Papaliolios et al. 1989) and for the mystery spot (14$`\pm `$3o to 194$`\pm `$3o; Nisenson et al. 1987; Meikle, Matcher, & Morgan 1987; Nisenson & Papaliolios 1999). The major axis of the ejecta lies roughly in the direction of the minor axis of the inner nebular ring and of Napoleon’s Hat, but is noticeably offset by about 15o.

It is natural and expected that the 56Ni clumps of the Bochum event and the distribution of 44Ti that powers the late-time ejecta be associated. This is another argument that the 56Ni clumps share the PA of the images and the mystery spot, and in turn that these features all share the same orientation with respect to the line of sight, about 45o to the plane of the sky, away from the observer in the North.

## 4 The Structure of the Ejecta

The schematic structure shown in Figure 2 may consistently account for most of the observed features summarized above. Figure 2a shows the perspective from Earth. Figure 2b shows the structure projected onto the plane formed by the line of sight and the line from the origin to the North. The ring plane is inclined by about 44o to the line of sight. The 56Ni clump that can consistently explain the redshifted component of the Bochum event is shown in Figure 2b as well. Note that although this 56Ni clump can be consistently located on the symmetry axis of the ejecta by arguing that the 56Ni and 44Ti have similar distributions, the only direct observational constraint is that the 56Ni clump is located at an angle to the line of sight of $``$ 45o.

To test whether or not the ejecta morphology observed with the HST is consistent with early spectropolarimetry, we estimated the expected polarization PA assuming either an oblate or prolate spheroid for the scattering geometry. We examined two cases, one with the symmetry axis of the spheroid placed along the symmetry axis of the inner ring and one with the symmetry axis along the axis of the HST image of the ejecta (and of the mystery spot). For the oblate case, the PA of the polarization is co-aligned with the symmetry axis of the ellipsoid; for the prolate case, the symmetry axis implied by the polarization is at 90o from the PA of the polarization. Table 1 lists the PAs of the major axis of the inner circumstellar ring and of the ejecta and the predicted polarization PAs for the two cases of oblate and prolate ellipsoids. The polarimetry position angle defined by early broad-band observations from Cropper et al. (1987) is 129o$`\pm 6.5`$, and from Bailey (1988) is 110o, whereas the spectropolarimetry on June 3 for the 5800 Å-region defines a line at close to 0o (Cropper et al. 1988). When all the spectropolarimetric data are plotted on a single Q-U digram they define a dominant axis corresponding to a position angle of 16.5o–106.5o$`\pm 3`$ (Cropper et al. 1988). The polarization PA corrected for interstellar polarization is around 110o; this is consistent only with a prolate spheroid oriented in the same direction as the HST images of the ejecta.

This is a remarkable concordance. The early spectropolarimetry probes the geometry of the outer hydrogen envelope, while the HST images show emission primarily from metal lines, a completely different region. The fact that these different measures show similar structure is a strong indication that both observations reveal the geometry of the ejecta, despite the fact that the ejecta may be heavily obscured by dust clumps synthesized about 530 days after explosion (Lucy et al. 1989). In particular, the constancy of the polarization PA shows that the scattering geometry did not, for instance, evolve from oblate to prolate. The simplest hypothesis is that the constant polarization angle reflected the geometry of the ejecta ultimately revealed by the direct images of the ejecta. This hypothesis is confirmed by this quantitive check on the polarization PA. From the earliest epochs, the polarization of SN 1987A revealed a prolate scattering geometry that was congruent with the physical asymmetry of the iron-rich ejecta.

## 5 Kinematics of the Ejecta

At first glance, if the schematic picture of Figure 2b is true, we should expect the northern part of the ejecta which points away from us to be redshifted and the southern part of the ejecta to be blueshifted since it points toward us. Spatially resolved spectroscopy obtained using the HST shows the opposite trend. Figure 3 shows the profile of the \[O II\], \[Ca II\] 7300 Å blend. The line peak of the southern part (the bottom slice) is obviously redshifted compared to that of the northern part (the top slice) by about 1270 km s-1.

This apparent inconsistency has two important lessons to teach. In the first place, this pattern is not consistent with a spherical geometry for the ejecta. Even with the effects of dust, one would expect the same Doppler shift in the north and in the south, not a gradient, as seen. We note that the shift of line centers in Figure 3 is not monotonic from the north to south, but the central slice, representing (within positioning errors, see below) the centroid of the geometry, is blueshifted in an absolute sense and with respect to the general trend from north to south. This departure from the trend could easily be explained by a clump of obscuring dust that is thickest near the center of the ejecta.

The second lesson of Figure 3 is that the apparent discrepancy with the kinematics of the jet-like structure in Figure 2 can be resolved in the following way. The ejecta are likely to be neither spherically symmetric nor homogeneously mixed and the \[O II\], \[Ca II\] emission is probably formed in a different zone compared to the Fe I and Fe II emitting zones that are responsible for the shape of the direct images. In particular, if the density of the oxygen and calcium-rich ejecta is enhanced in the equatorial plane, then the emissivity in the equatorial plane will be higher than at the poles. This effect, plus some central dust, can easily produce the gradient in velocity from the north to south. This behavior is expected from a dusty oblate ellipsoid or torus. Detailed analyses of the spatially resolved spectra is clearly needed to map the ejecta geometry and kinematics and the dust distribution, but is beyond the scope of this paper. Nevertheless, we show in the next section that the inferred geometry – nickel clumps along the axis and an expanding equatorial torus of oxygen and calcium – are predicted by a jet-induced model for the explosion.

While an equatorial torus is consistent with the overall gradient in velocity displayed in Figure 3, we note the curious fact that the whole system seems to have a rather large positive velocity (with the exception of the possibly dust-obscured central portion). The southernmost portion has a velocity of about 1770 km s-1 and the northernmost about 500 km s-1. A small portion of this overall redshift can be accounted for by the space motion of the LMC of 290 km s-1, but that leaves about 1000 km s-1to explain. Some of this, but surely not the whole amount, could be due to a space motion of the supernova itself. We suspect that the bulk of the effect is related to the fact that the feature in Figure 3 is actually a blend of \[O II\] $`\lambda `$7320 and \[Ca II\] $`\lambda `$7300. It is not clear which is more dominant, so the rest frame wavelength (7300 Å) used in making the plot may not be correct. A 20 Å uncertainty corresponds to about 800 km s-1. There may also be some error in the STIS wavelength calibrations. We have checked this hypothesis by examining the He I/Na I d blend. If this feature is dominated by Na D, then the data agree with a torus picture with a blueshift in the north and a redshift in the south, but there is a net offset if the feature is identified as pure helium. Another contribution could come from the uncertain positioning of the slit. The overall expansion speed of the inner portion of the ejecta (note the width of the features in Figure 3 is about 3000 km s-1. While an attempt was made to position the slit so that it was centered on the ejecta, an offset by only one or two pixels to the south could mean that the light down the slit was weighted toward southern, redshifted ejecta.

The important aspect of these observations in the current context is the clear gradient in velocity from north to south. In the next section we show that this is consistent with the expectations of a jet-induced explosion.

## 6 Implications of Explosion Models of Core-Collapse Supernovae

The observed bipolar structure of the ejecta and the ejection of 56Ni at velocities above 5,000 km s-1 argue strongly that the engine that caused the supernova event was highly aspherical. A model assuming bipolar, non-relativistic jets formed at the center of the explosion that can successfully eject the supernova envelope can help us to understand how the asymmetry was produced (Khokhlov et al. 1999; Khokhlov & Höflich 2000; Höflich, Khokhlov, & Wang 2001).

As an example, in Figure 4 we give the resulting chemical distribution of the inner layers at 250 s after core collapse for an explosion of a Type II supernova triggered by a non-relativistic jet with a velocity of 11,000 km s-1 that operates over $`2.5`$ and provides an explosion energy of about $`2\times 10^{51}`$ erg in a progenitor of 15 $`M_{}`$. At about 200 s, the jet material “stalls” in the expanding material. The density contours in the subsequent phases of the explosion become almost spherical with the asymmetric chemical profiles frozen into the envelope (Khokhlov & Höflich 2000; Höflich, Khokhlov, & Wang 2001). During the phase of free expansion, the bow shock shows velocities of about 5000 km s-1. The resulting chemical composition structures are highly aspherical. Plumes of 56Ni will be selectively ejected along the jet axes and helium, oxygen, and calcium (not shown) will be concentrated in expanding tori on the equator.

The nickel plumes can possibly explain the Bochum event and the apparent prolate asymmetry of the ejecta image. In addition, the axial distribution of excitation by radioactive decay will produce an aspherical photosphere even in a nearly spherical density distribution (Wang et al. 2001; Höflich, Khokhlov, & Wang 2001). Thus if 56Ni is distributed along a bipolar axis, we expect the polarization to reflect that axis with an orientation fixed in time, as observed. The higher density of calcium in the expanding equatorial torus will give a greater emission measure near the equator and hence account for the observed northern blueshifted component.

We feel it is highly unlikely that a basically spherical model involving only Rayleigh-Taylor instabilities could reproduce the observed geometry. Jet-driven models provide an interesting alternative, although the mechanism of jet formation during the core-collapse process remains unclear. Possible mechanisms of jet formation during the collapse to form a neutron star have been discussed recently (Wheeler et al. 2000, 2002; Akiyama et al. 2002) as have asymmetries associated with neutrino flow (Shimizu, Yamada, & Sato 1994; Shimizu et al. 2001; Fryer & Heger 2000). Black-hole formation is another possibility (MacFadyen & Woosley 1999).

## 7 Conclusions

We have shown that SN 1987A has revealed a common bipolar axis of symmetry from the early photospheric phase to the late imaging of the ejecta. We have shown for the first time that the images and the polarimetry are consistent with a prolate geometry oriented nearly (but not exactly) along the axis of the inner CSM ring. We have also deduced the orientation of the nickel clumps that were observed as the “Bochum event” with respect to the line of sight and argue that these clumps are likely to be co-aligned with the bipolar ejecta revealed in late HST images. The matter that shows up in the images is precisely that heated by the long-term input of radioactive decay of 44Ti. Nagataki (2000) argues that the required abundance of 44Ti itself may be evidence for explosion by a jet-like phenomenon. Emission from \[O II\] and \[Ca II\] does not contribute substantially to the HST images, but is shown in resolved spectroscopy to be blueshifted in the North where the bipolar flow and one component of the Bochum event are directed away from the observer. The existence and orientation of the Bochum event, the common axis of polarization PA with the late-time image, and the northern blueshifted oxygen and calcium are all consistent with, and even predicted by, jet-induced supernova models. In these models, radioactive matter is ejected primarily along the poles and elements synthesized in the progenitor (He, O, Ca) are ejected in expanding equatorial tori by shocks that propagate away from the bow shock of the jets and converge on the equator. The synthesis of observations we have presented here is consistent with a model in which SN 1987A underwent a bipolar jet-induced explosion along a common axis that ran from the photosphere deep into the ejecta.

The small difference of 15o between the symmetry axis of the ejecta and the inner CSM ring may provide some additional clues to the overall evolution and explosion of SN 1987A. If the axes of symmetry are all related to the rotation of the progenitor, then we should expect a single axis during all these processes in order to conserve angular momentum. What has changed the symmetry axis, and therefore carried away angular momentum from the system? Could it be that the explosion was dominated by a magnetic field and the inner ring was dominated by rotation? Or vice versa? Finally, we note again that the mystery spot closely aligns with the ejecta and the prolate axis of the polarization geometry. The mystery spot still begs explanation (Rees 1987; Piran & Nakamura 1987). Perhaps the geometry we present here can provide new clues to that problem as well. These are all open questions to be pursued with the strong suspicion that SN 1987A was a bipolar explosion.

Acknowledgements

Support for this work was provided by NASA through grants GO 08243 and 09114 from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26255. Additional support was provided by NASA through grant NAG5-10766 and by NSF through grant AST-0098644 to the University of Texas.


text:
# Caloric Curves and Nuclear Expansion

## I Acknowledgements

This work was supported by the U S Department of energy under Grant DE-FG03- 93ER40773 and by the Robert A. Welch Foundation.