text:
# Dark Matter Constituents

## 1 INTRODUCTION

The question of the nature of the average composition of the Universe is getting more interesting than ever. As new observations have arrived, the picture emerging has become quite puzzling, with the favoured model now containing both dark energy (vacuum energy, “cosmological constant” or perhaps a slowly evolving scalar field, “quintessence”) and dark, non-baryonic matter.

The recent supernova discovery of acceleration of the cosmic expansion  taken in combination with the spectacular new WMAP microwave background measurements  has led to a “concordande model” ($`\mathrm{\Lambda }`$CDM):

$$\{\begin{array}{ccc}\mathrm{\Omega }_M& =& 0.3\\ \mathrm{\Omega }_\mathrm{\Lambda }& =& 0.7\end{array}$$
(1)
This is in nice agreement with the inflationary prediction of the matter and cosmological constant contributions adding up to unity, $`\mathrm{\Omega }_{tot}=1`$, but produces the mystery of having a non-zero but exceedingly small value of the cosmological constant $`\mathrm{\Lambda }`$ (compared to the natural scale of gravity, which is the Planck scale). However, the supernova data (of two independent groups) and CMBR measurement in combination with the measured range of the Hubble constant ($`H_070`$ km/s Mpc-1, $`hH_0/100=0.7\pm 0.1`$) has more or less forced this model upon us.

The WMAP data and Big Bang nucleosynthesis (BBN) deuterium data  agree on the baryonic contribution $`\mathrm{\Omega }_b`$ to $`\mathrm{\Omega }_M`$:

$$\mathrm{\Omega }_bh^20.023\pm 0.002.$$
(2)
This means that non-baryonic dark matter is required beyond any doubt, since $`h=0.7\pm 0.1`$ gives the baryonic contribution $`\mathrm{\Omega }_b`$ less than around $`0.08`$, so much less than the measured $`\mathrm{\Omega }_M`$.

Weakly interacting massive particle candidates (WIMPs) for dark matter are generically found to decouple at a temperature that is roughly $`m_{\mathrm{WIMP}}/20`$, which means that they are non-relativistic and behave as CDM when structure forms.

WIMPs have the advantage that their relic density is determined by the calculable departure from thermal equilibrium in the early universe, and give without difficulties a contribution to $`\mathrm{\Omega }`$ which may be close to the observed value.

## 2  DARK MATTER CANDIDATES

### 2.1 Neutrinos

With the masses of neutrinos in view of observed oscillations being definitely non-zero, we may for the first time say with certainty that non-baryonic dark matter exists. However, the neutrino contribution to $`\mathrm{\Omega }_M`$ is most likely very small. For small mass neutrinos (hot dark matter, HDM), the relic abundance is given by

$$\mathrm{\Omega }_\nu h^2=\frac{\underset{i}{}m_{\nu _i}}{94\mathrm{eV}}.$$
(3)
For the kind of mass values indicated by the Super-Kamiokande measurements, 0.01 – 0.1 eV (unless there are unnatural neutrino mass degeneracies), we see that although neutrinos indeed contribute to non-baryonic dark matter, they do so at the level or $`10^2`$ or less. This may seem tiny from the cosmological point of view, but is still about as much as all visible material like stars, gas, dust, etc contribute to the energy density of the Universe.

The WMAP results, especially when augmented by large galaxy survey data such as SDSS  or 2DF , also give very little room for the standard neutrino contribution to dark matter. Since the distribution of matter is quite different in CDM and HDM models (with the latter lacking small-scale power due to free-streaming), one may use the observed distribution to put a direct upper bound to the sum of all neutrino masses which turns out to be of the order of an eV or so (the WMAP limit is 0.23 eV for all three neutrinos degenerate, but this includes less understood smaller scale data from the Lyman alpha forest; for a discussion see also ). The advantage of the cosmological limit is that it is sensitive to the absolute mass scale of neutrinos whereas neutrino oscillations only probe mass differences.

In principle, one could imagine having a sterile neutrino as dark matter, if it is non-thermally produced, e.g., generated through mixing with the active neutrinos. Generally, this candidate will have to have a mass in the keV to MeV range and would act as something inbetween cold and hot dark matter (sometimes named “warm dark matter”, WDM). An unpleasant feature of these models is a necessary, delicate finetuning of the mixing angle versus mass to get the right abundance, but models of this kind have been constructed which may perhaps evade experimental constraints .

### 2.2 Other Solutions to the Dark Matter Problem

A class of dark matter candidates, which should be kept in mind, are super-heavy relics which were produced by some non-thermal mechanism in the early Universe. Examples are given by “wimpzillas” .

At the other end of the mass spectrum are axions, light pseudoscalar mesons that appear in all string theories and which are needed to explain the weakness of CP violation in strong interactions. There are experiments which search for these particles in mass ranges interesting for cosmology .

There has recently been interest in dark matter in the form of topological, extended objects, “Q-balls” which may or may not be supersymmetric and may or may not be stable on cosmological time scales (see  for a recent review). In certain scenarios, supersymmetric Q-balls may decay into the lightest neutralino, giving a similar scenario as in the MSSM, but with different phenomenology. For instance, the usual bounds on the relic density for thermally produced neutralinos may be evaded.

## 3  SUPERSYMMETRIC DARK MATTER

One of the prime candidates for the non-baryonic component is provided by the lightest supersymmetric particle, plausibly the lightest neutralino $`\chi `$.

The lightest neutralino $`\chi `$ is a mixture of the supersymmetric partners of the photon, the $`Z`$ and the neutral part of the two Higgs doublets present in the minimal extension of the supersymmetric standard model (for reviews see, e.g., ). The attractiveness of this dark matter candidate stems from the fact that its generic couplings and mass range naturally gives a relic density in the required range to explain halo dark matter. Besides, its motivation from particle physics, which was originally based on solving the so-called hierarchy problem (the puzzling discrepancy between the mass scales of electroweak interactions and gravity), has become stronger due to the apparent need for 100 GeV - 10 TeV scale supersymmetry to achieve unification of the gauge couplings in view of LEP results , and the prediction that the lightest Higgs boson should be below 135 GeV, as may be favoured by LEP data .

Recently, there has been discussion about the constraints on the MSSM which follow from the measurements of $`(g2)_\mu `$, the anomalous magnetic moment of the muon . The new set of data recently released shows a mild discrepancy with the Standard Model (but only at the 2 – 3 $`\sigma `$ level), which could in principle be due to supersymmetry .

The relic density calculation in the MSSM for a given set of parameters is nowadays accurate a few % or so, thus matching the precision given by cosmological measurements. A recent important improvement is the inclusion of coannihilations, which can change the relic abundance by a large factor in some instances . Much of the effort that has gone into this field has resulted in publicly available computer program packages, for instance DarkSUSY .

### 3.1 Detection methods

#### 3.1.1 Direct detection

If neutralinos are indeed the CDM needed on galaxy scales and larger, there should be a substantial flux of these particles in the Milky Way halo. Since the interaction strength is essentially given by the same weak couplings as, e.g., for neutrinos there is a non-negligible chance of detecting them in low-background counting experiments. Due to the large parameter space of MSSM, even with the simplifying assumptions above, there is a rather wide span of predictions for the event rate in detectors of various types.

The experimental situation is becoming interesting as several direct detection experiments after many years of continuing sophistication are starting to probe interesting parts of the parameter space of the MSSM, given reasonable, central values of the astrophysical and nuclear physics parameters.

The most enigmatic finding is that of DAMA , which after seven years of running in the low-background environment in Gran Sasso, now reports a visible seasonal variation  of the scattering rate which is formally at the $`6\sigma `$ level.

This seems to be difficult to reconcile with results from the CDMS  and EDELWEISS experiments , which have recently improved their limits by an order of magnitude. The data released essentially exclude all the parameter range for supersymmetric dark matter indicated by the DAMA experiment . However, it is difficult to compare the experiments directly, since they use different materials and very different techniques. For instance, it is not excluded that a large spin-dependent interaction may explain the DAMA results .

Although neutralino dark matter may, with careful choice of physical and astrophysical parameters , explain the DAMA results, the rates in minimal supersymmetry are generally lower. Even the CDMS and EDELWEISS experiments are so far only scratching the surface of the parameter space that is provided by the MSSM. However, a new generation of experiments is being prepared with an order of magnitude better sensitivity . The urgent challenge to the experimental dark matter community is now to build an independent, and preferably more sensitive detector which is similar enough to DAMA to enable a direct comparison .

### 3.2 Indirect detection

There is also the possibility of indirect detection through neutralino annihilation in the galactic halo or in astrophysical objects like the Sun and the Earth. This is becoming a promising method thanks to very powerful new detectors for cosmic gamma rays and neutrinos planned and under construction.

There has recently been balloon-borne detection experiments , where an excess of positrons over that expected from ordinary sources has been found. However, since there are many other possibilities to create positrons by astrophysical sources, e.g. near the centre of the Milky Way, the interpretation is not yet conclusive. In particular, it has been difficult to reproduce the apparent “bump” seen at 7 – 8 GeV in any model based on WIMP annihilation (for recent attempts, see ).

Antiprotons, $`\overline{p}`$, from neutralino annihilations were long hoped to give a useful signal , and there have been several balloon-borne experiments  performed and a very ambitious space experiment, AMS, to search for antimatter is under way . For kinematical reasons, antiprotons created by pair-production in cosmic ray collisions with interstellar gas and dust are born with relatively high energy, whereas antiprotons from neutralino annihilation populate also the sub-100 MeV energy band.

However, it was found  that the cosmic-ray induced antiprotons may populate also the low-energy region to a greater extent than previously thought, making the extraction of an eventual supersymmetric signal much more difficult. Recently, a thorough reanalysis  has arrived at the same conclusion.

Another problem that plagues estimates of the signal strength of both positrons and antiprotons is the uncertainty of the galactic propagation model and solar wind modulation.

Even allowing for large such systematic effects, the measured antiproton flux gives, however, rather stringent limits on MSSM models with the highest annihilation rates. One can also use the experimental upper limits to bound from below the lifetime of hypothetical $`R`$-parity violating decaying neutralinos . There may in some scenarios with a clumpy halo  (which enhances the annihilation rate) be a possibility to detect heavy neutralinos through spectral features above several GeV , something which will be probed by the upcoming PAMELA satellite project .

A very rare process in proton-proton collisions, antideuteron production, may be less rare in neutralino annihilation . However, the fluxes are so small that the possibility of detection seems marginal even in the AMS experiment.

With these problems of positrons and antiprotons, one would expect that problems of gamma rays and neutrinos are similar, if they arise from secondary decays in the annihilation process in the halo. For instance, the gamma ray spectrum arising from the fragmentation of fermion and gauge boson final states is quite featureless and gives the bulk of the gammas at low energy where the cosmic gamma ray background is severe. Also, the density of neutralinos in the halo is not large enough to give a measurable flux of secondary neutrinos, unless the dark matter halo is very clumpy . However, neutrinos can escape from the centre of the Sun or Earth, where neutralinos may have been gravitationally trapped and therefore their density enhanced.

Also, gamma-rays may be emitted from places where the dark matter density is substantially higher than average, e.g. in dark matter clumps , or near the galactic center . For a recent attempt to interpret the various anomalies in the measured fluxes of gamma rays, positrons and antiprotons in terms of neutralino dark matter, see . As a general comment, one should observe that the systematic uncertainties that come from unknown astrophysical parameters and propagation functions are much larger than the statistical errors used in .

### 3.3 Gamma ray lines

Gamma rays may also result from loop-induced annihilations  $`\chi \chi \gamma \gamma `$ or $`\chi \chi Z\gamma `$.

This would give monoenergetic photons with $`E_\gamma =m_\chi `$ or $`E_\gamma =m_\chi (1m_Z^2/4m_\chi ^2)`$ from the halo. The detection probability of a gamma line signal depends on the very poorly known density profile of the dark matter halo. N-body simulations have given a clue to the final halo profile obtained by hierarchical clustering in a CDM scenario . It turns out that the universal halo profile found in these simulations has a rather significant enhancement $`1/r`$ near the halo centre. If applicable to the Milky Way, this would lead to a much enhanced annihilation rate towards the galactic centre, and also to a very characteristic angular dependence of the line (and continuous gamma) signal. This would be very beneficial when discriminating against the galactic and extragalactic $`\gamma `$ ray background, and Air Cherenkov Telescopes (ACTs) would be eminently suited to look for these signals, if the energy resolution is at the $`1020`$ % level.

The models which give the highest rates should be within reach of the new generation of ACTs presently being constructed. These will have an effective area of almost $`10^5`$ m2, a threshold of some tens of GeV and an energy resolution approaching 10 %. For low-mass models, the space-borne telescope GLAST will have a better sensitivity. (See  for details.)

Another possibility to detect dark matter in gamma-rays has recently been investigated . If N-body simulations of structure formation are taken seriously, it appears that the average enhancement of the integrated signal from all cosmic structure in the universe would be several orders of magnitude compared to the case when the dark matter density only scales with the cosmic dilution factor $`(1+z)^3`$. The signature would be a continuum from neutralino annihilations plus a characteristic redshift-smeared line with a very rapid fall-off beyond the energy corresponding to the neutralino mass.

### 3.4 Indirect detection through neutrinos

More model-independent predictions, where essentially only the relatively well-determined local halo dark matter density is of importance, can be made for neutrinos from the centre of the Sun or Earth, where neutralinos may have been gravitationally trapped and therefore their density enhanced. (However, also the less well known velocity distribution also does affect the capture rate, weakening the predictions somewhat, especially for the Earth.) As neutralinos annihilate, many of the possible final states give after decays and perhaps hadronization energetic neutrinos which will propagate out from the interior of the Sun or Earth. In particular, the muon neutrinos are useful for indirect detection of neutralino annihilation processes, since muons have a quite long range in a suitable detector medium like ice or water. Therefore they can be detected through their Cherenkov radiation after having been produced at or near the detector, through the action of a charged current weak interaction $`\nu _\mu +A\mu +X`$.

Detection of neutralino annihilation into neutrinos is one of the most promising indirect detection methods, and will be subject to extensive experimental investigations in view of the new neutrino telescopes (AMANDA/IceCUBE, Baikal, NESTOR, ANTARES) planned or under construction . The advantage shared with gamma rays is that neutrinos keep their original direction.

The detector threshold energy, for “small” neutrino telescopes like Baksan, MACRO (now discontinued) and Super-Kamiokande is around 1 GeV. Large area neutrino telescopes in the ocean or in Antarctic ice typically have thresholds of the order of tens of GeV, which makes them sensitive mainly to heavy neutralinos (above 100 GeV) .

The detector response is weighted towards high neutrino energies, both because the cross section $`\sigma _\nu `$ rises approximately linearly with energy and because the average muon energy, and therefore the range, also grow approximately linearly with $`E_\nu `$. Therefore, final states which give a hard neutrino spectrum (such as heavy quarks, $`\tau `$ leptons and $`W`$ or $`Z`$ bosons) are usually more important than the soft spectrum arising from light quarks and gluons.

The capture rate, and therefore the signal rate, of the Earth is very similar to the scattering rate in the materials used for spin-independent detection. It seems therefore that the new generation of direct detection experiments generally have a greater sensitivity than the new neutrino telescopes searching for signals from the Earth . For the Sun, however, neutrinos provide a complementary and in some cases superior method .

## 4  NON-SUSY CANDIDATES

The phenomenology of supersymmetric dark matter (neutralinos) may be very similar for other types of weakly interacting massive particles (WIMPs). However, one can also imagine models where the WIMP only couples to leptons. These leptonic WIMPs, or LIMPs, may at first seem essentially undetectable in present-day experiments. It may be shown, however, that in most cases they necessarily give energetic gamma rays in their annihilations, due to higher-order processes .

Of course, the current interest among theoretical physicists in exotic, but at least partly realistic models like branes and large extra dimensions has also led to a number of dark matter models. Just as an example, models where the Standard Model fields live in compactified, “universal” extra dimension has a rich dark matter phenomenology, where an orbifold parity plays the role of R-parity in supersymmetric models. Detection rates turn out to be similar to the SUSY case but some interesting differences in the detailed phenomenology .

It finally should be noted with caution, that there are many models for dark matter where the candidate only interacts very weakly, e.g. gravitationally, with other matter. For some models, there still may exist cosmological signals from early decays , but it may not be excluded that detection by other methods is impossible. Still one may hope to learn more about the detailed distribution of dark matter, for instance by studying their effects on gravitational lensing. There are in fact studies  that interpret flux ratios of multiple gravitational images in terms of clumpiness of galactic halos, in rough agreement with what is found in N-body simulations.

## 5  CONCLUSIONS

To conclude, non-baryonic cold dark matter is required to explain new cosmological data, in particular on the microwave background. The fact that the favoured value of $`\mathrm{\Omega }_M`$ has gone down during the last decade from near 1 to around 0.3 is good news for detection, since larger cross sections generally means lower relic density. In particular this is true for the main particle physics candidate, the neutralino, which we have presented in some detail here. Direct detection experiments are rapidly developing their discovery potential. Indirect detection methods may be very useful complements. In particular, new gamma-ray and neutrino telescopes may have the sensitivity to rule out or confirm the supersymmetry solution of the dark matter problem. And of course, when the LHC starts taking data, we may get direct clues to the nature of the particle that most likely makes up the dark matter.

## 6  ACKNOWLEDGEMENTS

I wish to thank my collaborators, in particular Ted Baltz, Joakim Edsjö, Paolo Gondolo, Mia Schelke and Piero Ullio, for many helpful discussions. This work has been supported in part by the Swedish Research Council (VR).


format: text

text:
# Epitaxially strained [001]-(PbTiO3)1(PbZrO3)1 superlattice and PbTiO3 from first principles

## I INTRODUCTION

Perovskite oxide superlattices are expected to have properties distinct from those of the bulk constituents. Due to lattice mismatch, constituent layers below the critical thickness generally are in highly strained states. As can be confirmed by first-principles studies of bulk materials, the lattice instabilities and properties of bulk perovskites are extremely sensitive to the strain state. The degree of strain can be further tuned by varying the lattice constant of the substrate on which the superlattice is grown.

This discussion is of practical interest as it has been shown that perovskite oxide superlattices can indeed be grown with atomic scale precision Schlom ; Eckstein ; Christen , opening up a highly varied family of artificially structured materials for investigation. Thus, first-principles prediction of superlattice properties can directly lead to productive interactions with experimentalists and progress in materials designBT ; BTST .

In this work, we consider the simplest possible (and as yet hypothetical) -(PbTiO3)m(PbZrO3)n superlattice with $`m`$=1 and $`n`$=1.

As shown in Fig. 1, this consists of single unit cell layers of PbZrO3 (PZ) and PbTiO3 (PT) alternating along . The equilibrium cubic lattice constant of bulk PZ is significantly greater than that of PT, with the computed values being $`a_{0,PZ}`$=7.77 a.u. and $`a_{0,PT}`$=7.37 a.u., leading to a 5% mismatch (the corresponding experimental values are 7.81 a.u., 7.50 a.u, and 4%). At their equilibrium lattice constants, both cubic PT and cubic PZ have an unstable three-fold degenerate polar phonon mode at the zone center. The lattice mismatch means that in the superlattice, PZ is under compressive in-plane stress, while PT is under tensile in-plane stress. If the structure is allowed to relax without breaking the central mirror plane symmetry (i.e. in space group P4/mmm), the in-plane lattice constant takes an intermediate value of 7.57 a.u.. The PZ unit cell elongates, giving a local tetragonal unit cell with $`c/a`$ $`>`$ 1, and conversely, the PT unit cell will shorten, resulting in a local tetragonal unit cell with $`c/a`$ $`<`$ 1. This lowering of unit cell symmetry leads to a splitting of the three-fold degenerate unstable mode in each layer, with the lower frequency mode(s) along the long direction(s). Thus, on qualitative grounds, one expects at least two unstable $`\mathrm{\Gamma }`$ modes of different symmetry, one associated with polarization along the normal in the PZ layer, and a second, two-fold degenerate, mode producing polarization in the plane of the PT layer. DFPT calculations of the phonon dispersion of the relaxed centrosymmetric superlattice show two unstable zone-center polar modes, one with in-plane polarization confined in the PT layer and the other with polarization along the normal confined in the PZ layer, quantitatively confirming this picture Bungaro . The energy of the system can be lowered by either of these modes alone; the actual ground state might involve choosing one, or the other, or coupling both together.

First principles ground-state-structure calculations have previously been performed for this system, regarding the 1:1  ordered PZT 50/50 supercell as an approximate realization of the PZT solid solution KrakauerWu . Under the constraint of experimental volume, $`c/a`$ = 1.035 and $`\alpha =\beta =`$ 90 degrees, it was found that the relaxed structure has monoclinic $`Cm`$ symmetry, with internal coordinates in reasonable agreement with those experimentally determined for the monoclinic PZT phase Noheda . Connecting this result to our discussion in terms of lattice instabilities, we see that the ground state is in fact obtained by freezing in the two modes together, giving a total polarization $`𝐏=(P,P,P_z)`$ rotated away from the c-axis. Piezoelectric coefficients computed for this structure in Ref. KrakauerWu,  have very large values. This is attributed to the low energy for polarization rotation; this mechanism was first suggested to explain the high piezoelectric response of single crystal PMN-PT FuCohen ; Shrout  and accounts for the enhancement not only of the piezoelectric response in the monoclinic phase of the superlattice, but also, as we will show below, for the enhancement of the dielectric response.

Under an epitaxial constraint imposed by a cubic  substrate, the in-plane lattice vectors are equal in length and make an angle of 90 degrees, while the system is otherwise free to relax. Our qualitative picture of the superlattice instabilities suggests that the two modes should respond differently to changes in the substrate lattice constant. As the substrate lattice constant increases, the mode with polarization along the normal should become higher in frequency, while the mode with in-plane polarization should soften. This change in relative stability is expected to affect the direction of the polarization in the ground state structure, in a manner consistent with previous calculations in which a different structural constraint, on $`c/a`$, was imposed.KrakauerWu

In this paper, we present a quantitative first-principles analysis of the effect of varying in-plane strain on the epitaxially constrained -(PbTiO3)1(PbZrO3)1 superlattice. In section II, we describe the details of the computations. In section III, we present results of the optimized structural parameters for in-plane lattice constants $`a_{}`$ ranging from well below $`a_{0,PT}`$, at which the PT unit cell is expected to be approximately cubic, to above $`a_{0,PZ}`$. This wide range of strain allows us more fully to investigate the effects of strain on structure and properties. As we will see, there is a sequence of phases, from the tetragonal c-phase at the smallest values of $`a_{}`$, through a monoclinic r-phase with increasing $`a_0`$ and finally to an orthorhombic aa-phase. The locations of the phase boundaries are determined by computing the Hessian matrix of each optimized structure, comprised of the zone-center force-constant matrix, the elastic constants and the quadratic-order coupling between zone-center atomic displacements and homogeneous strain, and extrapolating to zero eigenvalue. Near the phase boundaries and throughout the r-phase the dielectric and piezoelectric responses are high, the sensitivity to applied fields and stresses resulting from the ease with which the polarization can be changed through rotation. We relate these results to the lattice instabilities of the centrosymmetric reference structure obtained from DFPT computation of the $`𝐪=0`$ frequencies. The behavior of the superlattice is compared with that of pure PT, which is seen to exhibit a qualitatively different behavior. In section IV, we discuss this distinction using a two-mode vs one-mode picture, and suggest some ways that these results can be used in materials design.

## II Method

To predict the ground state structure and zone-center phonon frequencies, we use density functional theory (DFT) and density functional perturbation theory (DFPT)hks ; Dfptrev  within the local density approximation (LDA).CAPZ  Calculations in this work have been done using the PWscf package. PWSCF  We use Vanderbilt ultrasoft pseudopotentials ultrasoft  to describe the interaction between ionic cores and valence electrons, and plane-wave basis set with kinetic energy cut-off of 35 Ry. The augmentation charges, required by the use of ultrasoft pseudopotentials, are expanded with an energy cutoff of 350 Ry. Brillouin zone (BZ) integrations are performed with a 6$`\times `$6$`\times `$3 Monkhorst-Pack mesh.

The polarization P is computed using the linearized expression Pα= $`_{i\beta }`$Z$`{}_{}{}^{}{}_{i,\alpha \beta }{}^{}`$ u, where ui are the atomic displacements of the optimized structure relative to the centrosymmetric reference structure (described below), and Z are the atomic Born effective charge tensors. The values of polarization reported below are computed with the averaged effective charge tensors of pure PT and pure PZ given in Table II of Ref. Bungaro, , which have been shown to be a good approximation to those of the superlattice.

The epitaxial constraint from a cubic substrate is treated implicitly by constraining the in-plane lattice constant of the system. We consider a range of in-plane strains from -6% to +4%. As the reference in-plane lattice parameter to define the in-plane strain, we use the equilibrium in-plane lattice parameter computed for the optimized polar structure, $`a_0`$ = 7.73 a.u., for the superlattice and $`a_0`$ = 7.32 a.u. for PT. Thus, the system at zero in-plane strain has the minimum elastic energy.

For each value of the in-plane strain, we perform the structural optimization of the superlattice and of pure PT in two stages. First, the ideal double-perovskite 10-atom unit cell of the superlattice is relaxed without any symmetry breaking beyond that associated with the Zr/Ti ordering (tetragonal space group $`P4/mmm`$, $`\mathrm{\#}`$123). Similarly, the perovskite 5-atom unit cell of PT is relaxed in space group $`P4/mmm`$. These will be referred to as the centrosymmetric reference structures. The zone center phonon modes of these structures are computed using DFPT. This shows the presence or absence of unstable modes that provide a guide to energy-lowering distortions.

Then, the symmetry is broken to allow relaxation into a lower symmetry phase, which we label following the notation of Pertsev et al..Pertsev  For the whole range of in-plane strain, we consider relaxations into the c-phase (tetragonal space group $`P4mm`$, $`\mathrm{\#}`$99) with polarization $`P\widehat{z}`$, the aa-phase (orthorhombic space group $`Cmm2`$, $`\mathrm{\#}`$35), with polarization $`P(\widehat{x}+\widehat{y})`$, and the r-phase (monoclinic space group $`Cm`$,$`\mathrm{\#}`$8), with polarization $`P(\widehat{x}+\widehat{y})+P_z\widehat{z}`$. We also investigated the ac-phase (monoclinic space group Pm, $`\mathrm{\#}`$6), with polarization $`P\widehat{x}+P_z\widehat{z}`$, and the a-phase (orthorhombic space group $`Pmm2`$, $`\mathrm{\#}`$25), with polarization $`P\widehat{x}`$; these two phases are found to be energetically unfavorable and do not appear in the present phase diagram. In principle, an equilibrium c-phase can be obtained by optimizing structural parameters from a starting r-phase or ac-phase (and an equilibrium aa-phase from a starting r-phase) though in practice more accurate results can be obtained by studying each of the space groups separately. Compatible k-point sampling allows us directly to compare the energies of these structures. Although the angle $`\gamma `$ between the c axis and the plane can in principle vary, in the total energy calculations it is held fixed at 90 degrees to simplify the calculation. The main effect of this approximation is to reduce the range of stability of the r-phase as a function of in-plane strain; for the systems considered here, we describe below how we compute the effect of strain through analysis of the Hessian matrix.

The zone-center dynamical matrices are computed for the optimized structures using DFPT. The zone-center phonon modes, obtained by diagonalizing the dynamical matrices, are used according to the formalism of Ref. dielectric,  to compute the phonon contribution to the zero-temperature constant-strain static dielectric tensor and to study its dependence upon epitaxial strain. We also used these results for a systematic study of the stability of the various phases as a function of in-plane strain. The full Hessian matrix of the second derivatives of the energy with respect to lattice-periodic atomic displacements and homogeneous strain is constructed with DFPT for the atomic displacements and finite difference stress and force results for the derivatives involving strain. Tracking the lowest eigenvalue of the Hessian matrix, we can accurately identify the instability point of a given phase as the in-plane strain where the lowest eigenvalue becomes negative Garcia . This allows a more accurate determination of the position of the phase boundary than total energy calculations alone, which are limited by computational precision, and also allows us to include the effects of the coupling to shear strain $`\eta _4`$ and $`\eta _5`$ in the r-phase, equivalent to allowing $`\gamma `$ to relax from 90 degrees.

## III Results

This section is organized as follows. In A, we present the calculations of the equilibrium structural parameters and polarization as a function of in-plane strain. Calculations of zone-center phonon modes, elastic constants, and coupling between the atomic displacements and homogeneous strain are performed to carry out the stability analysis and precisely locate the phase boundaries. In B, we present the results of the zone-center phonons in the centrosymmetric reference structure as a function of in-plane strain and use them to interpret the phase diagram. In C, we present the in-plane-strain dependence of the static dielectric tensor. Results in each subsection are presented first for the superlattice and then, for comparison, for pure PT.

### III.1 Structural parameters and phase diagram

In Fig. 2(a) we show the structural phase transitions of the superlattice with varying in-plane strain. The symmetry is evident from the nonzero components of the polarization. There is a continuous (second-order) transition from the c-phase, stable at large compressive in-plane strain, to the r-phase. In the r-phase, the polarization has almost constant magnitude and rotates smoothly in the ($`1\overline{1}0`$) plane from $`\widehat{z}`$ to $`\widehat{x}+\widehat{y}`$, followed by a second continuous transition from the r-phase to the aa-phase, stable at large tensile in-plane strain. The in-plane-strain dependence of the energy in each of these three symmetries is shown in Fig. 3(a). The a-phase and ac-phase energies computed at $`a_{}`$ = 7.57 a.u. were found to be significantly higher than that of the r-phase and thus are not included in further discussion. Accurate positions of the second-order c-r and r-aa phase boundaries can be best obtained not from intersection points of the total energy curves, but by interpolating the dependence of the quadratic-order energy around the equilibrium phase at each in-plane strain to obtain the in-plane strain at which the instability occurs.

We compute the Hessian matrix (phonons plus tetragonal strain, the latter coupling being zero for the c-phase and the aa-phase) and plot the lowest eigenvalue in Fig. 3(c). The zero-crossing of the eigenvalues can be seen to correspond to the points where the r-phase becomes more stable in Fig. 3(a); in the case of the r-aa transition the critical in-plane strain is clearly much easier to determine accurately in Fig. 3(c). To test the effect of fixing $`\gamma `$ to 90 degrees, we also compute the full Hessian matrix, including coupling to shear strains $`\eta _4`$ and $`\eta _5`$, and plot the lowest eigenvalue in Fig. 3(c). It can be seen that the effect of the additional strain relaxation is to expand the region of stability of the r-phase, the effect being slight for the c-r boundary and more significant for the r-aa boundary.

The corresponding results for pure PbTiO3 present an informative contrast to those for the superlattice. In Fig. 2(b) we show the structural phase transitions of pure PT with varying in-plane strain. The system makes a weakly first-order direct transition from the c-phase at large compressive in-plane strain to the aa-phase at large tensile in-plane strain. The intermediate r-phase present in the superlattice has been eliminated. There are small discontinuities in the $`c/a`$ ratio (inset of Fig. 3(b)) and the magnitude of the polarization at the phase transition. The optimized energies for each symmetry as a function of in-plane strain are shown in Fig. 3(b). Optimization in the monoclinic space group corresponding to the r-phase produces either the c-phase or the aa-phase depending on in-plane strain; these energies lie very slightly higher than those obtained by optimizing in the tetragonal space groups, reflecting the convergence criterion for the relaxation. We can confirm the weakly discontinuous nature of this transition by computing the Hessian matrix (phonons plus tetragonal strain) and plotting the eigenvalues in Fig. 3(d). At the transition in-plane strain, each of the two phases is marginally stable. To test the effect of fixing $`\gamma `$ to 90 degrees, we also compute the full Hessian matrix, including coupling to shear strains $`\eta _4`$ and $`\eta _5`$, and plot the lowest eigenvalue in Fig. 3(d). It can be seen that the effect of the additional strain relaxation introduces a narrow region of stability of the r-phase, between in-plane strains of 0.86$`\%`$ and 1.48$`\%`$. This window of 0.62$`\%`$ is about three times smaller than the T=0 width of 2$`\%`$ obtained by Pertsev et al.,Pertsev  the latter being an extrapolation to low temperature of a Landau theory fit near the paraelectric-ferroelectric transition temperatures above 700 K.

### III.2 Lattice instabilities of the centrosymmetric reference structure

Next, we investigate the relation of the ground state structure to the zone-center phonon modes in the centrosymmetric reference structure. In the unstrained superlattice ($`a_{}`$=7.57 a.u.), there are two unstable, nearly degenerate, polar zone-center modes, an $`E_u`$ mode with displacements confined in the PT layer, and an $`A_{2u}`$ mode confined in the PZ layer.Bungaro ; IrrepSymbols  In Fig. 4 we see that these behave differently with increasing in-plane strain, the PZ-like mode being stabilized by increasing in-plane strain, while the PT-like mode becomes increasingly unstable. This is consistent with the qualitative picture discussed in the introduction, with the mode with polarization along the normal being raised in frequency by increasing in-plane lattice constant, while the in-plane mode softens. The ground state structure reflects these opposing trends in the strength of the two unstable modes. At the smallest in-plane lattice constant, it is energetically favorable to freeze in only the $`A_{2u}`$ mode, producing the c-phase. As the in-plane lattice constant increases, the in-plane $`E_u`$ mode softens to the point that it also contributes to the ground state, producing the r-phase. The frequencies of the two modes cross at an in-plane lattice constant roughly in the middle of the r-phase. With further increase of the in-plane lattice constant, the $`A_{2u}`$ mode stiffens to the point that the $`E_u`$ mode is the only mode contributing to the ground state, producing the aa-phase. The analysis of the Hessian matrix in the ground state structures shows that the freezing-in of the polar modes eliminates the nonpolar instabilities also seen in Fig. 4, most notably the $`B_{2u}`$ mode.

The behavior of pure PT is qualitatively similar. There are two unstable polar zone-center modes $`A_{2u}`$ and $`E_u`$, increasing and decreasing in frequency with increasing lattice constant, respectively. The crossing of the two modes is likewise very close to the c-aa phase boundary, where the r-phase appears when shear strain relaxation is allowed.

For comparison with the unstable modes in the superlattice, in Fig. 4(c) we show the in-plane strain dependence of the zone-center unstable modes of pure PZ. For the unstrained superlattice we have previously shown Bungaro  that the atomic displacements of the lowest Eu mode, mainly confined in the PT layer, are very close to the displacements of the polar unstable mode of PT bulk. Similarly the atomic displacements of the lowest A2u mode, mainly confined in the PZ layer, are very close to the displacements of the polar unstable mode of PZ bulk. By comparing the phonon frequencies in Fig. 4, we see that certain modes in the superlattice can be identified with modes in pure strained PT and PZ. The PT-like mode of the superlattice, the lowest Eu mode, has a strain-dependent frequency very similar to that of the lowest Eu mode of pure PT. The PZ-like mode of the superlattice, A2u mode, compares well with the A2u mode of pure PZ. The B2u mode in the superlattice is a PZ-like mode with O-displacements mainly confined in the PZ layer and its strain-dependent frequency closely corresponds to that of the B2u mode of pure PZ.

### III.3 Static dielectric tensor

The phonon contribution to the static dielectric tensor computed as a function of epitaxial strain is presented for the superlattice in Fig. 5(a) and (c).

In (a), we plot the component $`ϵ_{zz}`$, as for thin films this is relevant to the dielectric screening of an electric field applied across the film. It can be seen that this component diverges on both sides of the r-aa phase boundary, and is large in a substantial fraction of the r-phase. In (c), we plot the inverses of the three eigenvalues of the dielectric tensor; thus a value approaching zero corresponds to a divergence of $`ϵ`$. In the c and aa phases, all three eigenvectors are uniquely determined by symmetry. In the c-phase $`ϵ_3=ϵ_{zz}`$ is the “collinear” response, i.e. for a field applied along the direction of the spontaneous polarization, and $`ϵ_1=ϵ_2=ϵ_{xx}=ϵ_{yy}`$. In the aa-phase, $`ϵ_3`$ is the collinear component, along $`[110]`$, $`ϵ_2`$ is the in-plane component perpendicular to the polarization, and $`ϵ_1=ϵ_{zz}`$ is along the direction normal to the layers. Two eigenvalues of the dielectric tensor ($`ϵ_1=ϵ_2`$) diverge at the continuous c-r transition, and one ($`ϵ_1`$) at the continuous r-aa transition, through there the second noncollinear eigenvalue $`ϵ_2`$ is also fairly large.

In PT (Fig. 5(b) and (d)), the behavior of the dielectric constant in the c and aa phases is rather similar to that in the superlattice, with a near-divergence in $`ϵ_{zz}`$ as the boundary of the aa-phase is approached. However, comparison of Fig. 5(a) and (b) show that the enhancement of the r-phase in the superlattice gives a large $`ϵ_{zz}`$ over a much larger range of epitaxial strain.

## IV Discussion

The most striking features of the phase diagram of the epitaxially constrained -(PbTiO3)1(PbZrO3)1 superlattice are the width of the r-phase, the separation of the minimum energy in-plane lattice constants of the c and aa-phases, and the linearity of the energy with in-plane strain in the r-phase. These have significant implications for the coherence of lattice-mismatched thin films and for the dielectric response.

From Fig. 3(a), it can be seen that the elastic energy of the superlattice as a function of in-plane lattice constant is not parabolic. Rather, as the in-plane lattice constant decreases from 7.73 a.u., the energy rises linearly with a small slope to the minimum of the c-phase at 7.50 a.u., behaving parabolically only for smaller $`a_{}`$. This makes the film compliant to a broad range of substrates; the estimated critical thickness of fully coherent films is close to constant in this range. For tetragonal PbTiO3, in comparison, within each phase (c or aa), the energy function is parabolic, so the critical thickness behaves as usual, increasing inversely to the deviation of the in-plane lattice constant from the equilibrium value for that phase.

Based on the computed energies for the various phases of the superlattice in Fig. 3(a), one can predict the occurrence of structural phase transitions with increasing thickness, as the system changes from fully coherent to partially relaxed and finally to fully relaxed. For systems under compressive stress from the substrate, there should be a transition from the c phase to the r phase with increasing thickness, while for systems under tensile stress, a transition from the aa phase to the r phase is expected. For pure PT, systems under tensile stress are predicted to undergo a transition with increasing thickness from the aa-phase through a narrow r-phase into the c-phase.

The present analysis identifies the most stable superlattice structure with the ten-atom cell shown in Fig. 1. The possibility remains that there may be instabilities to structural distortions that double or further expand the unit cell, in the plane, along the stacking direction, or both. This may be relevant for making comparisons with experimental observations.

Certain components of the dielectric tensor diverge as the phase boundaries are approached. While the polarization induced by an electric field along the direction of the spontaneous polarization is rather modest, the noncollinear responses are of great interest. This polarization rotation mechanism, through the coupling of polarization and strain, has been invoked to explain the large piezoelectric response in single-crystal relaxors FuCohen ; Shrout , and to explain the large piezoelectric response computed for ordered PZT (with the same structure as the superlattice considered here) in Ref. KrakauerWu, . The r-phase, with its continuously variable direction for the spontaneous polarization, is expected to have particularly large responses. In particular, the $`zz`$ component of $`ϵ`$, the normal polarization induced by a field along the normal, will be enhanced by a noncollinear contribution throughout the r-phase.

Finally, we consider the origin of the different behavior of the two systems. First, we discuss the difference in equilibrium lattice constant of the c-phase and the aa-phase. In a simple perovskite, such as pure PT, the difference arises from coupling between the strain and the unstable polar mode. PT has an unusually strong strain coupling, leading to a $`c/a`$ of about 1.06 in the tetragonal ground state; the c-phase and aa-phase equilibrium in-plane lattice constants are 7.32 a.u. and 7.45 a.u., respectively. Analogous information about the behavior of pure PZ can be obtained from the published parameterization of first-principles resultsKingSmith  with the c-phase and aa-phase equilibrium in-plane lattice constants being 7.71 a.u. and 7.83 a.u., respectively. A simple elastic-energy model for the superlattice can be constructed as the sum of the elastic energies of the two layers. For each constituent, the elastic constants of the c and aa phases turn out to be similar, with the values for PZ about 70% larger than those for PT. In the superlattice, the first-principles equilibrium in-plane lattice constant of the c-phase is 0.07 a.u. less than that of the elastic-energy model, and in the aa-phase it is 0.05 a.u. greater, leading to a separation of the first-principles lattice constants of the two phases about twice as large as in the model. This deviation from a simple elastic energy model might be due to the thinness of the individual layers and/or to the importance of inter-layer interactions beyond lattice matching.

The stability of the r-phase is favored by in-plane and perpendicular unstable polar modes having both a similar strength and positive anharmonic coupling. In a simple perovskite, such as pure PT, tetragonal strain splits the three-fold degenerate mode to give a dominant unstable $`A_{2u}`$ mode for $`c/a`$ less than 1, and a dominant in-plane $`E_u`$ mode for $`c/a`$ greater than 1. If the ground state is obtained by freezing in only the dominant unstable mode, we expect a discontinuous change in character at the in-plane lattice constant corresponding to the lowest-energy cubic structure, as seen in PT. Even if the strengths of the two modes are similar, the fourth order anharmonic term coupling them must be large enough to favor this combination over the higher-symmetry structures (z, as in the c-phase, or along  as in the aa-phase). Its magnitude can vary greatly from material to material, being enhanced by coupling to particular types of strain, especially shear strain.

In the superlattice, it seems that separating the two components of the polarization in two spatially distinct regions helps to stabilize the r-phase structure. Even if the strengths of the unstable polar modes in the two layer are different, a low-energy structure can be obtained by combining the perpendicular mode with polarization along z in one layer (e.g. the PZ layer) with a mode of different character, with polarization along , in the other layer (e.g. the PT layer). This idealized picture will be somewhat modified by electrostatic considerations, which will tend to minimize the divergence of P by polarizing the PT layer along the z direction.

The usefulness of this simple picture lies in its applicability to longer period PT/PZ superlattices; such systems are also of interest as they are accessible to experimental investigation.

## V Conclusions

The phase diagrams for the -(PbTiO3)1(PbZrO3)1 superlattice and pure PT with varying in-plane strain on a cubic substrate have been obtained from first-principles density-functional total-energy calculations and a stability analysis involving the computation of zone center phonons and their coupling to strain. At the smallest $`a_{}`$, both systems are tetragonal with pseudo-perovskite $`c/a`$ greater than 1 and $`𝐏`$ along  (the c-phase), and at the highest $`a_{}`$, both are pseudo-tetragonal orthorhombic with pseudo-perovskite $`c/a`$ less than 1 and $`𝐏`$ along  (the aa-phase). However, at intermediate in-plane strain, the transitional behavior is quite different. In the superlattice, $`𝐏`$ rotates away from the normal to give a monoclinic structure (the r-phase), and at a higher value of $`a_{}`$, completes the rotation into the plane of the layers. In the monoclinic phase, the dielectric and piezoelectric response is high, the sensitivity to applied fields and stresses resulting from the ease with which the polarization can be changed through rotation. In contrast, in pure PT changing epitaxial strain results in a discontinuous transition from $`𝐏`$ along  to $`𝐏`$ in plane if only tetragonal strain is allowed, with full coupling to shear strain opening a narrow range of r-phase stability. While the characteristic divergence of certain components of the dielectric tensor is seen as the phase boundaries are approached, there is no significant contribution from polarization rotation. For both systems, we relate the observed ground states to the instabilities of a relaxed centrosymmetric structure computed via DFPT. From these results we conclude that in the superlattice, the presence of two unstable modes confined in different layers, of different symmetry and comparable strength significantly enhances the window for the monoclinic structure, which suggests new approaches for designing artificially structured materials with high dielectric and piezoelectric response.

###### Acknowledgements.

We thank D. Vanderbilt, J. Neaton, O. Dieguez, A. Antons, S. Tinte, M. H. Cohen and H. Krakauer for valuable discussions. This work was supported by ONR N00014-00-1-0261. The work of K.M.R. was performed in part at LMCP, Université de Paris VI and Materials Department, University of California at Santa Barbara. The majority of the computations were performed at the Department of Defense High Performance Computing Centers NAVO and ERDC.


text:
# Early Disk Evolution

## 1. Introduction

Disks are the natural byproduct of the formation of stars: the resting place of material with too high an angular momentum to participate fully in the initial gravitational collapse that forms the central protostar. Accretion of an appreciable fraction of the stars mass may (or may not) occur through the disk, and disks appear to be an essential component of the engine that powers protostellar jets. A somewhat bewildering array of physical processes come into play to determine disk structure and must be addressed if we wish to understand disk evolution, forming planetary systems, and especially to interpret essentially unresolved observations of disks.

To zeroth order the material in a disk is in orbit about the central object: in other words, the radial component of the gravitational field from the protostar is balanced primarily by rotation. Even moderately flattened structures supported by other means tend to be unstable to interchange instabilities - whether the support is by gas pressure, radiation pressure or magnetic stresses (e.g. Stehle & Spruit 2001), and must therefore be evolving rapidly. This implies that the internal signal speed in a disk that survives for many orbits is much less than the Keplerian speed. A slowly-evolving disk is thin because its scale height is determined by the balance between vertical pressure gradient and the tidal field of the central object.

One consequence is that the local dynamical time scale (i.e. the orbital time scale) is much shorter than the time scale on which the disk evolves. This makes physical sense, as a well-defined disk cannot otherwise be expected to exist for more than a few orbits. There is however, a flip side to this: the range of orbital time scales across the radial extent of the disk makes globally steady evolution most unlikely. It is unreasonable to expect that the radial variation of the processes that determine disk structure conspires to force the local behaviour to be just that needed to guarantee that $`\dot{M}`$ is constant all radii. Disk models that rely on this assumption to tie together the behaviour at different radii should be regarded with caution, though they *might* represent the disk in a time-averaged sense.

A variety of processes play a role in the evolution of protoplanetary disks. Some factors are well-defined and operate on a scale large enough to be amenable to observation, such as tidal fields and irradiation by UV from companion stars or nearby stars, the X-ray irradiation from the central star, and self gravity if the disk mass $`>`$30% of the stellar mass. Gravitational disturbance by planets forming within the disk is, of course, critical in the later stages of disk evolution.

The factors that operate on small-scales are not directly observable, requiring a resolution equivalent to the disk thickness or less. They do however control transport processes, the internal dynamics of the disk, and chemical evolution of the disk material. (This latter point may yield a way of teasing apart the processes that occur within disks.) Although much of the theoretical work on the physics operating on these scales is directed towards the identification of mechanisms responsible for angular momentum transport, note that protoplanetary disks need not necessarily be *accretion* disks. Although magnetospheric accretion and emission from boundary layers are observed in some systems, the total mass accreted onto the star may be quite small compared to the total disk mass, and so accretion may not be the defining characteristic of the disk structure.

Convection was an early proposed mechanism for the origin of angular momentum transport in protostellar disks (e.g. Lin & Papaloizou 1980). The idea was that convective cells couple different disk radii, thereby acting as an effective viscosity. The convection would be driven by internal heating due to the disipation associated wioth the energy lost by the accreting material. However, simulations show that angular momentum is transported inwards rather than outwards (e.g. Cabot et al. 1987a,b), and in any case the convective cells are thin in the radial direction (ironically because of angular momentum conservation) and would be relatively inefficient as a mechanism for driving accretion.

The leading contender is now the magnetic field, which can transport angular momentum through the generation of turbulence via the magnetorotational instability (Balbus & Hawley 1991,2000; Balbus 2003), or through the centrifugal acceleration of material from the disk surfaces and the formation of a jet (e.g. Blandford & Payne 1982; Königl & Pudritz 2000). At present, it is not known what the initial radial distribution of mass and magnetic flux is. However, the initial poloidal field will have a strength of at least several milligauss, even if the ambient molecular cloud field is not particularly dragged in during the collapse that forms the disk. Even this weak field is sufficient to act as a seed for the magnetorotational instability at 1 AU (Salmeron & Wardle, in preparation).

Protostellar disks are weakly ionized, so the role of the magnetic field is complicated by its dependence on the trace charged species in the gas. Thus the evolution of protostellar disks is linked to chemistry, which determines the thermal and ionisation structure of the disk, and also to the population of dust grains which evolves through the accumulation of ices, sticking or shattering as a result of grain collisions, advection by turbulence, grain drag effects, and settling to the midplane (e.g. Weidenschilling & Cuzzi 1993).

## 2. Magnetic Coupling

The ability of the magnetic field to drive turbulence via the MRI or to produce a disk-driven wind is determined by its coupling to the weakly-ionised fluid. This is achieved by the collisions of charged species drifting in response to the electric and magnetic fields. The coupling is sensitive to ionisation rate, density and the very different mobilities of electrons, ions and grains.

The magnetic field in a weakly-ionised medium evolves according to the induction equation

$$\frac{𝐁}{t}=\times (𝐯\times 𝐁)𝐜\times 𝐄^{}$$
(1)
where $`𝐄^{}`$ is the electric field in the rest frame of the fluid, which is related to the current density $`𝐉=c/4\pi \times 𝐁`$ by

$$𝐄^{}=\frac{𝐉}{\sigma _{}}+\frac{\sigma _\mathrm{H}}{\sigma _\mathrm{H}^2+\sigma _\mathrm{P}^2}\frac{𝐉\times 𝐁}{B}\left(\frac{\sigma _\mathrm{P}}{\sigma _\mathrm{H}^2+\sigma _\mathrm{P}^2}\frac{1}{\sigma _{}}\right)\frac{(J\times B)\times B}{B^2}.$$
(2)
Here $`\sigma _{}`$, $`\sigma _\mathrm{H}`$ and $`\sigma _\mathrm{P}`$ are the field-parallel, Hall and Pedersen conductivities and the three terms involving $`𝐉`$, $`𝐉\times 𝐁`$ and $`(𝐉\times 𝐁)\times 𝐁`$ correspond to Ohmic, Hall, and ambipolar diffusion respectively.

Generally there is a trend from ambipolar to Hall and then Ohmic diffusion with increasing density and decreasing ionisation rate, which is why calculations of core collapse and disk formation have generally assumed that ambipolar diffusion dominates whereas solar nebula studies traditionally employ a resistive approximation. However, the Hall regime spans a surprisingly wide range of intermediate conditions and cannot be neglected. Furthermore, Hall diffusion qualitatively changes the vector evolution of the magnetic field, exhibiting a dependence on the strength of the magnetic field and a sensitivity to handedness that breaks many simple planar geometries that could otherwise be maintained in the presence of ambipolar or Ohmic diffusion.

The broad range of ionisation rates and densities obtaining over the radial and vertical extent of protostellar disks means that each form of diffusion will dominate in different regions. The ionisation rate is dominated by x-rays from the central low-mass star in the upper few g cm-2 of the disk (Glassgold, Najita & Igea 1997) and then by cosmic rays to depths of a few hundred g cm-2 (Hayashi 1981). The increase in ionisation rate and decrease in density with increasing height produces a trend from Ohmic through Hall to ambipolar diffusion with increasing height. If sufficient electrons are present Ohmic diffusion may never become completely dominant even at the midplane. The biggest uncertainty in estimates of the degree of coupling of the magnetic field to the disk is the grain population, as ions and electrons stick to grain surfaces and the mobility of the typical charge carriers is much reduced. If, on the other hand, grains have settled to the disk midplane, Hall diffusion dominates at the midplane over disk radii from 0.1 to 100 AU (Sano & Stone 2002a).

By way of example, the resulting vertical profile of the conductivity tensor in a minimum solar nebular model at 1AU is plotted in the presence and absence of 0.1 $`\mu `$m grains in the top panel of Fig. 1. The much larger conductivities in the no-grain model within three scale heights of the midplane reflects the lack of grains – the charge is carried by more mobile free ions and electrons. The decoupling of ions from the magnetic field for densities $`>10^9`$ cm-3 means that Hall component is larger than the Pedersen component beween 1.5 and 4 scale heights, and is 80% of $`\sigma _\mathrm{P}`$ even at the midplane. In the 0.1 $`\mu `$m grain model, the Hall conductivity dominates the Pedersen conductivity within 4 scale heights of the midplane.

The variation of magnetic coupling with height above the midplane determines how deeply magnetic activity penetrates the disk and therefore the thickness of the “dead zone” (Gammie 1996) which is not affected by magnetic activity, and does not participate in accretion.

## 3. Magnetorotational Instability

One measure of whether the conductivity is sufficient for the magnetic field to interact with the disk material is whether the field is unstable to the magnetorotational instability. This is determined by comparing the coupling parameter $`\chi =\omega _c/\mathrm{\Omega }`$ to the ratio of Alfvén speed to sound speed, $`v_A/c_s`$. Here $`\omega _c`$ is the frequency above which ideal MHD breaks down and $`\mathrm{\Omega }`$ is the Keplerian frequency. If Hall diffusion is unimportant the criterion is $`\chi >v_A/c_s`$, if it is dominant (and the magnetic field has the correct orientation), then $`\chi >(v_A/c_s)^2`$ is necessary.

The coupling parameter is plotted as a function of height in the lower panel of Fig. 1. The entire thickness of the disk is magnetically active in the absence of grains, whereas in the single-size grain model the layers above 2.5 scale heights are active. In both cases, Hall diffusion is important throughout the active regions.

If grains have settled to the disk midplane, Hall diffusion dominates at the midplane over disk radii from 0.1 to 100 AU (Sano & Stone 2002a). The implications for the magnetorotational instability are twofold. First, in the presence of Hall diffusion the instability operates under much weaker coupling than in the resistive or ambipolar diffusion cases, and therefore the potentially magnetically turbulent region is much more extensive than previously thought (Wardle 1999): the magnetic “dead zone” in which the field does not couple to the disk material is more restricted. This is illustrated in Fig. 2, which compares the linear development of the instability in models with and without Hall diffusion (Salmeron & Wardle 2003). If Hall diffusion is neglected, the instability grows in the upper layers of the disk and there is a substantial dead zone ($`z/h<2`$). With Hall diffusion, the growth rate of the most unstable mode is increased somewhat and the mode penetrates deeper into the disk, to within a scale height of the midplane.

Second, in the nonlinear regime the turbulence grows and decays cyclically: as the magnetic field is amplified through the nonlinear growth of the instability, the charged particles become more tightly coupled to the magnetic field, Hall diffusion becomes less important and growth shuts off. The field then decays until Hall diffusion becomes important again, at which point the instability restarts and drives turbulence again (Sano & Stone 2002b).

## 4. Magnetically-driven Jets

So far, I have focussed on the dynamics of an initially weak field lying within the disk. However, coupling of a large-scale poloidal magnetic field to the rotation of the disk may produce a centrifugally-driven outflow. This mechanism has been suggested as the origin of protostellar jets, though it is unclear whether the outflow occurs just at the inner edge of the disk or extends over a significant fraction of the disk (Königl & Pudritz 2000). The poloidal field could plausibly be the remains of the partial dragging in of the original field threading the parent molecular cloud during the core collapse and disk formation.

The centrifugal acceleration mechanism relies on coupling the magnetic field to the disk and to the accelerated material, the field transferring angular momentum from the former to the latter. Clearly, as with the magnetorotational instability, there must be sufficient coupling for this mechanism to work. In addition, there must be some slippage so that the magnetic field is not dragged inwards by the accreting material. Thus there are constraints on the magnitude of the diffusivity of the magnetic field. What appears to be unappreciated is that the *nature* of the coupling (i.e. whether ambipolar, Hall or Ohmic diffusion is dominant) determines the magnetic field direction as the lines emerge from the disk surface. This in turn, controls the tendency of material to slide along the field lines by the rotation (e.g. if the magnetic field lines are swept back, the acceleration is less efficient).

Simulations are currently unable to solve the complete disk+protostellar jet problem because of the large dynamic range in density between the disk and wind. Therefore disk-wind simulations are generally restricted to the wind region, with the magnetic field, density and fluid velocity specified at the base of the jet (roughly equivalent to the disk surface). Physically, these quantities are not independent but are related through the diffusion of the field within the disk, which simultaneously controls the bending of the field lines and how matter is loaded onto them and lifted away from the disk (Wardle & Königl 1993; Wardle 1997). This is still very much an open problem.

## 5. Summary

The examples discussed here – magnetorotational instability and disk winds – show how the coupling between a magnetic field and the disk material likely plays a crucial role in the structure and evolution of protostellar disks. It should be borne in mind that these two examples need not be mutually-exclusive: within the disk the large-scale poloidal field may be the mean of a field which is tangled (cf. Fig. 3 of Blandford & Payne 1982) by turbulence driven by the magnetorotational instability.

The magnetic activity in protostellar disks is stratified because of the strong dependence of diffusion on the level of ionisation. This implies that when grains are present, the bulk of the disk is magnetically inactive within a few AU of the centre, where it is self-shielded from ionisation sources. Much of the magnetic activity, and hence accretion, occurs in the surface layers of the disk (Gammie 1986). When grains have settled to the midplane, Hall diffusion of the magnetic field dominates over much of the disk (Balbus & Terquem 2000, Sano & Stone 2002b). This tends to increase the extent of the magnetically-coupled region within the disk, and reduces the extent of the dead zone.

The separation of the disk into turbulent and inactive layers affects planet-building, which is seeded by the coagulation and settling of solid material to the midplane (e.g. Weidenschilling & Cuzzi 1993). This is complicated by the intriguing coupling between the grain population and the magnetic field. The grains determine the field diffusivity and therefore affects the fluid dynamics, the fluid motions in turn feed back on the grain population through the competing processes of coagulation, shattering, mantle accumulation or evaporation and settling mediated through the stirring and energy dissipation by magnetically-driven turbulence in the disk.

The interdependence of the grains and fluid dynamics is, at present, unexplored territory. There is some hope for constraining theoretical models using observational signatures of disk chemistry, which also be dependent on the stratification of heating and the evolution of the grain population.

## References

Balbus, S. A. 2003, ARAA, 41, 555

Balbus, S. A. & Hawley, J. F. 1991, ApJ, 376, 214

Balbus, S. A. & Hawley, J. F. 2000, Sp Sci Rev, 92, 39

Balbus, S. A. & Terquem, C. 2001, ApJ, 552, 235

Blandford, R. D. & Payne, D. G. 1982, MNRAS, 199, 883

Cabot, W., Canuto, V. M., Hubickyj, O. & Pollack, J. B. 1987a, Icarus, 69, 387

Cabot, W., Canuto, V. M., Hubickyj, O. & Pollack, J. B. 1987b, Icarus, 69, 423

Gammie, C. F. 1996, ApJ, 457, 355

Glassgold, A. E., Najita, J. & Igea, J. 1997, ApJ, 480, 344

Hayashi, C. 1981, Prog Theor Phys Supp, 70, 35

Konigl, A. & Pudritz, R. E. 2000, in Protostars and Planets IV, ed. V. Mannings, A.P. Boss, & S. S. Russell (Tucson: U. Arizona), 759

Lin, D. N. C. & Papaloizou, J. 1980, MNRAS, 191, 37

Salmeron, R. & Wardle, M. 2003, MNRAS, 345, 992

Sano, T. & Stone, J. M. 2002a, ApJ, 570, 314

Sano, T. & Stone, J. M. 2002b, ApJ, 577, 534

Stehle, R. & Spruit, H. C. 2001, MNRAS, 323, 587

Wardle, M. 1997, in Proc. IAU Colloq. 163, Accretion Phenomena and Related Outflows, ed. D. Wickramasinghe, L. Ferrario & G. Bicknell (San Francisco: ASP), 561

Wardle, M. 1999, MNRAS, 307, 849

Wardle, M. & Königl, A. 1993, ApJ, 410, 218

Weidenschilling, S. J. & Cuzzi, J. N. 1993, in Protostars and Planets III, ed. E. H. Levy & J. I. Lunine (Tucson: University of Arizona Press), 1031


text:
# Some Recent Peculiarities of the Early Afterglow

## 0.1 Introduction

The early afterglow and the GRB/afterglow transition are among the unexplored regimes of GRBs. A great progress on this front was achieved during the last year when, following fast HETE II identifications, two afterglows (GRBs 021004, 030329) were followed from very early on showing remarkable variability and rich structure. Also, somewhat unexpectedly, a search within the BATSE/EGRET archives revealed a new very hard long lasting ($`200`$s) component of GRB 941017. This is most likely a manifestation of the early afterglow and of the GRB/afterglow transition. We discuss these developments and their implication within the internal-external shocks model.

## 0.2 The Prompt High Energy Emission from GRB 941017 as a GRB/afterglow transition

Recently, González et al. Gonzalez et al. (2003) discovered a high energy tail that extended up to $`200`$MeV in the combined BATSE and EGRET data of GRB 941017. The tail had a hard spectral slope ($`F_\nu \nu ^0`$) up to $`200`$MeV. It appeared $`10`$-$`20`$s after the beginning of the burst and displayed a roughly constant flux while the lower energy component decayed. At late times ($`150`$s after the trigger) the very high energy ($`10`$-$`200`$MeV) tail had a luminosity $`50`$ times higher than the “main” $`\gamma `$-ray energy band ($`30`$keV-$`2`$MeV).

Granot & Guetta Granot & Guetta (2003) where the first to suggest that we see here another manifestation of the very early afterglow. Sari Sari (1997) has shown that for long bursts the external shocks (and hence the afterglow) begin $`R/c\mathrm{\Gamma }^2`$ after the beginning of the burst while the internal shocks are still going on and hence the burst is still active. This behavior was seen in the transition from the harder initial burst to the softer early afterglow Burenin et al. (1999); Giblin et al. (1999); Piro00 . It was also seen in the prompt optical flash of GRB 990123 where the lower (optical) energy component did not trace the $``$MeV $`\gamma `$-rays and a pronounced hard to soft evolution was seen in the $`\gamma `$-ray signal.

The very high energy of this emission suggests that it is inverse Compton. There are two relevant emitting regions in the early afterglow: the forward shock and the reverse shock. With typical parameters Sari & Piran (2003) the energy of the synchrotron photons and the electrons’ Lorentz factor in the forward and reverse shocks are:

$$\begin{array}{ccc}\nu _{\mathrm{synch},\mathrm{F}}0.1\mathrm{MeV}(\mathrm{\Gamma }/300)^4;& \gamma _{e,\mathrm{F}}10^4(\mathrm{\Gamma }/300)& \\ \nu _{\mathrm{synch},\mathrm{R}}1\mathrm{eV}(\mathrm{\Gamma }/300)^2;& \gamma _{e,\mathrm{R}}300& \end{array}$$

There are four possible combinations of seed photons and scattering electrons:

$$\begin{array}{ccc}& \mathrm{Rev}.\mathrm{shock}\mathrm{electrons}& \mathrm{For}.\mathrm{shock}\mathrm{electrons}\\ \mathrm{Rev}.\mathrm{shock}\mathrm{photons}& 0.1\mathrm{MeV}(\mathrm{\Gamma }/300)^2& 100\mathrm{MeV}(\mathrm{\Gamma }/300)^3\\ \mathrm{For}.\mathrm{shock}\mathrm{photons}& 10\mathrm{GeV}(\mathrm{\Gamma }/300)^4& 10\mathrm{TeV}(\mathrm{\Gamma }/300)^5(\mathrm{within}\mathrm{Klein}\mathrm{Nishina})\end{array}$$

While these approximate results are very sensitive to $`\mathrm{\Gamma }`$, they indicate that inverse Compton scattering of the reverse shock photons on the forward shock electrons yeilds the right energy range. The detailed calculations of Pe’er & Waxman Pe’er & Waxman (2003) confirm these naive estimates. The main problem, however, is not to explain the location of the spectral peak but to explain the spectral slope ($`F_\nu \nu ^0`$) and the temporal slope ($`F_\nu t^0`$). Pe’er & Waxman Pe’er & Waxman (2003) reproduce the spectral slope by requiring that the synchrotron self absorption frequency of the reverse shock emission would be high enough to effect the observed spectrum. Granot & Guetta Granot & Guetta (2003) reproduce both the temporal and spectral behaviors by considering a slightly different scenario where the high energy component is produced by Synchrtron self Compton within the reverse shock. They require a slightly higher external density with a somewhat unusual profile, $`R^1`$, for a uniform ejecta shell (which may explain the rareness of the event). In both models the high energy component is a clear manifestation of the onset of the afterglow and the GRB/afterglow transition.

## 0.3 Afterglow Light curve variability

Theory: The different scenarios that lead to afterglow temporal variability can be distinguished according to their characteristic features. Density variations produce only weak fluctuations above the cooling frequency, $`\nu _c`$, and cannot produce sharp changes in the light curve. Energy variations produce variability both above and below $`\nu _c`$, and can arise either due to refreshed shocks or due to a patchy shell structure. These two mechanisms produce very different light curves. While the former produce a step-wise increase in the light curve, the later produces random fluctuations with a decreasing amplitude.

Rees & Mészáros Rees & Meszaros (1998), Kumar & Piran Kumar & Piran (2000a) and Sari & Mészáros Sari, R. & Meszaros (2000) suggested that slow shells take over the slowing down matter behind the afterglow shock and produce refreshed shocks. Slow shells with $`\mathrm{\Gamma }_s`$ emitted from the source right after the fast ejecta, catch up and collide with the slowing down ejecta at an observer time $`t0.25(\mathrm{\Gamma }_s/10)^{8/3}(E_{\mathrm{iso},52}/n_0)^{1/3}`$days (where $`E_{\mathrm{iso},52}`$ is the isotropic kinetic energy in units of $`10^{52}`$ergs and $`n_0`$ is the external density in cm<sup>-3</sup>) when the ejecta’s Lorentz factor drops slightly below $`\mathrm{\Gamma }_s`$. The clearest feature of refreshed shocks is a monotonous increase in the overall energy. Therefore the observed flux can only increase (relative to the expected decay). The light curve has a step wise form with each step produced by the arrival of a single shell. This step wise structure is seen both above and below the cooling frequency with a similar amplitude. Each step (in the optical light curve) should be accompanied by a flare in low frequencies that is produced by the reverse shock which propagates back into the slow shell Kumar & Piran (2000a); Sari, R. & Meszaros (2000). The time scale, $`\mathrm{\Delta }t`$ of the steps and the corresponding flares depends on their timing relative to the jet break. Before the jet break the refreshed shocks are “locally” spherically symmetric and therefore the angular time imposes $`\mathrm{\Delta }tt`$ Kumar & Piran (2000a). The intensity of the reverse shock flare in this regime is calculated in Kumar & Piran (2000a) , and the decay after the peak is $`t^2`$ Sari & Piran (2003). In the post-break case the cold slow shells may not expand sideways (if cold enough). Then they keep their original angular size, $`\theta _j`$, which is smaller than $`\mathrm{\Gamma }_s^1`$. In this case $`\mathrm{\Delta }tt_j<t`$, where $`t_j`$ is the jet break time Granot, Nakar & Piran (2003) and the transition is fast. A reverse shock flare is expected in this case as well. However, its frequency, intensity and temporal decay (which is expected to be steeper than in the spherical case) are much harder to calculate.

Kumar & Piran Kumar & Piran (2000b) suggested, in the Patchy shell model, that the shells have an intrinsic angular structure. As the blast wave decelerates the angular size of the observed region ($`1/\mathrm{\Gamma }`$) increases. The effective (average) energy of the observed region and hence the observed flux, relative to the expected decay, varies with time depending on the angular structure. The variability time scale is $`\mathrm{\Delta }tt`$ Nakar & Oren (2003). The averaging over a larger and larger random structure leads to a decay of the envelope as $`t^{3/8}`$ Nakar & Oren (2003); Nakar, Piran & Granot (2003). An important feature of this scenario is the break of the axial-symmetry and therefore the production of a linear polarization. The variation of the polarization, both in degree and in angle, are correlated with the light curve variations Nakar & Oren (2003); Granot & Königel (2003). The variability will be observed both above and below the cooling frequency $`\nu _c`$ with a similar amplitude.

Wang & Loeb Wang & Loeb (2000), Lazzati et al. Lazzati et al. (2002) and Nakar et al. Nakar, Piran & Granot (2003) considered External density variations. Such variations may result from ISM turbulence or from a variable pre-burst stellar wind. Wang & Loeb Wang & Loeb (2000) analyzed the light curve resulting from mild density fluctuations due to ISM turbulence. They show that these density fluctuations can produce short time scale ($`\mathrm{\Delta }t<0.1t`$) and low amplitude ($`10\%`$) fluctuations in the light curve. Lazzati et al.Lazzati et al. (2002), Nakar et al. Nakar, Piran & Granot (2003) and Nakar & Piran Nakar & Piran (2003) considered large amplitude spherical density fluctuations. A basic feature of the resulting light curve that distinguishes it from energy variations is that in the former the light curve is different above and below the cooling frequency, $`\nu _c`$. Density variations produce only weak fluctuations above $`\nu _c`$. The amplitude of the fluctuations above $`\nu _c`$ is at most tens of percents and it is much smaller than the amplitude of the fluctuations below $`\nu _c`$ Nakar & Piran (2003).

A second feature of density fluctuations is their inability to produce a sharp variation (either increase or decrease) in the light curve Nakar Piran & Granot (2003). First, we note that because of angular spreading, spherical density drops cannot produce decays sharper than $`t^{2.6}`$ Nakar & Piran (2003); Kumar & Panaitescu (2000) and even this decay is reached very slowly. More interesting is the fact that even a sharp density enhancement cannot produce a steep increase in the light curve. The earlier calculations Nakar, Piran & Granot (2003); Lazzati et al. (2002); Nakar & Piran (2003) assumes that the ejecta can be described by a Blandford-McKee solution whose density profile varies instantaneously according to the external density. These calculations do not account, however, for the reverse shock resulting from density enhancement and its effect on the blast-wave. Thus the above models are limited to slowly varying and low contrast density profiles. Now, the observed flux depends on the external density, $`n`$, roughly as $`n^{1/2}`$. Thus, a large contrast is needed to produce a significant re-brightening. Such a large contrast will, however, produce a strong reverse shock which will sharply decrease the Lorentz factor of the emitting matter behind the shock, $`\mathrm{\Gamma }_{sh}`$, causing a sharp drop in the emission below $`\nu _c`$ and a long delay in the arrival time of the emitted photons (the observer time is $`\mathrm{\Gamma }_{sh}^2`$). Both factors combine to suppresses the flux and to set a strong limit on the steepness of the re-brightening events caused by density variations. Note that while non spherical density fluctuations may lead to a steeper decline they usually do not lead to a steeper increase in the flux.

Implications: The early afterglow of GRB 021004 showed clear deviations from a smooth power law decay, lasting from 0.04 days to 3 days. The fluctuations in the light curve where accompanied by fluctuations both in the degree and in the angle of the polarization (Nakar & Oren (2003) and references therein).

The steep decays after each bump imply that the variations do not result from refreshed shocks. Thus, variable external density variations Nakar, Piran & Granot (2003); Lazzati et al. (2002); Heyl & Perna (2003) and the patchy shell model Nakar & Oren (2003); Nakar, Piran & Granot (2003) were considered as possible explanations. Unfortunately, the X-ray observations are not detailed enough to clearly distinguish between density and energy variations (although the former are favored by Heyl & Perna (2003)). However, the sharp decays cannot be produced by “locally” spherical density variations Nakar & Piran (2003). Furthermore, the first bump requires, using the instantaneous Blandford-McKee approximation, an increase in the external density by a factor of $`10`$ over $`\mathrm{\Delta }R/R0.05`$ Lazzati et al. (2002); Nakar & Piran (2003). Such a density contrast produces a mildly relativistic reverse shock which reduce $`\mathrm{\Gamma }_{sh}`$ by a factor of $`2`$, making the approximation inconsistent. Preliminary results Nakar Piran & Granot (2003) of detailed numerical simulations (including both hydrodynamics and synchrotron radiation) suggest that this bump cannot be produced by density enhancement (due to the suppressed forward shock emission, caused by the reverse shock). These results leave the patchy shell as the only viable explanation. Indeed, Nakar & Oren Nakar & Oren (2003) show that patchy shell can reproduce the light curve (including the sharp rise and steep decay of the first bump). They show that angular energy profiles which produce the observed light curve produce also a polarization curve that fits the observed polarization. We conclude that angular energy fluctuations are the dominant process that produce the observed fluctuations in GRB 021004.

In addition to the remarkable supernova signature, the optical afterglow GRB 030329 has shown also a unique variability (Granot, Nakar & Piran (2003) and references therein). Several step-wise bumps, at $`t=1.5,2.6,3.3`$ and $`5.3`$days, are seen after the jet break at $`t_j0.5`$days. The first bump was the largest and best monitored, but even with the less dense monitoring of the later bumps the step-wise profile (where after each bump the original decay slope is resumed) is clear. All the bumps had a short rise time $`\mathrm{\Delta }t0.4`$-$`0.8\mathrm{d}<t`$. The step-wise profile seems like a clear signature of post-break refreshed shocks Granot, Nakar & Piran (2003), where $`\mathrm{\Delta }tt_j<t`$ because the later slower shells did not expand sideways before colliding with the faster earlier ejecta. Moreover, the energy injected in these shocks is 10 times the energy in the original blast-wave. This late (or rather slow) energy injection explains an additional peculiarity of this GRB: the low energy output in $`\gamma `$-rays and in the early X-ray afterglow.

Berger et al. Bergeret al. (2003) suggested that the first bump and the energy deficiency can be explained by a two component jet: A slow and energetic component with a wide half-opening angle ($`17^{}`$) dominates the afterglow after $`1.5`$days, and a fast component with a narrow half-opening angle ($`5^{}`$) dominates the afterglow before $`1.5`$days. The slow component is observed only after $`1.5`$days$`t_{\mathrm{dec}}0.5(\mathrm{\Gamma }_s/10)^{8/3}(E_{\mathrm{iso},52}/n_0)^{1/3}`$days, since only at this time its reverse shock consumes the slow shell. However, this model, which received a great publicity, predicts that the rise time, $`\mathrm{\Delta }t`$, of the first bump should be of the same order as the observed time, $`t_{\mathrm{dec}}`$. Furthermore, it predicts a smooth light curve after the first bump. Both predictions are contradicted by the observations.

Millimeter observations of GRB 030329 Sheth et al. (2003) show that at 100 & 250 GHz the flux is rather constant during the first week. Most surprising are two measurements at $`100`$GHz, one before the first bump ($`0.6`$-$`1`$d) and one after ($`1.7`$-$`1.9`$d), which show a constant flux. Sheth et al. Sheth et al. (2003) and Berger et al. Bergeret al. (2003) claim that these results support the two component jet model and reject the refreshed shocks model due to the lack of radio flares. However this analysis overlooks the fact that the encounter of the slow component in the two-components jet model with the external matter produces a reverse shock which should produce a radio flash. This flash is analogous to the optical flash produced by the deceleration of the fast component Sari & Piran (2003), and to the radio flare expected in the refreshed shocks model. The timing of this flash is $`t_{\mathrm{dec}}=1.5`$d and its magnitude is easily calculated. The contribution from this reverse shock at $`t_{\mathrm{dec}}=1.5`$d is larger by a factor of up to $`\mathrm{\Gamma }_s^{5/3}15`$ than the flux from the forward shock. Thus a very bright and fast fading millimeter flash is expected in the two component jet model. A more detailed calculation (using Sari, R. & Meszaros (2000) and the parameters of the wide jet presented in Bergeret al. (2003)) shows that both $`\nu _a`$ and $`\nu _m`$ of the reverse shock at $`t_{\mathrm{dec}}`$ are around 100-200 GHz and that the flux at $`\nu _m`$ is $`500`$mJy which is an order of magnitude larger than the expected flux from the forward shock and than the observed fluxes at 100 & $`250`$GHz ($`50`$mJy) at this time.

The main argument of Bergeret al. (2003) in favor of a two-components jet is the existence of a second jet break in the radio after $`t_{j,2}10`$d. However, even this argument is not strongly supported by the data. According to this model the flux below $`\nu _m`$ should rise as $`t^{1/2}`$ before $`t_{j,2}`$ and decay as $`t^{1/3}`$ after $`t_{j,2}`$. However, the data of Sheth et al. (2003) contradict this prediction. At 100 & $`250`$GHz the flux is constant before the passage of $`\nu _m`$ and the decay after this passage (at $`t=6`$ & $`8\mathrm{d}<t_{j,2}`$, respectively) is steeper than $`t^{1.7}`$ in both bands. This looks like a clear signature of a post break behavior in the radio at $`t>1`$d. Now the radio observations at lower frequencies ($`22`$GHz) do not conform with the simple post-break model ($`F_\nu t^{1/3}`$). The flux at these frequencies rises with time before the passage of $`\nu _m`$. Interestingly enough, this radio behavior is exactly the one predicted as the post break radio behavior by Granot et al. (2001), using a 2D relativistic hydrodynamical simulations. We find the striking similarity between the radio observations and Fig. 2 of Granot et al. (2001) as a very strong support that the broad band data are totally consistent with a single jet.

Sheth et al., Sheth et al. (2003) emphasize the fact that radio flare Kumar & Piran (2000a) was not detected during the first bump at $`t1.5`$ days and argue that this rules out the refreshed shocks model. However, in the post-break refreshed shocks model, the radio flash is expected to be fainter and to decay faster than the radio flash in the two-components jet model (due to the lower energy and the lateral spreading of the slow shell as opposed to the more energetic and “locally” spherical wide jet). Thus, by assuming spherical symmetry, Sheth et al., Sheth et al. (2003) over estimate the expected flash in the post-break refreshed shocks scenario. It may be that the flash was missed by the sparse measurements due to its lower intensity and faster decay. A detailed (and highly non trivial) calculations should be done in order find out.

## 0.4 Summary

With an increasing flow of new observations we discover that GRBs and afterglows are richer than what was previously thought. The simple spherical theory had to be modified, first with jets and now with additional angular structure (patchy shells) and more extended velocity structure (refreshed shocks). The simple synchrotron theory has to be modified with inverse Compton scattering. One can worry, are we adding epicycles trying to revive a wrong theory? We don’t believe so. Complications and variation are common in astrophysics and are found everywhere in nature. Moreover, the three main themes that have been introduced here: patchy shells, refreshed shocks and inverse Compton, were not invoked aposteriori to explain the new observations. On the contrary, all three have been suggested long ago. It is just natural and even reassuring to discover them when better data become available.

The research was supported by US-Israel BSF.