text:
# Phase diagram of heavy fermion systems

## 1 Monochalcogenides of Sm and Tm: conduction, valence, and magnetism

In the monochalcogenides of anomalous rare earths (RE) like Sm or Tm the equilibrium between two valence states (Sm2+/Sm3+ or Tm2+/Tm3+) governs the release of a 5$`d`$ electron . Let us start with our new guided tour of the canonical example of SmS. The full line $`T_{BG}`$ of the ($`TP`$) phase diagram (see Fig.1) represents the first order isostructural transition between the black (B) insulating (I) Sm2+ phase and the so called gold (G) phase where the Sm ions jumps to an intermediate valence state $`v2.7`$. In this gold phase, smoothly increasing the pressure must lead to recover a trivalent Sm3+; it is a Kramers ion, long range magnetic order must appear above this threshold $`P=P_c`$.

Furthermore a sound experimental fact is that up to $`P_\mathrm{\Delta }=20`$ kbar in the low pressure gold phase a many body insulating phase is smoothly formed on cooling. Whatever the sample quality and the pressure conditions, $`P_\mathrm{\Delta }`$ is a fixed pressure above which the sample is clearly metallic i.e. its resistivity $`\rho `$ decreases on cooling after reaching a maximum $`\rho _{max}`$ at $`T_m`$ . The recent combination of macroscopic measurements by transport and ac specific heat and of the microscopic probe of the hyperfine field on the Sm nuclei by nuclear forward scattering (NFS) using ESRF facility demonstrates that long range antiferromagnetic (AF) order appears at $`P_cP_\mathrm{\Delta }`$. As the pressure variation of $`T_m`$ now identified as the Néel temperature $`T_N`$ is huge, $`P_c`$ is presumably a first order point associated with a lattice parameter discontinuity i.e. a tiny valence jump at $`P_v=P_c`$ towards the trivalent state.

By contrast to the Sm case where the valence fluctuations occur between the Sm2+ state with a zero angular momentum ($`J=0`$) and Sm3+ with a $`J=5/2`$ angular momentum, for the Tm ions the mixing is between two finite angular momentum states, respectively $`J=7/2`$ (Tm2+) and $`J=6`$ (Tm3+). Whatever the pressure i.e. the valence and even the electronic conduction, the ground state is magnetic. However, the similarity with SmS is a link between magnetic properties and 4$`f`$ localisation. For TmTe as shown on Fig. 2, the $`P`$ sequence of ground states is AF/I up to $`P_{BG}`$, ferromagnetic (F) and metallic (M) up to $`P_v`$ and finally the entrance to a AF metallic state ; the rapid collapse of $`T_{Curie}`$ just below $`P_v`$ and the fast increase of $`T_N`$ just above $`P_v`$ suggests a first order transition. For TmSe, already in the gold phase at $`P=0`$, the switch at $`P_v`$ occurs between an AF insulating phase $`P<P_vϵ`$ and another AF metallic state above $`P_v`$ . TmS itself appears already at $`P=0`$ as a quiet quasitrivalent metallic state, in principle the ideal Kondo lattice .

## 2 Cerium heavy fermion systems: valence and magnetism

As pointed out already, the cerium heavy fermion problem can be summarized by two major ($`T,P`$) phase diagrams. The wellknown $`T_{\alpha ,\gamma }`$ first order transition of the cerium metal and the up to date magnetic phase diagram $`T_N`$ or $`T_{Curie}`$ versus pressure where the ordering temperature may collapse at $`P_c`$ . For the cerium metal, the treatment of the magnetism can be crude as all characteristic temperatures $`T_{\alpha ,\gamma }`$ are high by comparison to any hypothetical $`T_N`$ or any fine structure as the crystal field splitting $`C_{CF}`$ of the trivalent $`\gamma `$ phase. At $`P=0`$, $`T_{\alpha ,\gamma }(0)`$ is already positive, finite and near 100 K and the critical end points $`T_{cr}=700`$ K, $`P_{cr}=22`$ kbar are large. For typical cerium heavy fermion compounds $`P_c`$ is lower or comparable to $`P_{cr}`$ and of course $`T_{\alpha ,\gamma }(P=0)`$ and $`T_{cr}(P_{cr})`$ will be negative. However, the 4$`f`$ electron will feel the distance from $`P_{cr}`$ and cause the crossing through a characteristic pressure $`P_v`$ where the local susceptibility loses the anisotropy given by its reaction to the crystal field (that corresponds to the well known condition $`k_BT_K`$ (Kondo energy) greater than $`C_{CF}`$). Above $`P_v`$, the 4$`f`$ electron stays rigidly linked to the full degenerancy of the $`J=5/2`$ angular momentum. At $`P_v`$, the sensitivity of the 4$`f`$ electron to the crystal field has vanished; $`P_v`$ is characteristic of a valence crossover. The relative location of $`P_c`$ and the hidden parameter $`P_v`$ is one of the key issue. The idea that singular properties is underpinned by a quantum critical point i.e. pressure is one of the trends of strongly correlated electronic systems (see references  for the hidden order of URu2Si2 or  for the metamagnetic quantum criticality in Sr3Ru2O7 and  for the conversion of a Kondo gas to heavy electron Kondo liquid).

Recently special focus has been given on systems just on the paramagnetic verge of $`P_c`$ i.e. at $`P_c+ϵ`$ . The Non Fermi liquid (NFL) label refers to the difficulty to reach the standard Fermi liquid regime where usual laws are found, even for non-interacting Fermi gas (specific heat $`C=\gamma T`$ linear in temperature $`T`$, susceptibility in $`T^2`$, scattering rate proportional to $`T^2`$); the NFL label guarantee “vintage” describes a very large temperature crossover before reaching the Fermi liquid graal or a new state of matter.

To clarify the complexity of this strongly correlated electronic system, an old fashioned approach via thermodynamic languages will be carried out . In the same spirit as the Clapeyron and Ehrenfest relations on the $`P`$, $`T`$ feedback for first or second order transitions, the Grüneisen parameter

$$\mathrm{\Omega }^{}(T)=\frac{\alpha }{C}\frac{V_0}{\kappa }$$
(1)
(where $`\alpha `$, $`V_0`$ and $`\kappa `$ are respectively the volume thermal expansion, the molar volume and the isothermal compressibility) is a unique experimental gauge of interfering hierarchies. It will be reduced to a constant $`\mathrm{\Omega }^{}(0)`$ only if the free energy $`F`$ can be expressed by a unique parameter $`T^{}`$ i.e.

$$F(T)=T\mathrm{\Phi }(\frac{T}{T^{}})\mathrm{with}\mathrm{\Omega }^{}(0)=\frac{\mathrm{log}\mathrm{T}^{}}{\mathrm{log}\mathrm{V}}$$
(2)
Figure 3 represents the temperature variation of the Grüneisen parameter of CeRu2Si2 at $`P=0`$ i.e. roughly 3 kbar above the critical pressure $`P_c=3`$ kbar. The two singular points are: (i) the huge extrapolated value of $`\mathrm{\Omega }^{}(0)=+190`$ and (ii) the slow entrance in a simple regime ($`T1`$ K) where $`\alpha `$ and $`C`$ reach their proportionality . The pressure dependence of $`\mathrm{\Omega }_0^{}`$ is shown in the insert of fig.3. At low pressure close to $`P_c`$, $`\mathrm{\Omega }^{}(0)`$ is weakly pressure dependent changing from 190 to 220 for a negative shift of $`P=2`$ kbar  while $`\mathrm{\Omega }_0^{}`$ drops to 80 at $`P_v`$ . Other values of $`\mathrm{\Omega }_P^{}(0)`$ can be listed crudely for other heavy fermion system. In a double logarithmic representation of the volume dependence of the $`A`$ coefficient of the well known $`AT^2`$ Fermi liquid resistivity assuming the usual Kadowaki-Woods rule $`A\gamma ^2`$, $`\mathrm{\Omega }^{}(0)`$ is directly linked to the slope of the curve .

## 3 Magnetic field switch from AF correlation to F interactions

The compound CeRu2Si2 is a quasi-ideal Kondo lattice with a nice ordered tetragonal lattice, a robust Ising character of the local susceptibility, the realization of a clean limit condition i.e. a mean free path higher than 1000 Å far larger than an hypothetical superconducting coherence length with temperature $`T_s=200`$ mK and good enough to obtain the complete determination of the Fermi surface (FS) at least above $`P_c`$ (see references in ).

The application of a magnetic field $`H`$ along the easy $`c`$ axis leads to a switch through a drastic continuous crossover at $`H=H_m`$ from a low field paramagnetic phase (Pa) dominated by AF correlations at a finite wavevector $`k_0`$ to a highly polarized state (PP) dominated by the low wavevector $`q`$ excitation. The microscopic vision of CeRu2Si2 was nicely given by inelastic neutron scattering experiments, (i) with at $`H=0`$, the balance between the local response and the intersite coupling, (ii) with, under magnetic field, the continuous spreading of the AF response with the same characteristic glumly energy $`\omega _{AF}1.6`$ meV up to $`H_m`$  and the observation just in the vicinity of $`H_m`$ of an inelastic ferromagnetic signal at far lower energy transfer $`\omega _F0.4`$ meV than $`\omega _{AF}`$ (see Fig.4) .

A skilful mechanism characteristic of the Kondo lattice CeRu2Si2 controles the “ping-pong” between the magnetic interactions. Quite remarkably, the field and temperature response is well reproduced by a simple form of the thermodynamic as the entropy.

$$S=S(\frac{T}{T^{}},\frac{H}{H_m})$$
(3)
with equal Grüneisen parameters $`\mathrm{\Omega }_T^{}=\mathrm{\Omega }_{H_m}=190`$. With this hypothesis one can predict the field variation of $`\gamma _H`$ knowing the field variation of the linear thermal expansion $`\alpha _v=a_HT`$. A satisfactory agreement is found between this scaling derivation of $`\gamma _H`$ and the bare determination either by specific heat or by magnetization. Physically, the adaptive formation of the many body Kondo lattice bands is characterized by a pseudogap . The Ising spin character of the bare local magnetic ion plays a capital role in the sharpness of this electronic substructure and correlatively for the pseudo-metamagnetism . The pseudogap ingredient is a key input in any theory which needs to respect the scaling behavior . Reminiscent of the $`\alpha \gamma `$ collapse of the cerium metal, the associated spectacular lattice softening ($`30`$%) is again in agreement with the push pull between magnetic, electronic and lattice instabilities  and the underpinned vicinity of a quantum critical end point.

## 4 To be or not to be? Speculations on superconductivity

After a lecture describing the “beauty” CeRu2Si2 to quantum criticality, the first attack is to produce the proofs that CeRu2Si2 is close to a magnetic instability. Looking back to the primitive results one decade ago on the CeRu2Si2 series either sweeping electronic density (doping or pressure) or magnetic field, a NFL increase of $`C/T`$ on cooling in a given temperature window justifies the proximity to $`P_c`$ . Furthermore, this variation is not drastically different from that popularized in the canonical example of CeCu5.9Au0.1  after a rescaling of the temperature by some phenomenological parameter describing the proximity to $`P_c`$ and also differences in physical ingredients. As on approaching $`P_c`$, no drastic pressure dependence of $`\mathrm{\Omega }^{}(0)`$ occurs. Collapsing to $`P_c`$ on the paramagnetic side does not lead to an “earthquake” even if of course $`T^{}`$ will collapse when $`T0`$. The deviation from any singular scaling law of the imaginary part of the dynamical susceptibility $`\chi ^{\prime \prime }T^\alpha \mathrm{\Phi }(\frac{\omega }{T})`$ becomes also quite difficult to point out. The yet unanswered dilemma is if a new state of matter occurs right at $`P_c`$.

The obsessive focus just at $`P_c`$ or on the paramagnetic (Pa) side ($`P>P_c`$) leads to neglect the disappearance of long range magnetism. As many tunings were realized by doping i.e. at the cost of an increase of disorder, $`P_c`$ appears as a second order quantum phase transition. Comparison with the previous chalcogenide cases, recent NMR experiments on CeIn3 close to $`P_cϵ`$ with evidence of AF and Pa magnetic phase separation  as well as drastic Fermi surface changes at $`P_c`$ in CeRh2Si2  point out that at the limit of a pure lattice at $`P_c`$ the transition may be first order. Today one key question is to determine precisely the Fermi surface on both sides of $`P_c`$ knowing the two limit conditions of a localized 4$`f`$ picture at low pressure ($`PP_c`$) and an itinerant 4$`f`$ picture at high pressure ($`PP_c`$). Progress will imply a deep interdisciplinary collaboration between experimentalists and theoreticians eager to track carefully the $`P`$ and $`H`$ variations including the modification of the Brillouin zone and the differentiation between spin up and spin down effective mass carriers.

For unconventional superconductivity, the nature of the phase transition at $`P_c`$ is not critical as proved by the occurrence of superconductivity in CeIn3 and CeRh2Si2. There are examples like CePd2Si2, CeIn3 and CeRh2Si2 where the superconducting pocket is stocked to $`P_cP_v`$; there are cases like CeCu2Si2 , CeCu2Ge2  and CeNi2Ge2  where two superconducting domains are centered on distinct pressures $`P_c`$ and $`P_v`$. We will test soon if superconductivity appears on the verge of $`P_v20`$ kbar in CeRu2Si2; near $`P_c=3`$ kbar, the Ising character of the spin allows only longitudinal magnetic fluctuations which does not favor the observation of $`d`$ wave pairing  while near $`P_v`$ this superconductivity can be promoted by the emergence of transverse magnetic modes as well as valence fluctuations .

It is time to concentrate on the complete bundle of curves on each side of $`P_c`$ in CeRu2Si2 as well as in the other cerium tetragonal lattices and Yb homologues following our Dresden colleagues .

J. F. wants to thank the Lorentz center for its invitation to the workshop “Non Fermi liquid behavior and quantum phase transitions” (May 2003) and FERLIN ESF support during four years.


format: text

text:
# Reply to Comment on: “Rashba precession in quantum wire with interaction”

## Abstract

It is shown that the recent Comment  by Y. Yu on the above article is not substantiated.

Recently a Comment appeared on the ArXiv preprint server  on my article ‘Rashba precession in quantum wires with interaction’ . Its main statement is to question the dependence of the Rashba spin-orbit generated persistent spin current on the interaction. I demonstrate now that Comment  contains numerous faults, one of which leading to the above (mis)conclusion.

The Author begins by stating that non-integer valued $`J_\nu `$ would contradict to the seminal work by Haldane . In this regard  carefully distinguishes between eigenvalues and expectation values at non-zero $`\alpha `$ and it are the latter for which Eq. (2) of  is correct (as clearly stated before Eq. (7) in ). It should be pointed out that the expression $`\stackrel{~}{J}_\sigma `$ found after Eq. (9) in  is non-integer either ($`q_R=m\alpha `$ is tuned externally), contrary to what is stated one sentence later.

Secondly, the Fermi velocities $`v_{Fa}`$ in  are taken as spin dependent which in the effective mass description is erroneous at given Fermi energy. The effective mass description holds up to the perturbative order $`𝒪(\alpha ^5)`$, as is carefully derived in , cf. also . The $`v_{Fa}`$ in  describe a non-equilibrium situation.

The main error, however, occurs in Eq. (10) of  where the shifted current $`\stackrel{~}{J}_\sigma `$ is inserted into the interaction, instead of the current $`J_\sigma `$. For not explicitly spin dependent interactions the $`J_\sigma =0`$ state is of course the state of lowest interaction energy and not the state $`\stackrel{~}{J}_\sigma =0`$. Correct insertion would immediately lead to the (correct) result, Eq. (2) of . In Eq. (10) of  the interaction depends on the Rashba coupling $`q_R`$ while the Rashba spin-orbit term depends explicitly on the interaction strength. Both is clearly not representing the model under consideration, the Rashba spin-orbit coupling being a single particle operator, cf. Eq. (1) of .

In Eq. (14) the prefactor $`v_F/\lambda _\sigma `$ has simply been written ad hoc in front of the $`J_\sigma `$-linear term, without further justification. Correct would have been this term proportional to $`q_R`$ precisely as in Eq. (14) with $`\lambda _\sigma =1`$, describing indeed SU(2) symmetric spin sectors (a case which has not been considered in ). It is not the scope of this Reply to correct for the inconsistencies between Eq. (14) of  and the work by Y.-S. Wu , where the Author of  is Coauthor.

In conclusion, the persistent spin current and the Rashba length both being ground state properties (and as such unrelated to $`H_q`$) do depend on the interaction. The Comment  by Y. Yu lacks substantiation.


text:
# Long Term and Short Term Effects of Perturbations in a Immune Network Model

## 1 Introduction

The main task of the immune system is to protect the organism against dangerous elements: antigens (virus, bacteria, poison, cell residues, etc). Depending on the antigen the immune system may elicit different kinds of responses: the cell-mediated immune response or the humoral response. The models discussed in this paper are related to the humoral responses generated by B cells (one of the main classes of lymphocytes), which are the cells that produce the antibodies. The antibodies produced by a given population of B cells are replicas of its molecular receptor. Each molecular receptor exhibited by a given B cell population recognize different recognition sites (epitopes) of the antigen by lock-key interactions.

The clonal selection theory  is the most accepted theory about the B cells and was proposed by Burnet in $`1959`$. It states that the antigen chooses by pattern recognition the clones of B cells (population of B cell and antibodies) that will proliferate. In order to recognize any foreign (or dangerous) element the immune repertoire must be complete. According to estimates the human immune system is able to express at least the order of $`10^6`$ different receptors . Due to the completeness of the repertoire, the immune system would be able to recognize and be recognized (recognizing epitopes of its own antibodies), therefore the same mechanism of recognition should work for both antibody-antigen and antibody-antibody reactions. In 1974 Jerne , taking into account these different mechanisms, suggested that when the antigen is presented to the organism it will activate a set of clones of B cells, leading to the production of specific antibodies; those antibodies will in turn be recognized by a second set of clones activating them, and so on. Due to the interplay of the mechanisms of activation and suppression, this chain of reactions will be finite, preventing the “information” from “percolating” through the entire system. Therefore the immune system functions as a functional network , with complexity comparable to the nervous system. Since its proposal, only few evidences support the existence of the immune network theory , and those evidences suggest that if the network exists only $`20\%`$ of the lymphocytes will be activated, while the rest of the clones will form a pool of immunocompetent lymphocytes that are able to recognize any antigen.

In what follows we define in section 2 a model introduced by Stauffer and Weisbuch in 1992 , which incorporates the main concepts of Jerne’s immune network theory. In section 3 we review the main results obtained in the literature, focusing on both the physical  and biological aspects , since among many other attempts to model what an immune network could be, this was the first that could reproduce the behavior of a real immune system. In section 4 we present new results about the short and long time behaviors of the multi-connected network under small and large perturbations. We also discuss the results obtained about the aging effects in a very recent work  on the light of these new results. In our concluding remarks (section 5) we address the future prospects from the physical and biological points of view.

## 2 The Model

In 1992 Stauffer and Weisbuch (SW)  introduced a cellular automata model to describe the evolution of the B cell repertoire. This model was inspired in a previous one (based on a difference equation approach) introduced by de Boer, Segel and Perelson (BSP)  using a shape-space formalism to describe a large-scale immune network. Using a discrete shape-space formalism SW associate each point of a $`d`$-dimensional space to a different molecular receptor (or a clone of B cell). In this way each receptor is characterized by $`d`$ properties, and its neighbors in the shape space will correspond to molecular receptors that differ from it by one of these properties (according to the estimates based on a continuous approach , if the notion of shape space is relevant, then $`d5`$).

To each receptor they associate a three-state cellular automaton $`B(\stackrel{}{r},t)`$ that will describe the concentration of the population characterized by this receptor $`\stackrel{}{r}`$ over the time: low ($`B(\stackrel{}{r},t)=0`$), intermediate ($`B(\stackrel{}{r},t)=1`$) and high ($`B(\stackrel{}{r},t)=2`$). Following the ideas of the immune network theory, the interactions among different populations are based on complementarity (lock-key interaction) and are defined by deterministic rules. The time evolution of the cellular automaton is based in a deterministic non-local rule: population $`B(\stackrel{}{r},t)`$ at site $`\stackrel{}{r}`$ is influenced by the populations at site $`\stackrel{}{r}`$ (its mirror image or complementary shape) and its nearest-neighbors ($`\stackrel{}{r}+\delta \stackrel{}{r}`$) (representing defective lock-key interactions). The influence on the population at site $`\stackrel{}{r}`$ caused by its complementary populations is described by the field $`h(\stackrel{}{r},t)`$:

$$h(\stackrel{}{r},t)=\underset{\stackrel{}{r}^{}(\stackrel{}{r}+\delta \stackrel{}{r})}{}B(\stackrel{}{r}^{},t)$$
(1)
where for a given $`\stackrel{}{r}`$ the sum runs over the complementary shape $`\stackrel{}{r}^{}=\stackrel{}{r}`$ and its nearest neighbors. Due to the finite number of states of the population $`B`$, the maximum value of the field $`h(\stackrel{}{r})`$ is $`h_{max}=2(2d+1)`$. The rules are based on a window of activation for each population which is inspired in a log bell shaped proliferation function associated to the receptor cross-linking involved in the B cell activation . There is a minimum field necessary to activate the proliferation of the receptor populations ($`\theta _1`$), but for high doses of activation (greater than $`\theta _2`$) the proliferation is suppressed. The updating rule may be summarized as:

$$B(\stackrel{}{r},t+1)=\{\begin{array}{cc}B(\stackrel{}{r},t)+1\hfill & \mathrm{if}\theta _1h(\stackrel{}{r},t)\theta _2\hfill \\ B(\stackrel{}{r},t)1\hfill & \mathrm{otherwise}\hfill \end{array}$$
(2)
but no change is made if it would lead to $`B=1`$ or $`B=3`$. We define the densities of sites in state $`i`$ at time $`t`$ as $`B_i(t)`$ ($`i=1`$, 2, 3).

The initial configurations are randomly generated according to the parameter control $`x`$, which determines the initial concentrations: $`B_1(0)=B_2(0)=x/2`$, while the remaining $`L^d(1x)`$ sites are initiated with $`B(\stackrel{}{r},0)=0`$.

## 3 Results obtained from previous studies

The above definition of the model corresponds already to the modified version which has been proposed in Ref. . SW have shown that for the original version of the model , there is a stable-chaotic transition on the behavior of the discrete system for $`d4`$, when varying the parameter $`x`$. As pointed out by the authors, according to the characteristics of those behaviors, none of them would be appropriate to describe the evolution of the immune repertoire.

One of us has studied how the system attains the chaotic regime, by studying the transient times and periods of the original model close to the transition . It was observed that close to the transition, the system is trapped in cycles with large but finite periods. In 1995 Zorzenon dos Santos and Bernardes  introduced the modified version of SW model described above and have extensively studied the behavior of this model in the parameter space. Differently from SW, they obtained  a stable-chaotic transition for $`d2`$, and mapped all the behaviors on “phase” diagrams for different set of parameters, showing that there was a broad transition region between stable and chaotic behaviors.

In the following work  Bernardes and Zorzenon dos Santos have investigated the dynamical behavior of the model on this broad transition region. They have observed an aggregation-disaggregation dynamics, with clusters splitting and fusing along the time, as a multi-connected network. The authors have also studied the behavior of the system in the transition region when subjected to antigen presentation. The antigen presentation was simulated by flipping randomly chosen populations from the inactive state to the highly activated state, reproducing in this way the activation of the populations caused by the presence of an antigen. The results obtained by the standard spread of damage procedure  indicate that after the antigen presentation, the perturbation first increases and then relaxes after a few time steps, indicating that some information about this perturbation is incorporated by the system. These results suggested that this model would be a good candidate to model a real immune system behavior. Zorzenon dos Santos and Bernardes then applied it  to reproduce immunization and aging experiments performed with mice under multiple antigen presentations  obtained by the group of Nelson Vaz.

In the experiment of immunization  six mice of the same linage are subjected to the following protocol: the antigen presentations are produced by intra-venal injections of ova. The time interval between the first and second presentation is 14 days, but after that the time interval between consecutive presentations is fixed to 7 days. Before each antigen presentation (except for the first one) they measure the amount of specific antibodies in the blood. Among many experiments reported by this group, Zorzenon dos Santos and Bernardes chose the one showing a refractory (saturation) behavior related to the immunization protocol  that could not be explained by clonal selection theory. However, the aging experiment chosen  is one among others obtained by the group showing the same result: a reduction on the intensity of the response as the system ages.

In order to simulate the immunization protocol , for each sample they let the system evolve from its initial configuration to $`1000`$ time steps, and associate this time step to the birth of the simulated mouse. Then random small perturbations, in size and location in the shape-space (using different time intervals between them), are produced in order to simulate the noise to which the immune system is subjected due to the environment, food ingestion, etc. By using the arbitrary scale of 1 day corresponding to 5 time steps, it is possible to simulate independent antigen presentations without over-exciting the system . The immunization protocol is simulated by large perturbations (one order of magnitude greater than the small ones) being produced periodically every $`7`$ days. The protocol was simulated for different samples ($`10`$) of young and old mice (8 and 24 weeks, respectively). The results could be interpreted on the light of the immune network theory, the refractory behavior being associated to the saturation of the multi-connected network in incorporating information about the perturbation, and the aging effects related to the loss of plasticity of the system.

Since there are only few evidences that support the existence of the immune network, this model played a very important role in making the connection between some experimental results that could not be explained by clonal selection theory, and the immune network theory. Moreover the theoretical results indicate that under multiple antigen presentation the saturation is due not only to the number of specific antibodies produced (as shown in the Elisa experiments) but also the ability to incorporate information by including new populations in the network. Since the number of populations belonging to the network does not change significantly during the time evolution, in order to add information, part of the information already present in the network should be lost. Therefore the profile of expressed antibodies would change after each antigen presentation, reflecting such changes. This is an aspect that comes out of the theoretical results and should be investigated by the experimentalists when the appropriated tools become available. The techniques currently available do not allow to identify the difference among the specific antibodies. That would be a requirement to verify whether they change after multiple presentations as suggested by the model. However, there are some experimental evidences that it may happen .

Once the complex multi-connected network generated by the model can be used to reproduce the behavior of a real immune system, it becomes interesting to investigate its dynamical properties from the physical point of view. This investigation will allow to learn about the dynamics of a (real) complex system operating out of equilibrium. The characterization of such behavior started with the work of Bernardes and Zorzenon dos Santos , who analyzed the distribution of cluster sizes and the distribution of the permanence time (i.e. the time interval during which each population remains activated or belong to the multi-connected network). They found a characteristic cluster size and a power law behavior for the distribution of permanence times. The characteristic cluster size is associated to the loss of plasticity while the power law distribution indicate that the populations belonging to the network have no typical permanence time, reflecting the fact the dynamical memory is generated by incorporating information about the different antigens presented to the system.

In a very recent work in collaboration with D. Stariolo  we have used the auto-correlation functions (a common tool to study aging effects in e.g. glassy systems) in order to analyze the similarities between the aging effects in the multi-connected network and in glassy systems. While in glassy systems the aging effects results from the loss of plasticity generated by frustration mechanisms, in the network these effects are caused by multiple perturbations, since in order to incorporate information the system has to adapt to (or satisfy) the mechanisms of activation and suppression of the dynamics. The usual procedure adopted to study glassy systems consists in taking a “picture” of the system at time $`t_w`$ and calculating the number of sites that do not change during the following time steps. In this sense the auto-correlation function is equal to $`C_{tot}(t,t_w)=1hd(t,t_w)=1HD(t,t_w)/N`$, where $`hd(t,t_w)`$ \[$`HD(t,t_w)`$\] is the normalized (non-normalized) Hamming distance between the configurations at times $`t`$ and $`t_w`$. According to the results obtained, the system without any perturbation is driven to a long period attractor after a long transient time ($`10^410^5`$ time steps). All the biologically relevant effects are observed in the transient time of the system. In the purely deterministic case, the transient is simply the time it takes for the system to reach the attractor ($`10^4`$ time steps). When subjected to random small perturbations, however, the very notion of a transient becomes fuzzy: results in Ref.  show that the system approaches a cycle, but is deflected from it by the small perturbations after $`10^310^4`$ time steps. Therefore, the small perturbations will cause the system to change attractors from time to time due to their cumulative effects. This is reflected in the decrease of the auto-correlation functions, as can be seen Fig. 1, where the system leaves the cycle it had approached until $`t_w=10^5`$ time steps (these results are discussed in detail in Ref. ). The analogue observed in glassy systems corresponds to changes of the system to different minima of the energy landscape during the relaxation time. Contrary to one’s initial intuition that the large perturbations would accelerate the de-correlation process, the changes induced by large perturbations do not lead to this effect. Curiously, the large perturbations alone (which correspond to the immunizations) lead to a much weaker (slower) de-correlation: this is due to the fact that they are produced always at the same sites in shape space. Small perturbations can be more easily absorbed by the system than the large ones since they involve only local changes. The results from the autocorrelation functions therefore indicate that the system first approaches a cycle (phase space compression), but after some time the small perturbations eventually deflect it from the trajectory to its “natural” attractor .

## 4 Results

This scenario provides a possible explanation for the results regarding the spread of damage in the system, which managed to reproduce the experimental data from immunization experiments with mice . Fig. 2 shows the time evolution of the Hamming distance between a system which undergoes only small perturbations, and its copy, which undergoes the same small perturbations and the large periodic immunizations. All curves saturate after some time, just like the Elisa measurements . More important, however, is the fact that the value of the saturation depends on the age at which immunization started. As the age increases, the saturation value decreases, a result which was also observed in the mice and is interpreted as a sign of loss of plasticity with age . This might be related to the fact that the system is in its transient, while trying to reach its long-period attractor. The “older” the system is, the less configurations available there are.

If we now repeat the computer experiment of immunization (with small and large perturbations) for ages greater than $`10^3`$ time steps, we obtain a different result: for older systems the stationary Hamming distance may increase, as shown in the rightmost curve of Fig. 2. We argue that the cumulative effects of the small perturbations could be responsible for the change of behavior in the rightmost curve of Fig. 2: for sufficiently old systems, the “noise” of the small perturbations become more important than earlier in the transient, and the monotonicity of the stationary value of the HD with the age at which immunization starts, no longer holds. This result is summarized in Fig. 3. It is interesting to point out that this breakdown of the refractoriness induced by age, takes place on time scales which are comparable (typically larger than) the life time of the mice.

The large perturbations, on the other hand, propagate to the entire system and some of the active populations generated during the antigen presentation will remain as part of the network. This mechanism depends on the plasticity of the system or, in other words, will be controlled by the characteristic cluster size of the system . When we combine small and large perturbations, the system is driven to a cycle that incorporates information about the antigen, but will eventually be deviated from that trajectory as soon as the cumulative effects of the small perturbations dominate the dynamics . In order to observe the instantaneous effects of both types of perturbations, we have analyzed the normalized Hamming distance between the same system at two consecutive time steps, focusing only on the subspace of sites with $`B(\stackrel{}{r},t)=2`$. This quantity is denoted by $`HD_{22}(t,t+1)`$ and referred to simply as “instantaneous Hamming distance”. It is plotted in Fig. 4 as a function of time for different situations.

We start the analysis of Fig. 4 by noticing the behavior for the purely deterministic case (lower curve): the plot indicates that typically $`2.6\%`$ of the active network (sites with $`B(\stackrel{}{r},t)=2`$) changes at each time step. The end of the transient is clearly seen at $`t15000`$ time steps (for this particular realization): the curve then becomes periodic, with a very long period. Note that the periodicity of the dynamics can be confirmed by checking that the $`C_{tot}(t_w,t)=1`$ for some $`t>t_w`$ and for sufficiently large $`t_w`$ (see Fig. 1 and Ref. ). To illustrate the periodicity of $`HD_{22}(t,t+1)`$, we plot in Fig. 5 the return map of the time series from $`t=20000`$ to $`t=25000`$ time steps (upper plot). The extremely long period might lead one to suspect, based on the return map, that the system might be chaotic. This suspicion is ruled out when one observes the return map from $`t=25000`$ to $`t=30000`$ time steps (lower plot). We challenge the reader to find a discrepancy between the two return maps.

Returning to Fig. 4, notice that the size of the changes induced by the small perturbations changes near the end of the transient time. While the perturbed system no longer attains the cycle, it can still feel the effects of the attractor: when close to the cycle, the small perturbations induce larger changes in the active network at each time step. If subjected only to the large immunizations, the instantaneous Hamming distance take typically larger values. This might suggest a contradiction with the results in Ref. , according to which the system only with large immunizations de-correlates slower. But the contradiction is only apparent: $`HD(t_w,t)`$ measures a long term effect, which is dominated by the small perturbations. Fig. 4, on the other hand, is a measure of the instantaneous change occurring in the system, i.e. a short time effect. Finally, one can combine both types of perturbation, the result being even larger instantaneous changes.

Another aspect of Fig. 4 that should be noted is an interesting interplay between the perturbations. In what follows, we refer to the typical values of the four curves in Fig. 4 (for large $`t`$), from bottom to top, as $`d`$ (“deterministic”), $`s`$ (“small”), $`l`$ (“large”) and $`ls`$ (“large and small”) — therefore $`d<s<l<ls`$, as discussed above. The fact that $`lss>ld`$ indicates that the effect of the large immunizations (i.e. the instantaneous change it induces) is boosted by the presence of the small perturbations. The reverse is obviously also true: the previous inequality implies $`lsl>sd`$, indicating that the small perturbations induce larger changes in the network if the large immunizations are present. This signals a cooperative interaction between two kinds of perturbations which are significantly different in nature, leading to the necessary complexity of behavior which allow the model to reproduce experiments performed with real immune systems.

## 5 Concluding Remarks

In this paper we have reviewed the main results from the biological and physical point of view of a cellular automata model introduced by SW in 1992. Since it incorporates the main features of the immune network theory and models the functioning of a multi-connected network, it emerged as a candidate to reproduce the basic dynamics of what the real immune network could be. The model was then used to model real experiments performed with mice about immunization and aging. The immunization experiments are simulated by producing small and large perturbations or damages. In this paper we complement the results of a very recent work  with the study of the role of each kind of perturbation. The results obtained in the previous work, using auto-correlation functions, apparently suggested that the small perturbations would be more important than the large ones from the dynamical point of view, since their cumulative effects could change the attractor of the dynamics after a few thousand time steps. However this study was focusing in the long term behavior of the system. Here we have studied its short term behavior and shown that both types of perturbation are important: it is their cooperative effects that generate the necessary complexity which allowed to reproduce experiments performed with real immune systems.

Despite the understanding of the dynamics we have achieved up to this point, there are still some points that should be investigated: what happens when changing the frequency of presentation of the small perturbations for a fixed value of the large damages, and how these results will change when varying the large damage size? This study (already in progress) will draw a scenario that would allow us to try to apply this model to reproduce other experimental results concerning, for instance, cross-reactivity mechanisms observed in immunization protocols, and tolerance.


text:
# GRB 980425, SN1998BW AND THE EMBH MODEL

## INTRODUCTION

The aim of this talk is to present the application of the EMBH theory, previously successfully applied to GRB 991216 used as a prototype, to the case of GRB 980425 (Pian et al., 2000) and SN1998bw (Galama et al., 1998). This is a particularly important test for the validity of the EMBH theory over a range of energies of 6 orders of magnitude: as we will see, both sources appear to be spherically symmetric and the respective total energies are $`E_{tot}5\times 10^{53}\mathrm{ergs}`$ and $`E_{tot}10^{48}\mathrm{ergs}`$. We recall that the EMBH theory (see Ruffini et al., 2003b) depends only on three parameters, the energy of the dyadosphere $`E_{dya}`$, the $`B`$ parameter and the factor $``$ describing the interstellar medium (ISM) porosity. The theory, therefore, explains all the observed features of the bolometric intensity variations of the afterglow as well as the spectral properties of the source and, in the specific case of GRB 980425, it also allows to clarify the general astrophysical scenario in which the GRB actually occurs. In this system, in fact, we propose that GRB 980425 has been the trigger of a phenomenon of “induced gravitational collapse” (Ruffini et al., 2001c) originating the supernova explosion and we also witness the birth of a young neutron star out of the supernova event. This extraordinary coincidence of these three astrophysical events represents an unprecedented scenario of fundamental importance in the field of relativistic astrophysics. Using the EMBH theory we shall explore: a) the process of black hole formation in the event GRB 980425 (Pian et al., 2000), b) the concept of “induced gravitational collapse”, introduced in the GRB-Supernova Time Sequence (GSTS) paradigm (Ruffini et al., 2001c), and its link to a very special supernova type in SN1998bw (Galama et al., 1998), and finally c) the observation for the first time of the cooling of a hot newly formed neutron star.

The observational situation of this system is quite complex. In addition to the source GRB 980425 and the supernova SN1998bw, two X-ray sources have been found by BeppoSAX in the error box for the location of GRB 980425: a source S1 and a source S2 (Pian et al., 2000). Since the nature of the two sources S1 and S2 was not clear, a variety of slopes in the decaying part of the afterglow have been proposed (see Fig. 1). Kulkarni et al. (1998) have proposed to explain both the supernova SN1998bw and the GRB 980425 observations by a new class of GRBs, distinctly different from the cosmological ones, both originated by a single unusual supernova event. Similarly, Iwamoto et al. (1998) have tried to explain both the supernova event and the GRB with a new kind of supernova with an extremely large explosion energy, larger than $`10^{52}\mathrm{ergs}`$, which they identify with the “hypernovae” predicted by Paczyński (1998). In this approach a totally novel concept is introduced: the supernova itself is assumed to originate in the process of gravitational collapse to a black hole of a massive progenitor star ($`40M_{}`$) with a particularly large angular momentum and strong magnetic field. A large rotational energy of the black hole extracted with a strong magnetic field is called in, by these authors, to explain the successful explosion of this “hypernova” leading both to the GRB and the supernova.

Our approach is drastically different. We first interpret the GRB 980425 within the EMBH theory. This allows the computation of the luminosity, spectra, Lorentz gamma factors, and more generally all the dynamical aspects of the source. Having characterized the features of GRB 980425, we can gradually approach the remaining part of the scenario, disentangling the GRB observations from the supernova ones and from the sources S1 and S2. This leads to a natural time sequence of events and to their autonomous astrophysical characterization.

## THE ENERGETICS, DYNAMICAL PARAMETERS AND SPACE-TIME PARAMETRIZATION OF GRB 980425

Our approach has focused on identifying the energy extraction process from the black hole (Christodoulou & Ruffini, 1971) as the basic energy source for the GRB phenomenon. The distinguishing feature is a theoretically predicted source energetics all the way up to $`1.8\times 10^{54}\left(M_{BH}/M_{}\right)\mathrm{ergs}`$ for $`3.2M_{}M_{BH}7.2\times 10^6M_{}`$ (Damour & Ruffini, 1975). In particular, the formation of a “dyadosphere”, during the gravitational collapse leading to a black hole endowed with electromagnetic structure (EMBH) has been indicated as the initial boundary conditions of the GRB process (Ruffini, 1998; Preparata et al., 1998). Our model has been referred as “the EMBH model”, although the EMBH physics only determines the initial boundary conditions of the GRB process by specifying the physical parameters and spatial extension of the neutral electron positron plasma originating the phenomenon created in the dyadosphere. The creation of this plasma is due to the vacuum polarization process occurring in a supercritical field by the Heisenberg-Euler-Schwinger process (see Heisenberg & Euler, 1935; Schwinger, 1951; Damour & Ruffini, 1975; Preparata et al., 1998).

Traditionally, following the observations of the Vela (Strong, 1975) and CGRO111see http://cossc.gsfc.nasa.gov/batse/ satellites, GRBs have been characterized by few parameters such as the fluence, the characteristic duration ($`T_{90}`$ or $`T_{50}`$) and the global time averaged spectral distribution (Band et al., 1993). With the observations of BeppoSAX222see http://www.asdc.asi.it/bepposax/ and the discovery of the afterglow, and the consequent optical identification, the distance of the GRB source has been determined and consequently the total energetics of the source has been added as a crucial parameter.

The observed energetics of GRBs, computed for spherically symmetric explosions, do coincide with the ones theoretically predicted in Damour & Ruffini (1975). This has been the major reason which has motivated us to reconsider and develop in full details the EMBH model. For simplicity, we have considered the vacuum polarization process occurring in an already formed Reissner-Nordström black hole (Ruffini, 1998; Preparata et al., 1998), whose dyadosphere has an energy $`E_{dya}`$. It is clear, however, that this is only an approximation to the real dynamical description of the process of gravitational collapse to an EMBH. In order to prepare the background for attacking this extremely complex dynamical process, we have clarified some basic theoretical issues, necessary to implement the description of the fully dynamical process of gravitational collapse to an EMBH (see Ruffini & Vitagliano, 2002,2003; Cherubini et al., 2002).

We have then given the constitutive equations for the five eras in the EMBH model (see for details Ruffini et al., 2003a and references therein). The Era I: the $`e^+e^{}`$ pairs plasma, initially at $`\gamma =1`$, self propels itself away from the dyadosphere as a sharp pulse (the PEM pulse), reaching Lorentz gamma factor of the order of 100 (Ruffini et al., 1999). The Era II: the PEM pulse, still optically thick, engulfs the remnant left over in the process of gravitational collapse of the progenitor star with a drastic reduction of the gamma factor; the mass $`M_B`$ of this engulfed baryonic material is expressed by the dimensionless parameter $`B=M_Bc^2/E_{dya}`$ (Ruffini et al., 2000). The Era III: the newly formed pair-electromagnetic-baryonic (PEMB) pulse, composed of $`e^+e^{}`$ pair and of the electrons and baryons of the engulfed material, self-propels itself outward reaching in some sources Lorentz gamma factors of $`10^3`$$`10^4`$; this era stops when the transparency condition is reached and the emission of the proper-GRB (P-GRB) occurs (Bianco et al., 2001). The Era IV: the resulting accelerated baryonic matter (ABM) pulse, ballistically expanding after the transparency condition has been reached, collides at ultrarelativistic velocities with the baryons and electrons of the interstellar matter (ISM) which is assumed to have an average constant number density, giving origin to the afterglow. The Era V: this era represents the transition from the ultrarelativistic regime to the relativistic and then to the non relativistic ones (see Ruffini et al., 2003a,2003b and references therein).

The EMBH model differs in many respects from the models in the current literature. The major difference consists in the following points:
a) The appropriate theoretical description of all the above mentioned five eras is implemented, as well as the evaluation of the process of vacuum polarization originating the dyadosphere. The description of the inner engine originating the GRBs has never been addressed in the necessary details in the literature;
b) The dynamical equations as well as the description of the phenomenon in the laboratory time and the time sequence carried by light signals recorded at the detector have been explicitly integrated (see e.g. Ruffini et al., 2003a, 2003b). In doing so we have also corrected a basic conceptual inadequacy, common to all the current works on GRBs, which led to an improper spacetime parametrization of the GRB phenomenon, preempting all these works from their predictive power: the relation between the photon arrival time at the detector and their emission time in the laboratory frame, expressed in our approach by an integral of a function of the Lorentz gamma factor extended over all the past source worldlines, has been in the current literature expressed as a function of an instantaneous value of the Lorentz gamma factor. These two approaches are conceptually very different and lead to significant qualitative differences (Ruffini et al., 2003a, 2003b and references therein).
c) The treatment of the afterglow has been widely considered in the literature by the so-called “fireball model” (see e.g. Mészáros & Rees, 1992, 1993; Rees & Mészáros, 1994; Piran, 1999 and references therein). However, also in the description of the afterglow, there are major differences between the works in the literature and our approach (Ruffini et al., 2003b).

The equations of motion in our model depend only on two free parameters: the total energy $`E_{tot}`$, which coincides with the dyadosphere energy $`E_{dya}`$, and the amount $`M_B`$ of baryonic matter left over from the gravitational collapse of the progenitor star, which is determined by the dimensionless parameter $`B=M_Bc^2/E_{dya}`$. The best fit of GRB 980425 is reproduced in Tab. 1. It correspond to $`E_{dya}=1.1\times 10^{48}\mathrm{ergs}`$, $`B=7\times 10^3`$ and the ISM average density is found to be $`n_{ism}=0.02\mathrm{particle}/\mathrm{cm}^3`$. The plasma temperature and the total number of pairs in the dyadosphere are respectively $`T=1.028\mathrm{MeV}`$ and $`N_{e^\pm }=5.3274\times 10^{53}`$.

## THE GRB 980425 LUMINOSITY IN SELECTED ENERGY BANDS PREDICTED BY THE EMBH MODEL

Recently, within the EMBH model, we have developed an attempt to theoretically derive the GRB spectra out of first principles as well as the GRB luminosity in fixed energy bands (Ruffini et al. 2003c). We have adopted three basic assumptions: a) the resulting radiation as viewed in the comoving frame during the afterglow phase has a thermal spectrum and b) the ISM swept up by the front of the shock wave, with a Lorentz gamma factor between $`300`$ and $`2`$, is responsible for this thermal emission. c) We also adopt, like in our previous papers (Ruffini et al., 2001a,2001b,2002,2003b), that the expansion occurs with spherical symmetry. This three assumptions are different from the ones adopted in the GRB literature, in which the afterglow emission is believed to originate from synchrotron emission in the production of a shock or reverse shock generated when the assumed jet-like ejecta encounter the external medium (see e.g. Giblin et al., 2002 and references therein).

In the EMBH model the structure of the shock is determined by mass, momentum and energy conservation: the constancy of the specific enthalpy, which is a standard condition in shock rest frames (Zel’dovich & Rayzer, 1966) and have been used in our derivation (Ruffini et al., 2003b). The only free parameter of our model is the size of the “effective emitting area” in the shock wave front: $`A_{eff}`$. Since the determination of this free parameter is performed here by empirically fitting the observational data, we avoid ambiguities due to the absence of relevant theoretical and laboratory results on relativistic shocks for Lorentz factor $`\gamma 300`$.

The temperature $`T`$ of the black body in the comoving frame is then

$$T=\left(\frac{\mathrm{\Delta }E_{\mathrm{int}}}{4\pi r^2\mathrm{\Delta }\tau \sigma }\right)^{1/4},$$
(1)
where

$$=\frac{A_{eff}}{A_{abm}}$$
(2)
is the ratio between the “effective emitting area” and the ABM pulse surface $`A_{abm}`$, $`\sigma `$ is the Stefan-Boltzmann constant and $`\mathrm{\Delta }E_{\mathrm{int}}`$ is the proper internal energy developed in the collision between the ABM pulse and the ISM in a proper time interval $`\mathrm{\Delta }\tau `$ (see Ruffini et al., 2003b,2003c). The ratio $``$, which is a priori a function that varies as the system evolves, is evaluated at every given value of the laboratory time $`t`$.

All the subsequent steps are now uniquely determined by the equations of motion of the system. The basic tool in this calculation involves the definition of the EQuiTemporal Surfaces (EQTS) for the relativistic expanding ABM pulse as seen by an asymptotic observer. The key to determining such EQTS (see Fig. 1 in Ruffini et al., 2002) is the relation between the time $`t`$ in the laboratory frame at which a photon is emitted from the ABM pulse external surface and the arrival time $`t_a^d`$ at which it reaches the detector. We have instead adopted the equations (see Ruffini et al., 2002,2003b):

$`t_a^d=\left(1+z\right)\left(t{\displaystyle \frac{r\left(t\right)}{c}}\mathrm{cos}\vartheta +{\displaystyle \frac{r_0}{c}}\right)=\left(1+z\right)\left(t{\displaystyle \frac{_0^tv\left(t^{}\right)𝑑t^{}+r_0}{c}}\mathrm{cos}\vartheta +{\displaystyle \frac{r_0}{c}}\right)`$
$`=\left(1+z\right)\left(t\mathrm{cos}\vartheta {\displaystyle _0^t}\sqrt{1{\displaystyle \frac{1}{\gamma ^2\left(t^{}\right)}}}𝑑t^{}+{\displaystyle \frac{r_0}{c}}\left(1\mathrm{cos}\vartheta \right)\right),`$ (3)
where $`z`$ is the cosmological redshift of the source, $`r_0r\left(t=0\right)`$ and $`\vartheta `$ is the angle subtended by the emission point of the photon on the ABM pulse external surface, having defined $`\vartheta =0`$ along the line of sight. For such a relation, different approximations exist in the literature (see e.g. Fenimore et al., 1996; Sari, 1997,1998, Rees & Mészáros, 1998; Fenimore et al., 1999; Granot et al., 1999 and see also the reviews by Piran, 1999; Mészáros, 2002). There are both quantitative and qualitative important differences in the use of Eq.(3) instead of the ones usually adopted in the literature (Ruffini et al., 2003b). It is a matter of fact that only the use of the EQTS defined by Eq.(3) allows to fit the observational data and that even a very minor departure from it leads to unacceptable results. It is also important to recall that this inadequacy in the current literature in the relation between the arrival time and the laboratory time has affected also all the estimates of the power-law slopes in the afterglow, preempting all current theoretical considerations in the literature of their predictive power. All the considerations about beaming in GRBs existing in the current literature have to be reformulated on the ground of the proper theoretical treatment and of Eq.(3) (see Ruffini et al., 2003b). Having so determined the EQTS we have been able to evaluate the source luminosity in a given energy band, in agreement with the above mentioned new assumptions (Ruffini et al., 2003c).

We can now proceed to the best fit of the GRB 980425 observed data. The best fit of the observed luminosity in selected energy bands has been obtained for the above mentioned values of the parameters $`E_{dya}`$, $`B`$ and for $``$ monotonically decreasing from $`4.81\times 10^{10}`$ to $`2.65\times 10^{12}`$. The results are given in Fig. 2 where the luminosity is computed as a function of the arrival time for three selected energy bands.

We can then conclude:

1. The best fit is obtainable under almost perfect spherical symmetry. This has been proven as a result of an analysis of unprecedented redundancy: the luminosity curves are obtained from an integration over almost $`10^7`$ different paths, relating the observer to the EQTS. This procedure tests, to a very high level of accuracy, any departure from spherical symmetry as well as any departure from the computed equations of motion of the source (Ruffini et al., 2003b). This same circumstance was encountered in GRB 991216 (Ruffini et al., 2003b).
2. Each luminosity curve as a function of the arrival time presents complex behavior, which could be erroneously interpreted as evidence for breaks in the power-law indexes leading to erroneous inferences on the possible existence of jets (Ruffini et al., 2003c).
3. Although each luminosity curve presents some special features, the bolometric luminosity has a very clear and simple power-law behavior with the “golden value” index $`n=1.6`$ (Ruffini et al., 2003b).

## THE NEW ASTROPHYSICAL SCENARIO AND THE NEWLY BORN NEUTRON STAR

In Fig. 3 the luminosities in the three bands are represented together with the optical data of SN1998bw (black dots), the source S1 (black squares) and the source S2 (open circles). It is then clear that GRB 980425 is separated both from the supernova data and from the sources S1 and S2.

While the occurrence of the supernova in relation to the GRB has already been discussed with the GRB-Supernova Time Sequence (GSTS) paradigm (Ruffini et al. 2001c), we like to address here a different fundamental issue: the possibility of observing the birth of a newly formed neutron star out of the supernova event, which in turn has been triggered by the GRB 980425.

In the early days of neutron star physics it was clearly shown by Gamow & Schoenberg (1940,1941) that the URCA processes are at the very heart of the supernova explosions. The neutrino-antineutrino emission described in the URCA process is the essential cooling mechanism necessary for the occurrence of the process of gravitational collapse of the imploding core. Since then, it has become clear that the newly formed neutron star can be still significantly hot and in its early stages will be associated to three major radiating processes (Tsuruta, 1964,1979,2002; Canuto, 1978): a) the thermal radiation from the surface, b) the radiation due to neutrino, kaon, pion cooling, and c) the possible influence in both these processes of the superfluid nature of the supra-nuclear density neutron gas. Qualitative representative curves for these cooling processes, which are still today very undetermined due to the lack of observational data, are shown in Fig. 3.

It is of paramount importance to follow the further time history of the two sources S1 and S2. If, as we propose, S2 is a background source, its flux should be practically constant in time and this source has nothing to do with the GRB 980425 / SN1998bw system. If S1 is indeed the cooling radiation emitted by the newly born neutron star, it should be possible to notice a very drastic behavior in its luminosity as qualitatively expresses in Fig. 3.

## CONCLUSIONS

It is particularly attractive, in conclusion, to emphasize some of the analogies and the differences between the case of GRB 991216 and the one of GRB 980425:

1. In both these sources a GRB and, independently, a supernova event are present. In the case of GRB 991216 the inferences of the supernova can be obtained only on the ground of the emission of iron lines (Piro et al., 2000; Ruffini et al., 2001c). In the present case of GRB 980425 we have a very fortunate circumstance: the GRB source is much weaker and is much closer to us ($`z=1.0`$ for GRB 991216 and $`z=0.00835`$ for GRB 980425). This situation is particularly important for obtaining detailed data on the supernova and on the possible occurrence of a newly born neutron star. This occurrence could not be observed in GRB 991216 due to the very large distance and to the overwhelming X-ray luminosity of the afterglow (see Ruffini et al., 2003b).
2. The energetics of GRB 991216 is $`E_{dya}=4.83\times 10^{53}\mathrm{ergs}`$, while the one of GRB 980425 is $`E_{dya}=1.1\times 10^{48}\mathrm{ergs}`$. It is very impressive that the EMBH model applies over a range of more than 6 orders of magnitude, giving important inferences on the source as well as on the structure of the ISM surrounding the source. It is significant that in both sources the condition of spherical symmetry appears to be strongly implemented and this fact is a very clear discriminant among all possible sources of energy for GRBs.
3. The main difference between GRB 991216 and GRB 980425 is then traced back, within the EMBH model, to the different parameters occurring in the dyadosphere and to the nature of the “effective” Reissner-Nordström geometry which we have used as a reliable estimate for the dynamical processes of gravitational collapse leading to the formation of the EMBH. For GRB 991216 we have a ratio between the electromagnetic energy and the total mass of the imploding core, described by the parameter $`Q/M`$, given by $`Q/M=0.23`$ while for GRB 980425 we have $`Q/M=6.5\times 10^4`$. In both cases a reasonable mass of the black hole is $`M_{BH}10M_{}`$ (Ruffini et al., 2003b).

A dedicated observational campaign, both with XMM and Chandra, to follow the cooling of the newly formed neutron star is needed in order to gain for the first time information on this extremely important astrophysical process.


text:
# CLUMP DEVELOPMENT BY THE NICKEL BUBBLE EFFECT IN SUPERNOVAE to be published on ApJ 10 June 2005, v 626 1 issue

## 1 INTRODUCTION

Observations of SN 1987A showed that the distribution of Fe in the ejecta is not what would be expected in the simplest models: it extended to higher velocity than expected and had a large filling factor for its mass of $`0.07M_{}`$ determined from the supernova light curve (McCray 1993; Li, McCray, & Sunyaev 1993). A plausible mechanism for the large filling factor is the Ni bubble effect, in which the radioactive progenitors of the Fe expand relative to their surroundings because of the radioactive power deposition (Woosley 1988; Li et al. 1993; Basko 1994). This effect is important during the first $`10`$ days after the supernova, when the radioactive power is significant and the diffusion of energy is not yet an important process.

An expected effect of the Ni bubble expansion is to create clumps in the nonradioactive gas. There is widespread evidence that the ejecta of core collapse supernovae are clumpy. The oxygen line profiles in the nearby Type II supernovae SN 1987A (Stathakis et al. 1991) and SN 1993J (Spyromilio 1994; Matheson et al. 2000) showed evidence for structure, implying that the gas is clumped. The velocity range for the emission extends to $`1,500\mathrm{km}\mathrm{s}^1`$ in SN 1987A and $`4,000\mathrm{km}\mathrm{s}^1`$ in SN 1993J. Similar evidence for clumping has been found in the Type Ib supernova SN 1985F (Filippenko & Sargent 1989). Among young supernova remnants, Cas A is the prototype of the oxygen-rich SNRs, which show evidence for freely expanding, oxygen-rich ejecta in clumps. Puppis A, with an age $`3,700`$ yr, is a more elderly example of such a remnant (Winkler et al. 1988).

Wang & Chevalier (2001, 2002, hereafter WC01, WC02) investigated the role of clumps in the evolution of supernova remnants. In Tycho’s remnant, believed to be the remnant of a Type Ia supernova, the presence of ejecta knots near the outer shock front requires a density contrast $`\chi 100`$ relative to the surrounding ejecta (WC01). The remarkable protrusions in the X-ray image of the Vela remnant are likely to be caused by ejecta clumps (Aschenbach et al. 1995) and WC02 found that $`\chi 1000`$ is needed to create the structures. The free expansion velocities for the clumps were estimated at $`3,000\mathrm{km}\mathrm{s}^1`$.

Our aim here is to investigate whether the Ni bubble effect can cause the clump structure apparently needed in the Tycho and Vela remnants. We compute one-dimensional hydrodynamical calculations, building on the work of Basko (1994), suitable for core collapse supernovae. We extend Basko’s work by including the effects of radiative diffusion and by carefully examining the density structure of the shell swept up by the Ni bubble. We further consider Type Ia supernovae using the same methodology. We show our computational setup and methods in § 2. The solutions for radiation-hydrodynamical effects are given in § 3. We also draw on the analogy to power input in a pulsar bubble to provide insight into the shell structure and evolution. In § 4 and § 5, we discuss the radiative effects on the Rayleigh-Taylor instabilities of the Ni bubble shell and the inferred ejecta-clump properties. Our conclusions are in § 6.

## 2 INITIAL CONDITIONS AND METHODS

We consider all of the $`{}_{}{}^{56}\mathrm{Ni}`$ synthesized in a supernova explosion initially resides in an isolated sphere in the center of supernova ejecta, with a mass $`M_{Ni}`$ and a density contrast $`\omega `$ relative to the surrounding diffuse ejecta. The total radioactive energy released in the $`{}_{}{}^{56}\mathrm{Ni}^{56}\mathrm{Co}^{56}\mathrm{Fe}`$ decay sequence is deposited at a rate $`q(t)`$ (Basko 1994)

$$q(t)=\frac{Q_{Ni}}{t_{Ni}}\mathrm{exp}\left(\frac{t}{t_{Ni}}\right)+\frac{Q_{Co}}{t_{Co}}\frac{\mathrm{exp}(t/t_{Co})\mathrm{exp}(t/t_{Ni})}{1t_{Ni}/t_{Co}},$$
(1)
where $`t_{Ni}`$,$`t_{Co}`$ and $`Q_{Ni}`$,$`Q_{Co}`$ are, respectively, the mean life and specific decay energy of $`{}_{}{}^{56}\mathrm{Ni}`$,$`{}_{}{}^{56}\mathrm{Co}`$:

$$\left\{\begin{array}{cc}t_{Ni}=7.6\times 10^5\mathrm{sec}=8.8\mathrm{days}^{}\hfill & {}_{}{}^{56}\mathrm{Ni}^{56}\mathrm{Co}\hfill \\ Q_{Ni}=2.14\mathrm{MeV}/\mathrm{decay}=3.69\times 10^{16}\mathrm{ergs}\mathrm{g}^1\hfill & \end{array}\right\}$$

$$\left\{\begin{array}{cc}t_{Co}=9.64\times 10^6\mathrm{sec}=111.5\mathrm{days}\hfill & {}_{}{}^{56}\mathrm{Co}^{56}\mathrm{Fe}\hfill \\ Q_{Co}=4.57\mathrm{Mev}/\mathrm{decay}=7.87\times 10^{16}\mathrm{ergs}\mathrm{g}^1.\hfill & \end{array}\right\}$$

(Lide 1992; $``$ 8.5 days, Fireston & Shirley 1996). The accumulation of the total radioactive energy rises from $`10^4`$ sec and saturates at $`10^7`$ sec (about 25% of the total radiative energy is deposited between $`10^7`$ and $`10^8`$ sec, Fig. 1). Before the saturation point is reached, the radioactive energy can leak out of the bubble by radiative diffusion. This can occur while the gas remains optically thick to the $``$ 1 MeV $`\gamma `$-rays.

We assume that the diffuse supernova ejecta are initially freely-expanding so that each gas element moves with a constant velocity $`v=r/t`$. The part of the ejecta unaffected by the radioactivity has its density drop as $`t^3`$, and the pressure drops as $`t^4`$ for a $`\gamma =4/3`$, radiation dominated fluid. We first refer our simulation initial conditions to Basko’s parameters for SN 1987A (Basko 1994 and the references therein):

$$M_{Ni}=0.075M_{},\omega =3$$

$$t_0=10^4\mathrm{sec},\rho _{a0}=10^4\mathrm{g}\mathrm{cm}^3,P_\mathit{0}=2.5\times 10^{11}\mathrm{dyn}\mathrm{cm}^2,$$

where $`t_0`$, $`\rho _{a0}`$ and $`P_0`$ are the initial age, ejecta substrate density (surrounding the bubble), and background pressure, respectively. The initial velocity at the bubble edge $`r=R_0`$ is given by the initial $`{}_{}{}^{56}\mathrm{Ni}`$ density contrast and age:

$$U_0=\frac{R_0}{t_0}=\left(\frac{3M_{Ni}}{4\pi \rho _{Ni}}\right)^{1/3}\frac{1}{t_0}=7\times 10^7\omega ^{1/3}\mathrm{cm}\mathrm{s}^1.$$

For $`\omega =3`$, $`R_0=5.0\times 10^{11}\mathrm{cm}`$. The initial pressure is assumed uniformly distributed inside and outside the bubble, independent of the bubble size. The pressure in the Ni bubble is determined by the radioactive power input at later times.

We note that the above parameters give an ejecta structure comparable to the inner flat-density component of the power-law ejecta density model that has a power index $`n=8`$ ($`\rho r^n`$ in the outer parts), an ejecta mass $`M=10M_{}`$, and an explosion energy $`E=10^{51}`$ ergs (Chevalier & Liang 1989). The transition between the flat and power law parts of the density profile occurs at $`v=3162\mathrm{km}\mathrm{s}^1`$ for these parameters. This set of explosion conditions is suitable for core collapse supernovae and had been used in our two-dimensional hydrodynamical simulations for the Vela supernova remnant (WC02). The energy density is initially dominated by the energy deposited by the supernova shock. Shigeyama & Nomoto’s (1990, their Fig. 2) model for SN 1987A indicates an adiabat $`\kappa =6.3\times 10^{15}`$ (cgs units) where $`P=\kappa \rho ^\gamma `$; for $`t_0=10^4`$ sec, $`\rho _0=9.4\times 10^5\mathrm{g}\mathrm{cm}^3`$ in our density model, and so $`P_0=2.7\times 10^{10}`$ dynes cm-2. We take $`P_0=2.7\times 10^{10}`$ dynes cm-2 as the appropriate value at $`10^4`$ sec. For $`t_0=100`$ sec, $`P_0=2.7\times 10^{18}`$ dynes cm-2. The time to set up the free expansion phase depends on the radius of the star because of the effects of reverse shock waves. For a blue supergiant progenitor star, like that of SN 1987A, the timescale is $`10^4`$ sec; for a red supergiant progenitor, the timescale is $`10^5`$ sec. It is needed to initiate the simulations at $`t10^2`$ sec to acquire a numerically saturated solution at $`t10^4`$ sec (the age of interest). The use of a freely-expanding velocity prior to this stage should not affect the 1-D symmetrical solution, considering that the radioactive heating is strong. Models of supernova explosions predict smooth 56Ni profiles in both thermonuclear and core collapse supernovae, with the density contrast $`\omega `$ below a factor of 3 (Nomoto, Thielemann & Yokoi 1984; Woosley 1988). We vary $`\omega `$ in the range of $`0.1\omega 3`$, with the other initial parameters scaled accordingly to the standard values described above.

Next we consider the initial parameters for Type Ia supernovae. We assume a powerlaw density profile $`n=8`$ as in core collapse supernovae, with $`M=1.4M_{}`$, $`E=10^{51}`$ ergs, and $`M_{Ni}=0.5M_{}`$. The transition between the flat and power law parts of the density profile occurs at $`v=8452\mathrm{km}\mathrm{s}^1`$ for these parameters. The background pressure in Type Ia ejecta is not well documented. To estimate, we make use of the entropy change $`dS=dQ/T`$ during the nuclear burning from C/O into Ni. During the process, $`kT0.6`$ MeV, and the energy released per nucleon is

$$dQ\left[2m(\mathrm{C})+2m(\mathrm{O})m(\mathrm{Ni})\right]c^2/560.79\mathrm{MeV}/\mathrm{nucleon},$$
(2)
where $`m(\mathrm{X})`$ is the mass of element X; $`m(\mathrm{C})=12.0000\mathrm{amu}`$, $`m(\mathrm{O})=15.9949\mathrm{amu}`$, and $`m(\mathrm{Ni})=55.9421\mathrm{amu}`$ (Audi et al. 2003), and $`c`$ is the speed of light. Then $`dS1.3k/\mathrm{nucleon}=1.8\times 10^{16}\mathrm{erg}/\mathrm{K}/\mathrm{nucleon}`$, where $`k`$ is the Boltzmann constant (Shigeyama, private communication). The adiabat $`\kappa `$ was obtained by relating the entropy per nucleon $`S=[(4aT^3/3)(m(\mathrm{Ni})/56)/\rho ]`$, where the initial entropy before the burning is ignored, to the radiative pressure $`P=aT^4/3`$ and $`P=\kappa \rho ^{\frac{4}{3}}`$. Subsituting $`P`$ and $`\kappa `$ into $`S`$, we find

$$\kappa =\left[(\frac{3}{a})^{1/4}(\frac{S}{4(m(\mathrm{Ni})/56)})\right]^{4/3}6.1\times 10^{14}(\mathrm{cgs}\mathrm{units}).$$
(3)
For $`t_0=10^4`$ sec, $`\rho _0=6.92\times 10^7\mathrm{g}\mathrm{cm}^3`$ in our density model, and so $`P_0=3.7\times 10^6`$ dynes cm-2. The adiabat (and so the background pressure) could be higher, considering there is an initial entropy distribution that acts to freeze the deflagration front. The background energy density has the profile of the mass density because of the assumed uniformity.

We first approach the problem using hydrodynamical simulations with $`\gamma =4/3`$ in both the bubble interior and the ejecta substrate. We follow the bubble-shell interface on a spherical expanding grid, and add the specific radioactive energy to the internal energy within the bubble uniformly and locally, i.e., in proportion to the local mass density. A reflecting condition is used for the inner boundary, and a non-zero gradient outflow condition is invoked on the outer boundary to eliminate the spurious shock raised by the grid expansion.

We further include radiation hydrodynamics (RHD) in the simulations using the ZEUS2D code (Stone, Mihalas, & Norman 1992). The radioactive power is deposited into the radiation instead of the material and $`\gamma =5/3`$ for the matter. The code is based on finite difference and finite volume on an Eulerian grid and uses an artificial viscosity to smooth shock transitions. The RHD algorithm employs nonrelativistic full radiative transport, under the assumption of LTE (local thermal equilibrium) and gray opacity, and is thus equally applicable to both optically thin and thick media. The flow evolves towards the optically thin regime where radiation tends to stream away across the Ni bubble shell. The full transport scheme ensures the proper evolution while the radioactive energy is continually added within the bubble.

At the age of interest, the opacity is dominated by electron scattering. A maximum absorption 1/1000 times the magnitude of scattering is applied in our calculation. For completely ionized heavy elements, the Thompson scattering gives an opacity $`0.2\rho (r,t)\mathrm{cm}^1`$, where $`\rho (r,t)`$ is the material density. The initial energy is equally partitioned between radiation and material. Considering that a Type Ia has lower temperature and smaller radius, $`P_{rad}r^4`$ and $`P_{gas}r^5`$, the equi-partition assumption may be more valid for core collapse supernovae.

In the RHD algorithm, the coupled radiation and material energy densities are solved using Newton-Raphson iteration for convergence. When the radiation pressure is in excess of the gas pressure by several orders of magnitude, convergence becomes extremely difficult to reach. We thus stop the simulations prior to the optically thin stage. In our situation, the radiation effects are most important when the gas is still optically thick, so that the diffusion approximation is adequate. The choice of the flux limiters from which the Eddington tensor is computed does not affect the diffusion result. For details of radiation hydrodynamics, we refer readers to Stone et al. (1992) and the references therein.

## 3 SHELL STRUCTURE AND EVOLUTION

### 3.1 Simulations

We show in Figs. 2, 3, and 4 the ejecta structure and evolution properties in the core collapse supernova model using $`n=8`$, $`E=10^{51}\mathrm{erg}\mathrm{s}^1`$, and $`M=10M_{}`$. These simulations were initiated with $`t_0=10^2`$ sec, $`M_{Ni}=0.075M_{}`$ and $`\omega =3`$. The diffuse ejecta have an initial density $`\rho _0=94.4\mathrm{g}\mathrm{cm}^3`$ and an accumulated radioactive energy density $`e_{rad0}=1.38\times 10^{15}\mathrm{erg}\mathrm{cm}^3`$. The initial material and radiation energy density in the background are equally divided, $`e_0`$ = $`e_{r0}`$ =$`1.35\times 10^{10}\mathrm{erg}\mathrm{cm}^3`$.

We first examined the evolution of the flow in the purely hydrodynamic simulations. The expansion of the bubble gives rise to a strong forward shock behind which the ambient gas is compressed into a dense shell (Figs. 2 and 3). The inner edge of the shell is a contact discontinuity where the gas comoves with the bubble-shell interface; it has a high density because of the acceleration of the shell. The shock front is continually accelerated outward (Fig. 4(a)). To describe the acceleration, we use $`a(dR_{sh}/dt)/(R_{sh}/t)`$, the expansion rate evaluated at the shock front $`R_{sh}`$, equivalent to the velocity ratio of the shell to the preshock diffuse ejecta. We find that $`a`$ rises to a maximum $`1.20`$ around $`10^6`$ sec, and subsequently tends to 1.0 at late times as the shell becomes comoving with the preshock gas (Fig. 4(b)). Because the background pressure has little influence on the dynamics, the shell structure becomes frozen into the flow.

The shell is very thin, with a thickness as a fraction of the radius $`\beta h_{sh}/R_{sh}0.02`$, where $`R_{sh}`$ and $`h_{sh}`$ are respectively the radius and thickness of the shell (Fig. 4(c)). The initial shock compression is set by the jump condition

$$\frac{\rho _1}{\rho _0}=\frac{(\gamma +1)^2}{(\gamma 1)^2+2},$$
(4)
where

$$=\frac{v_{sh}}{c_0}=\left(\frac{dR_{sh}}{dt}\frac{R_{sh}}{t}\right)\left(\frac{\rho _0}{\gamma p_0}\right)^{1/2}$$
(5)
is the Mach number of the shock, $`\rho _1`$ is the postshock density, $`v_{sh}`$ is the shock velocity relative to the preshock gas, and $`c_0`$ is the sound speed in the preshock medium. For a radiation dominated, strong shock, $`\rho _1/\rho _0=7`$. The accelerating shell has a higher density in the postshock region so that a typical compression factor in the shell (relative to the preshock gas) is $`\chi 10`$ (Fig. 4(d)). The compression factor in the bubble is $`\chi 0.1`$ (Fig. 4(e)). We note that $`\beta (3\chi )^1`$. The highest density occurs at the inner edge, where the density is limited by the numerical resolution (Fig. 3).

In the RHD case, the energy from the heated gas leaks out from the bubble energy reservoir as time progresses (Fig. 2). The gas is radiation dominated. Radiative diffusion allows the internal energy to propagate ahead of the shock wave, and so eliminates the temperature and the radiation pressure gradients. In order for diffusion to be dynamically important, it requires that the diffusion time scale, $`t_d`$, be smaller than the hydrodynamic time scale of the flow, $`t_h`$; i.e.,

$$t_d\frac{h_{sh}^2}{\lambda c}<t_h\frac{h_{sh}}{v_{sh}}$$
(6)
where $`\lambda `$ is the photon mean free path. We note that

$$\frac{t_d}{t_h}\left(\frac{h_{sh}}{\lambda }\right)\left(\frac{v_{sh}}{c}\right).$$
(7)
At $`t3\times 10^6`$ s, $`h_{sh}=10^{13}`$ cm, $`\lambda =1/\kappa \rho =5\times 10^{10}`$cm, so that $`t_d=h_{sh}^2/\lambda c=6\times 10^4`$ s. At this time, $`v_{sh}1\times 10^7`$ cm s-1, so $`t_d/t_h0.07`$. The gas is optically thick, but radiative diffusion is important. As a result of diffusion, the preshock gas is accelerated by the radiation; the shock is then weakened. The shock front develops a radiative precurso (Fig. 2). The acceleration of the bubble shell is reduced because of the radiative loss. The rise in the preshock sound speed and the acceleration of gas ahead of the shock reduced the shock compression (eqn. ). The shell becomes broader and less dense compared to the HD case (Fig. 4(c)). The sharp density jump at the shock front in the RHD case is also smoothed out before the gradient of $`e_r`$ diminishes to zero (Fig. 4(d)). The gas tends toward a freely expanding state in which the density profile is frozen into the flow.

The diffusion effect speeds up the freeze-out of the ejecta structure. For our parameters, the freezing occurs at a velocity $`v_f1500\mathrm{km}\mathrm{s}^1`$ and a time $`t_f10^7`$ s, significantly earlier than the HD case where $`t_f10^9`$ s and $`v_f1900\mathrm{km}\mathrm{s}^1`$ (Fig. 4(f)). In the HD case, the expansion sustains until pressure equilibrium is reached. Both cases show a density contrast $`\chi 10`$ between the shell and the preshock ejecta. The inner edge (contact discontinuity) of the shell is the densest with $`\chi 100`$ . The shell thickness is increased to $``$ 4% of the radius of the bubble. The swept-up mass in the shell remains about $`M_s1M_{}`$ (Fig. 4(g)). Note that we define the thickness as the distance from the contact discontinuity to the postshock outer dense region at the shock front, considering that the outermost postshock density profile flattens out into the ejecta substrate. At $`10^7`$ sec, the surface density of the shell drops below $`1733g/cm^2`$, the mean free path for the $`1MeV`$ $`\gamma `$-ray photons (Woosley, Pinto, & Hartmann 1989; Shigeyama & Nomoto 1990). We stopped the RHD calculations prior to this time (Fig. 4(h)). When the initial Nickel mass is increased to 3 times $`0.075M_{}`$, the onset of freezing delays to $`t_f10^8`$ s with $`v_f1800\mathrm{km}\mathrm{s}^1`$ and $`M_s3M_{}`$. For a 10 times increase, $`t_f10^8`$ s, $`v_f2400\mathrm{km}\mathrm{s}^1`$, and $`M_s5M_{}`$ (Fig. 5). It is therefore expected that the radiation relaxes before $`10^8`$ s, given any large Ni overabundance. The reason is that the photon mean free path is increased by the broadening of the shell, and so diffusion becomes more rapid with time.

### 3.2 Effect of Radiation

The effect of radiative diffusion does not increase the density contrast of the Nickel bubble shell, as in shell formation during shock break-out in a Type II supernova where radiative energy losses by diffusion lead to a large compression that is eventually limited by the gas pressure (WC01, Chevalier 1976). In the latter case, the flow behind the shock wave is being decelerated, and the immediate postshock region is being compressed. A drop in the dominant radiation pressure then is compensated by compression of the gas to the point where gas pressure becomes important, resulting in a large density enhancement. In our case, the pressure gradient has an opposite sign; the pressure decreases outward from the contact discontinuity to the shock front, and the radiative loss to the preshock gas is replenished by the energy diffusing from the bubble. The diffusion of energy ahead of the shock front leads to a weaker shock front and less compression in the shell.

### 3.3 Models for Type Ia supernovae

We investigated the ejecta structure for Type Ia SNe using the inner flat component of the power law model, $`n=8`$, $`E=10^{51}`$ ergs, and $`M=1.4M_{}`$. We present the evolution of the dynamical properties in Fig. 7. The simulations were our highest resolved runs initiated with $`t_0=10^2`$, $`M_{Ni}=0.5M_{}`$, and $`\omega =1.5`$. The initiated time gives an initial mass density $`\rho _0=6.92\times 10^1\mathrm{g}\mathrm{cm}^3`$, and an initial material and radiation energy density $`e_0`$ = $`e_{r0}`$ = $`5.5\times 10^{14}\mathrm{erg}\mathrm{cm}^3`$ for a standard adiabat $`\kappa =6.1\times 10^{14}(\mathrm{cgs}\mathrm{units})`$. Here a lower background pressure was used, $`e_0=1.3\times 10^{13}\mathrm{erg}\mathrm{cm}^3`$. The initial $`{}_{}{}^{56}\mathrm{Ni}`$ has an accumulated radioactive energy density $`e_{rad0}=3.15\times 10^9\mathrm{erg}\mathrm{cm}^3`$.

The evolution of the flow is similar to the case of core collapse supernovae. The Nickel bubble shell has a density increase toward the contact discontinuity, with the density contrast reaching $`\chi 100`$. The diffusion effect speeds up the freeze-out of the ejecta structure. A notable difference in the shell is that it is thinner, with a very prominent precurso (Fig. 6). Compared to core collapse supernovae, the ejecta has a higher ratio of background energy to mass ($`E/M`$), thus the acceleration of the bubble is weaker. The expansion rate $`a`$ rises to a maximum 1.04 before radiation relaxes, barely faster than the free expansion. The thickness of the shell reaches $`\beta 0.004`$ at the turnover. The freezing occurs at a velocity $`v_f7000\mathrm{km}\mathrm{s}^1`$ and a time $`t_f10^7\mathrm{s}`$. The gas becomes optically thin at $`10^6`$ sec. By this time, the RHD results do not diverge substantially from the HD results, but the velocity and the accumulated mass curves had started to flatten out, indicating that the radiation effect becomes dynamically less important.

We note that the swept-up mass in the RHD case may be underestimated due to the large precurso; we expect $`M_s0.5M_{}`$ as in the HD case. For a less amount of $`{}_{}{}^{56}\mathrm{Ni}`$, we found that $`M_s<0.3M_{}`$ and $`v_f5500\mathrm{km}\mathrm{s}^1`$ for $`M_{Ni}=0.225M_{}`$, and $`M_s<0.2M_{}`$ and $`v_f4000\mathrm{km}\mathrm{s}^1`$ for $`M_{Ni}=0.075M_{}`$ (Fig. 5).

### 3.4 Variation of the Initial Parameters

We carried out simulations with a variety of initial conditions. We found the solution in the saturated regime is insensitive to the initial background pressure in the ejecta (Fig. 8). In the case of a higher initial background pressure that exceeds the accumulated radioactive pressure, the sharp structure in the shell takes more time to build up (for which grid resolution also has an effect), but the solution gradually converges once the radioactive energy input becomes dominant. E.g., the case using a smaller background pressure $`P_0=2.7\times 10^{10}\mathrm{dyn}\mathrm{cm}^2`$ is no different from the standard case using $`P_0=2.7\times 10^{18}\mathrm{dyn}\mathrm{cm}^2`$ (for the core collapse SN model with $`t_0=10^2`$ sec), except that the former shows a more pronounced postshock oscillation trailing in the bubble. Likewise, the solution in our Type Ia model appears insensitive to the variation of background pressure. Our case using the standard $`\kappa `$ shows little differences from the case in Fig. 7. We expect that the insensitivity would sustain with a higher value of $`\kappa `$. Simulations initialized before $`10^4`$ sec allow the flow to saturate before the age of interest, and so the results are also insensitive to the initial distribution of energy between radiation and material, absorption rate, and details of radiative transport.

We also varied the initial $`{}_{}{}^{56}\mathrm{Ni}`$ density contrast. We assumed $`{}_{}{}^{56}\mathrm{Ni}`$ initially in pressure equilibrium with the ambient ejecta, so that its initial background energy increases with decreasing $`\omega `$. Thus for $`\omega =0.1`$ the shell expands to an outer radius. However, the change from $`\omega =3`$ to $`\omega =1`$ does not alter the shell properties we determined (Fig. 8).

## 4 ANALOGY TO PULSAR WIND NEBULAE

Chevalier (1984) has considered a constant-luminosity pulsar bubble interacting with constant density, freely-expanding supernova ejecta. We compare his self-similar solution with our hydrodynamical case, since the time evolution of the driving power is the only difference. The self-similarity exists for a power-law ejecta density $`\rho v^nt^3`$ and a time-varying pulsar luminosity $`Lt^l`$, when $`n3`$ and $`l1`$ (Jun 1998); the expansion law is $`rt^{(6nl)/(5n)}`$, as obtained from dimensional analysis. For a constant luminosity $`l=0`$ and flat ejecta $`n=0`$, $`rt^{1.2}`$. The shock cannot be slowed down in the ejecta more than the free expansion rate, so the self-similar solutions only exist for $`a1`$, i.e. $`l1`$.

In our solutions, the shell starts to accelerate with $`a1`$. The evolution of the density profile can thus be compared to the self-similar solutions with varying $`l`$. In the self-similar solutions, decreasing $`l`$ (i.e. increasing power input), or increasing the expansion rate, increases the shell thickness. For $`l=0`$ and $`a=1.2`$, we have $`\beta =0.02`$. The shell thickness in our simulations shows a more complicated evolution than indicated by the measured expansion rate; the maximum expansion rate measured occurs at $`6\times 10^5`$ sec, while the maximum thickness occurs at $`5\times 10^6`$ sec. The first turnover in the thickness evolution seems to correlate with the change in the luminosity power index (see Fig. 9) at $`5\times 10^6`$ sec; the second turnover occurs after $`10^7`$ sec when the self-similarity breaks down because $`l>1`$.

We overplot our solution at $`4\times 10^5`$ sec on the self-similar solution with $`l=0.03`$ in Fig. 10. The density distribution in the shell shows a sharp inward increase towards the contact discontinuity, with over 80% of the shell mass concentrated within the inner 50% region in radius (Fig. 11). The self-similar density is infinite at the contact discontinuity. The coarseness of the grid obviously sets limits on the highest density computed in the grid domain. Our results show a higher density contrast on fixed grids, because the numerical noise is smaller than on expanding grids. The shell can be distorted by instabilities (Basko 1994; Jun 1998), but we expect the density contrast to remain comparable.

Basko’s (1994, his Fig. 1) HD solutions suggest that the density contrast of the shell is only $`\chi 5`$ at $`4\times 10^5`$ sec. We reproduced his case and found that the shell broadens to $`\beta >0.1`$ at $`4\times 10^4`$ s, the turnover is much earlier than in our calculation using $`t_0=10^2`$ s. In his case using $`t_0=10^4`$ s, the radioactive energy accumulated to $`10^4`$ s was deposited all at one time in the beginning of the calculation, so an excessive acceleration and a large broadening of the shell were the result. Basko’s Fig. 6 for the shell surface density has the same earlier turnover, suggesting this heating effect.

We note that the shell characteristics such as the density contrast and the thickness ratio should only be weakly dependent on the initial $`{}_{}{}^{56}\mathrm{Ni}`$ mass, since these properties can be approximated by the self-similar solutions.

The deformation of the bubble-shell interface due to the Rayleigh-Taylor instability is likely to have a morphology similar to that of the pulsar bubble shell studied by Jun (1998) in two-dimensional simulations. The main difference is that magnetic fields constrain the instability in the pulsar bubble case, but are unlikely to be important for our case. The Rayleigh-Taylor spikes are directed inward, because of the shell acceleration. However, the instability should eventually end when the shell stops accelerating in the comoving stage. The effect of radiation diffusion is to smooth the pressure gradients, which ends the shell acceleration and the instability growth. The Jun (1998) simulations indicate that the instability does not cause radioactive material to have a significantly higher final velocity (see also Basko 1994).

## 5 NICKEL BUBBLE SHELL AS EJECTA CLUMPS

We now compare the properties of the Ni bubble shell with the initial properties of the clumps/bullets inferred from our previous clump-remnant simulations for Tycho’s remnant and the Vela supernova remnant. We consider that the ejecta clumps originate as components from the breakup of the shell, and examine if the inferred clumps (shell fragments) are able to survive the crushing in the remnant. The robustness of the clumps can be determined by three parameters: the initial density contrast, $`\chi `$, between the clump and supernova ejecta; the initial impact time with the reverse shock, related to the initial ejection velocity; and the initial size of the clump. The high density contrast created across the shell, $`\chi 100`$, is compatible with that of ejecta clumps needed to survive crushing in the remnant. In our model for the Vela remnant ($`M=10M_{}`$, $`E_{51}=1`$, and $`M_{Ni}=0.075M_{}`$), the shell’s frozen-in velocity, $`v_f1500\mathrm{km}\mathrm{s}^1`$, is within the range $`1000v30004000\mathrm{km}\mathrm{s}^1`$ we determined as the most likely origin for the Vela bullets (WC02). If the shell has not been greatly disturbed, the clumps should be present in bands at the restricted velocity space; they are expected to enter the remnant’s intershock region at a normalized time $`t^{}2.2`$, or $`2800`$ yrs (see Fig. 1 in WC02 and Fig. 12). The velocity-restricted nature may have been detected in the wavelike clumpy structure immediately below Tycho’s blast wave front (Reynoso et al. 1999, WC01). If the clumps attain a size about the thickness of the shell, $`\beta =0.04`$, then at the time of the initial reverse shock impact, the clumps have an initial size $`a_08\%`$ in terms of the intershock width, or $`r_03\%`$ in terms of the forward shock radius (Fig. 1 in WC02). Compared to the late clump-remnant interaction in WC02 that used a larger clump size $`a_01/3`$, the inferred clumps do not appear to have a sufficient strength to cause a protrusion on the forward shock. Nevertheless, since the computed highest density is limited by numerical resolution, a larger density contrast should be achieved in a small region. An increase in the density contrast, e.g., from 100 to 1000, can significantly improve the strength, considering that the cloud crushing time scale is proportional to $`t_{cc}r_c^{1/2}`$ for a fixed mass (WC02, Klein et al. 1994). Such small-sized but dense clumps are likely to form, given the steep density profile of the shell.

In our model for Type Ia supernovae ($`M=1.4M_{}`$, $`E_{51}=1`$, and $`M_{Ni}=0.5M_{}`$), the frozen-in velocity reaches $`v_f7000\mathrm{km}\mathrm{s}^1`$, resulting in a clump-remnant interaction initiated in an early state, $`t^{}1.02`$. Compared to the present state $`t^{}1.7`$ of Tycho’s remnant, the ejection has a strength as indicated for Tycho’s knots in our previous 2-D HD simulations (WC01). In addition, the swept-up mass $`0.5M_{}`$ in the shell is reasonably in excess of the mass 0.002 $`M_{}`$ for the Si+S clump and 0.0004 $`M_{}`$ for the Fe clump in Tycho’s SE outline (Hwang, Hughes, & Petre 1998). However, there is evidence for Tycho that it is subluminous (van den Bergh 1993; Schaefer 1996), implying a deficiency in $`{}_{}{}^{56}\mathrm{Ni}`$. Observations of SN 1006 suggest an initial Nickel mass in the range of $`0.0750.16M_{}`$ (Hamilton et al. 1997, Mazzali et al. 1997). For Sn 1991 bg, Mazzali et al. (1997) estimate $`0.07M_{}`$ of $`{}_{}{}^{56}\mathrm{Ni}`$ extending to a velocity of 5,000 $`\mathrm{km}\mathrm{s}^1`$. In our model, the velocity only extends to $`4,000\mathrm{km}\mathrm{s}^1`$ with $`M_{Ni}=0.07M_{}`$.

For core collapse supernovae, we consider that the object RX J0852–4622 superposed the boundary of the Vela remnant is a result of the clump-remnant interaction. WC02 estimated a mass of $`0.1\beta M_{}`$, where $`\beta <1`$, for the responsible clump, which gives a fraction $`10\%`$ in terms of the total swept-up mass in the shell. This mass fraction appears somewhat large compared to the Si+S clump in Tycho’s remnant. It is unclear how such a massive clump is created, if the initial $`{}_{}{}^{56}\mathrm{Ni}`$ is not overabundant.

In our scenario, the initial $`{}_{}{}^{56}\mathrm{Ni}`$ has a symmetric distribution in the center of supernova ejecta. It should be noted that dynamical processes prior to the Nickel bubble effect, e.g., neutrino-driven convection (Kifonidis et al. 2000), may give rise to an anisotropic distribution of elements. If $`{}_{}{}^{56}\mathrm{Ni}`$ is mixed outwards during the explosion, the amount of radioactive material inside each individual bubble could be orders of magnitude lower; the subsequent radioactive heating may cause merging between Nickel bubbles, and so give further compression. This process remains to be examined by future multidimensional calculations. Our result here should be noted as a lower limit on the ejecta-clump properties for the clump-remnant interaction.

## 6 CONCLUSIONS

We investigated the ejecta structure in supernovae due to the radioactive heating from the $`{}_{}{}^{56}\mathrm{Ni}^{56}\mathrm{Co}^{56}\mathrm{Fe}`$ decay sequence, under the assumption of spherically symmetric flow. Two approaches were used: hydrodyamics (HD) simulations and radiation hydrodynamics (RHD) simulations. The RHD calculation included a proper treatment of radiation transport prior to the age $`10^7`$ sec, when about 25% of the total radioactive energy is still to be released. Radiative effect becomes dynamically less important as the evolution tends towards the optically thin regime. Based on the mutual coherence in the HD and RHD calculations, we obtained the properties of the shell driven by the radioative expansion. We believe that going to later stages of diffusion would not change this conclusion.

The expansion of the Nickel bubble can sweep up a dense shell of $`1M_{}`$ shocked ambient ejecta in core collapse and Type Ia supernovae, with the highest density in the shell over 100 times the ambient ejecta density. The shell has an inward density increase toward the bubble-shell interface, and the highest density contrast computed is limited by numerical resolution. For the hydrodynamical solution, the structure of the shell can be approximately described by a self-similar solution. For the RHD solution, the main difference is that radiative diffusion gives rise to a broader and less dense shell that freezes out in the ejecta. Because the radioactive pressure eventually dominates, the saturated solution is not sensitive to the initial conditions such as the background pressure and the initial density contrast of Ni; the only crucial factor to the solution is the initial $`{}_{}{}^{56}\mathrm{Ni}`$ abundance. We presume that a higher background pressure due to the uncertainty of the adiabat would not affect the results.

We examined if the properties of the shell fragments are comparable to that of the ejecta clumps protruding the outlines of Tycho’s remnant and the Vela remnant. The density contrast created across the shell, $`\chi 100`$, is compatible with that of the ejecta clumps needed to survive the clump-remnant interaction. In our standard Type Ia model, the shell’s frozen-in velocity attained at $`10^6`$ sec appears sufficient to give an ejection as indicated for Tycho’s knots, in spite that most ($`75\%`$) of the radioactive energy is still to be released. For core collapse supernovae like the Vela remnant, however, the small size and the late ejection of the inferred clump result in a weak strength to resist crushing and further cause protrusions on the remnant outline.

I am grateful to Roger Chevalier for supporting this work, and Toshi Shigeyama for useful correspondence on Type Ia supernovae. I thank the referee for careful and constructive comments on the manuscript. The computations were carried out on the IBM SP2 at University of Virginia, and the PC cluster provided by Taiwan’s NSC grant NSC93-2112-M033-002. Support for this work was also provided by NSC grant NSC93-2112-M033-013..