text:
# Li intercalation effects on magnetism in undoped and Co-doped anatase TiO2

## Abstract

The effects of $`n`$-type carrier doping by Li intercalation on magnetism in undoped and Co-doped anatase TiO2 are investigated. We have found that doped $`n`$-type carriers in TiO2 are localized mainly at Ti sites near the intercalated Li. With increasing the intercalation, local spins are realized at Ti. In the case of Co-doped TiO2, most of the added $`n`$-type carriers fill the Co 3$`d`$ bands and the rest are localized at Ti. Therefore, Co magnetic moment vanishes by Li intercalation to have a nonmagnetic ground state.

Anatase TiO2 is a wide band gap (3.2 - 3.7 eV) semiconductor. This wide band gap property provides a wide range of applications, such as air and water purifications using photocatalysis which converts solar energy into electrochemical energy and batteries, electrochromic devices based on lithium intercalation. A theoretical study shows that Li intercalation is easier in anatase TiO2 than in rutile TiO2. This is due to the open structure of anatase TiO2. Recently, magnetic features are observed in the Li intercalated anatase TiO2. It is possible to intercalate Li atoms up to Li/Ti ratio $``$ 0.7. With increasing the Li intercalation, a structural transition occurs from tetragonal to orthorhombic one. Also the insulator-to-metal transition is observed for $`x`$$`>`$ 0.3. The localized moments were measured, 0.003 $``$ 0.004$`\mu _B`$ per Ti, through the Li intercalation.

On the other hand, room temperature ferromagnetism is observed in Co-doped anatase TiO2 thin film made by the combinatorial pulsed-laser-deposition molecular-beam epitaxy method. A sizable amount of Co, up to 8$`\%`$, is soluble in anatase TiO2. The measured saturated magnetic moment per Co ion was 0.32$`\mu _B`$ apparently in the low spin state and $`T_C`$ was estimated to be higher than 400 K. It is recognized that the carriers play role of inducing the ferromagnetism in dilute magnetic semiconductors.

To explore the carrier doping effects on magnetism, we have investigated electronic and magnetic properties of Li intercalated anatase TiO2 and Ti1-xCoxO2 ($`x`$=0.0625). Li/Ti intercalation ratios of 0.0625, 0.125, and 0.25 are considered for the undoped case, while 0.067 and 0.133 for the Co-doped case. We consider only the tetragonal anatase structure. We have used the linearized muffin-tin orbital (LMTO) band method in the local-spin-density approximation (LSDA). We have considered a supercell containing 16 f.u. in the primitive unit cell by replacing one Ti by Co or intercalating several Li ions ($`a`$=$`b`$=7.570, $`c`$=9.514 $`\mathrm{\AA }`$). Sixteen empty spheres are employed at the interstitial sites to enhance the packing ratio in the LMTO band calculation.

We have first calculated the electronic structure of Li intercalated TiO2. In all Li/Ti ratio cases, we have obtained metallic ground states (Fig. 1). For Li/Ti=0.0625, the paramagnetic ground state is obtained. However, in other cases, total energies of the paramagnetic and the ferromagnetic ground state are almost the same. Doped $`n`$-type carriers fill the Ti 3$`d`$ conduction band. Maximum localized magnetic moments at Ti are 0.029, 0.027 $`\mu _B`$ for Li/Ti=0.125 and 0.25, respectively. Thus the magnetic moment is not necessarily proportional to the number of localized electrons at Ti. The exchange splitting is clearly seen in the valence band top mainly with O 2$`p`$ characters (Fig. 1).

Now, we have performed band calculations for Co-doped case. For Li/Ti=0.067, we have obtained paramagnetic and insulating ground state in Fig. 2. Insulating ground state results from filling up Co 3$`d`$ band by the $`n`$-type carriers. Therefore, Co $`t_{2g}`$ band is fully occupied, and the total occupancy of $`d`$ states amounts to $`d^7`$($`t_{2g}^6e_g^1`$). As seen in Fig. 2, the position of occupied $`t_{2g}`$ band is different from the non-intercalated Co-doped case. In the latter case, the occupied Co $`t_{2g}`$ band is located near the valence band top, while, in the former case, the fully occupied Co $`t_{2g}`$ band is located near the conduction band bottom. The unoccupied Co $`e_g`$ band is hybridized with Ti 3$`d`$ conduction band.

For the Li/Ti=0.133 case, nearly paramagnetic and metallic ground state is obtained (Fig. 2). The carriers are mainly of Ti 3$`d`$ character. The extra electrons after filling up the Co $`t_{2g}`$ band occupy not Co $`e_g`$ band but Ti 3$`d`$ conduction band, because the carriers are localized at Ti sites near the intercalated Li. From these, one can expect that $`n`$-type carrier induced ferromagnet can be fabricated in TiO2 by simultaneously intercalating Li and doping some 3$`d`$ transition metals having high spin magnetic ground state such as Mn, Fe.

In conclusion, we have found that, by intercalating Li in TiO2, Ti atoms have localized magnetic moments, 0.029 and 0.027$`\mu _B`$ for Li/Ti=0.125 and 0.25, respectively. In the case of Co-doped TiO2, nonmagnetic ground states are obtained for Li/Ti=0.067 and 0.133.

Acknowledgments$``$ This work was supported by the KOSEF through the eSSC at POSTECH and in part by the BK21 Project.

Fig.1

Fig.2


format: text

text:
# 1 Scatter plot of Br(τ→μγ) against the lighter stau mass for the ansatz H₁. We take the SU(2) gaugino mass to be 200 GeV, A₀=0, μ>0, and tanβ=10. We assume both the normal and inverted hierarchies for the light-neutrino mass spectrum.

Lepton-flavor Violation (LFV) in the neutrino sector is observed by the oscillation experiments. The $`\nu _\mu `$-$`\nu _\tau `$ mixing was first discovered in the atmospheric neutrino , and is being confirmed furthermore by the K2K experiments . Also, the neutrino oscillation in the solar neutrino is confirmed by the SNO experiment , and it is considered to comes from $`\nu _e`$-$`\nu _\mu `$ mixing. The SNO and the superKamiokande experiments pin down the solar neutrino problem to the LMA solution almost. This will be confirmed soon by the Kamland experiments. Now, the neutrino oscillation experiments use human-made neutrinos, and we are going to next stage for LFV in the neutrino sector. Now we have a big question, that is, how large LFV in the charge lepton sector is. The neutrino oscillation comes from the small neutrino masses, and we need to know the origin of the neutrino mass in order to answer the question.

The seesaw model , in which the right-handed neutrinos are introduced, predicts the finite, but tiny neutrino masses in a simple way. In this model, the Yukawa coupling and Majorana mass terms are lepton-flavor violating. In addition to the small neutrino mass, this model is well-motivated from other theoretical points of view. First, in the SO(10) GUT, the unification of matter is realized in an economical and elegant way. Quarks and leptons in one generation are embedded in a 16 dimensional multiplet. Here, existence of the right-handed neutrinos is required. Second is the baryon asymmetry in Universe. Since the sphareron effect washed out the $`B+L`$ number in the early universe, the $`BL`$ generation is required. The Majorana masses of the the right-handed neutrino violate the lepton number, and thus, the leptogenesis  induced by the right-handed neutrinos is a good candidate of the baryogenesis.

It is considered that the seesaw model should be supersymmetrized since this model has a big hierarchy between the weak scale and the right-handed neutrino mass scale. The supersymmetry has a big role for LFV in the charged lepton sector. The radiative correction by the Yukawa coupling of the right-handed neutrino may generate the sizable LFV slepton masses , and the predicted $`Br(\mu e\gamma )`$ or $`Br(\tau \mu \gamma )`$ may be in the reach of the future experiments  .

In this paper we concentrate on the lepton-flavor violating tau decay in the supersymmetric seesaw model. Now the bounds on the lepton-flavor violating tau decay modes are being improved by the Belle experiment in the KEK B factory. The new bounds by the Belle experiment are following;

$`Br(\tau \mu \gamma )`$ $`<`$ $`5\times 10^7,`$ (1)
$`Br(\tau e^{}e^+e^{})`$ $`<`$ $`7.8\times 10^7,`$ (2)
$`Br(\tau \mu ^{}\mu ^+\mu ^{})`$ $`<`$ $`8.7\times 10^7,`$ (3)
$`Br(\tau e^{}\mu ^+\mu ^{})`$ $`<`$ $`7.7\times 10^7,`$ (4)
$`Br(\tau \mu ^{}e^+e^{})`$ $`<`$ $`3.4\times 10^7,`$ (5)
$`Br(\tau e^{}K^0)`$ $`<`$ $`5.8\times 10^7,`$ (6)
$`Br(\tau \mu ^{}K^0)`$ $`<`$ $`5.3\times 10^7,`$ (7)
for the processes with $`\mathrm{\Delta }L_\tau =1`$ and $`\mathrm{\Delta }L_{\mu /e}=1`$ . The Belle experiment is expected to obtain about 300$`fb^1`$ data two year after. If they can keep the present situation, the bound for $`Br(\tau \mu \gamma )`$ will be improved to be $`6\times 10^8`$. Also, if the super B factory starts, it can reach to $`2\times 10^8`$. Thus, now it is important to discuss impact of the LFV decay of tau lepton on the supersymmetric seesaw model.

First, we review the supersymmetric seesaw model briefly. The superpotential of the lepton sector is given as

$`W=N_i^c(Y_\nu )_{ij}L_jH_2+E_i^c(Y_e)_{ij}L_jH_1+{\displaystyle \frac{1}{2}}N_{}^{c}{}_{i}{}^{}(M_N)_{ij}N_j^c,`$ (8)
where the indexes $`i,j`$ run over three generations and $`(M_N)_{ij}`$ is the heavy singlet-neutrino mass matrix. While the light-neutrino mass matrix is given as

$`(_\nu )_{ij}`$ $`=`$ $`{\displaystyle \underset{k}{}}{\displaystyle \frac{(Y_\nu )_{ki}(Y_\nu )_{kj}}{M_{N_k}}}H_2^2`$ (9)
from this superpotential, the radiative LFV slepton masses induced by the Yukawa coupling of the right-handed neutrino are approximately proportional to an Hermitian matrix $`H`$,111 Here, we do not include the other LFV sources, such as in SUSY GUT . However, typical SUSY GUTs predict smaller $`Br(\tau \mu \gamma )`$ or $`Br(\tau e\gamma )`$ , which is suppressed by the KM mixing. Thus, our result for the LFV tau decay is expected not to be changed even if effect of the SUSY GUT is included.

$`H_{ij}`$ $`=`$ $`{\displaystyle \underset{k}{}}(Y_\nu ^{})_{ki}(Y_\nu )_{kj}\mathrm{log}{\displaystyle \frac{M_G}{M_{N_k}}},`$ (10)
where $`M_G`$ is the GUT scale.

Now we can parameterize the seesaw model by $`(_\nu )`$ and $`H`$. The seesaw model has 18 physical degrees of freedom in addition to 3 charged lepton masses, those are 6 eigenvalues, 6 mixing angles, 4 mixing CP phases, and 2 Majorana CP phases for the Yukawa coupling and the Majorana mass matrices. $`H`$ has 3 eigenvalues, 3 mixing angles, and 3 phases, while the light-neutrino mass matrix has 3 mass eigenvalues, 3 mixing angles, 1 mixing CP phase, and 2 Majorana CP phases. Thus, the degrees of freedom in the seesaw model match with those in $`(_\nu )`$ and $`H`$ . Unfortunately, the parameter space in $`(_\nu )`$ and $`H`$ is wider than that in the seesaw model, however, the parametrization works well to realize the seesaw model from the low-energy data including the neutrino-oscillation data . As will be shown explicitly later, the fact we can take $`(_\nu )`$ and $`H`$ almost independently means that the charged LFV search, such as $`\mu e\gamma `$ or $`\tau \mu \gamma `$, will give information independent of the neutrino oscillation for the seesaw model.

Now we show how large the branching ratios of the LFV tau decay modes can be, using our parameterization. Here, we adopt two ansatz for $`H`$. First is

$`H_1=\left(\begin{array}{ccc}a& 0& 0\\ 0& b& d\\ 0& d^{}& c\end{array}\right).`$ (14)
where $`a,b,c`$ are real and positive, and $`d`$ is a complex number. Here, $`H_{12}`$ and $`H_{12}H_{32}`$ are zero. If either of them is of the order of 1, $`Br(\mu e\gamma )`$ is much larger than the experimental bound. This ansatz may predict large $`Br(\tau \mu \gamma )`$ if $`H_{23}`$ is large. Second is

$`H_2=\left(\begin{array}{ccc}a& 0& d\\ 0& b& 0\\ d^{}& 0& c\end{array}\right).`$ (18)
This may predict large $`Br(\tau e\gamma )`$ when $`H_{13}`$ is large. From a viewpoint of model-building, $`H_1`$ may be more favored than $`H_2`$. If the right-handed neutrino mass matrix $`M`$ is proportional to 1, $`H_{23}`$ in $`H`$ may be comparable to $`H_{33}`$ and $`H_{22}`$ since $`\mathrm{tan}^2\theta _{23}1`$ in the MNS matrix. Also, we may need some conspiracy between the Yukawa coupling and the Majorana masses of the right-handed neutrino in order to realize the large mixing angles of the atmospheric and solar neutrino oscillation in $`H_2`$. However, both possibilities are not rejected phenomenologically at present, and thus, we explore the both cases.

In order to reconstruct the seesaw model, the light-neutrino mass matrix also has to be fixed. We use $`\mathrm{\Delta }m_{32}^2=3\times 10^3`$ eV$`^2,`$ $`\mathrm{\Delta }m_{21}^2=4.5\times 10^5`$ eV$`^2,`$ $`\mathrm{tan}^2\theta _{23}=1`$ and $`\mathrm{tan}^2\theta _{12}=0.4`$ corresponding to the LMA solution. Since the bound on the angle $`\theta _{13}`$ is quite stringent, our results depend very weakly on its actual value. We fix $`\mathrm{sin}\theta _{13}=0.1`$ and the CP phases in the MNS matrix $`\delta =\pi /2`$. We study both the normal and the inverse hierarchy for the light-neutrino mass spectrum, since the neutrino oscillations do not discriminate between these two cases. As input parameters, we then have the lightest effective neutrino mass $`m_1`$ (or $`m_3`$ for inversely ordered neutrinos), which we generate in the range $`(10^40.3)`$ eV, the two low scale Majorana phases $`\varphi _{1,2}`$ and the matrix $`H`$, which we generate randomly.

In Fig. 1 we present Br($`\tau \mu \gamma `$) for the ansatz $`H_1`$, assuming the normal or the inverted hierarchy for the light-neutrino mass spectrum. The horizontal axis is the lighter stau mass $`m_{\stackrel{~}{\tau }_1}`$, and the other supersymmetry-breaking parameters are determined by choosing the SU(2) gaugino mass to be 200 GeV, $`A_0=0`$, $`\mu >0`$, and $`\mathrm{tan}\beta =10`$. We see from these figures that behaviors of the branching ratios are similar for the normal and inverted hierarchies for the light-neutrino mass spectrum. In our parametrization the branching ratio is determined mainly by $`H`$ and the sparticle mass spectrum. In this figure the maximum value of the prediction is fixed by the perturbation condition, and it is larger than the current experimental bound. We sample the parameters in $`H_1`$ randomly in the range $`10^2<a,b,c,|d|<10`$, with distributions that are flat on a logarithmic scale, and require the Yukawa coupling-squared to be smaller than $`4\pi `$, so that $`Y_\nu `$ remains perturbative up to $`M_G`$.

In Fig. 2, we present Br($`\tau e\gamma `$) for the ansatz $`H_2`$, assuming the normal or the inverted hierarchy for the light-neutrino mass spectrum. The input parameters are the same in Fig. 1. The significant difference between the normal and inverted hierarchies for the light-neutrino mass spectrum does not exit, and the maximum values of the prediction is determined by the perturbation condition, again. This means that both searches for $`\tau \mu \gamma `$ and $`\tau e\gamma `$ are important and we have possibility to observe either of them in future experiments.

Fig. 1 and Fig. 2 show that the prediction of the LFV tau decay modes may be almost independent of the light-neutrino mass matrix. The search for the LFV tau decay is one of important pieces to reconstruct the seesaw model.

Finally, we comment on the other LFV tau decay processes. The LFV tau decay modes to three charged leptons are suppressed in the supersymmetric models if $`R`$ parity is conserved. In the almost parameter space the photonic penguin diagram dominates over the other diagrams.222 If $`\mathrm{tan}\beta `$ is very large and pseudoscalar Higgs mass is light, the Higgs penguin may give a sizable contribution .  Then, we can get the following approximate formulas,

$`{\displaystyle \frac{\mathrm{Br}(\tau 3\mu )}{\mathrm{Br}(\tau \mu \gamma )}}`$ $``$ $`{\displaystyle \frac{\mathrm{Br}(\tau e2\mu )}{\mathrm{Br}(\tau e\gamma )}}{\displaystyle \frac{1}{440}},`$ (19)
$`{\displaystyle \frac{\mathrm{Br}(\tau \mu 2e)}{\mathrm{Br}(\tau \mu \gamma )}}`$ $``$ $`{\displaystyle \frac{\mathrm{Br}(\tau 3e)}{\mathrm{Br}(\tau e\gamma )}}{\displaystyle \frac{1}{94}}.`$ (20)
The branching ratio for $`\mu e^+e^{}(ee^+e^{})`$ is larger than that to $`\mu \mu ^+\mu ^{}(e\mu ^+\mu ^{})`$, because the phase space is larger.

Acknowledgments
This work is partially supported by the Grant-in-Aid for Scientific Research from the Ministry of Education, Science, Sports and Culture of Japan (No. 13135207 and No. 14046225).


text:
# 1 Introduction

## 1 Introduction

Accreting black holes radiate in the two main spectral states which we refer to later as hard and soft (see Fig. 1). The hard state is characterized by a power-law–like spectrum which abruptly cuts off at $``$ 100 keV. The energy output is dominated by the 100 keV photons . Such a spectrum is interpreted as Comptonization by thermal electrons in the inner hot disk or active magnetic corona above the accretion disk (e.g. ). A weak MeV tail observed in Cyg X-1 is probably a signature of non-thermal electrons in the source .

In the soft state the emission is dominated by the black-body–like component peaking at a few keV with a power-law–like component above 30 keV extending up to 1 MeV or even higher . These spectra cannot be fitted with thermal Comptonization models and require the radiating electrons to have a significant non-thermal fraction .

## 2 Spectral Modeling

### 2.1 Cygnus X-1

A weak MeV tail observed in Cyg X-1 in its hard state already shows us that the emitting electrons cannot be purely thermal (see left panel of Fig. 2). An obvious generalization is to assume that there is a non-thermal tail in the electron distribution that is produced by some acceleration process (e.g. ). The soft state data (even below 1 MeV) could not be fitted at all with thermal models. Modeling of the spectral transitions with a generalized hybrid thermal/non-thermal model eqpair  predicted a stronger power-law–like tail in the soft state extending up to 10 MeV. CGRO/COMPTEL observations confirmed the existence of this tail: the decrease of the hard X-ray luminosity was accompanied by the increase of the soft $`\gamma `$-ray luminosity . In the context of the hybrid Comptonization model, the power-law is a result of single Compton scattering off non-thermal population of electrons (see right panel in Fig. 2). A cutoff at about 10 MeV (depending on the compactness of the source) should appear in the spectrum due to the absorption of the $`\gamma `$-rays by softer photons resulting in pair production.

### 2.2 GRS 1915+105

In spite of the fact that microquasar GRS 1915+105 show dramatic variability pattern, almost all its hard X-ray spectra are remarkably similar. During eight (out of nine) observations with the CGRO/OSSE the source showed a simple power-law–like spectrum in the 50–500 keV band with photon index $`\mathrm{\Gamma }3`$ . Only in one occasion (when the X-ray flux was very high), the hard X-ray spectrum was much harder $`\mathrm{\Gamma }2.3`$ and the flux was low. We note here that in all observations the source has a much softer spectrum than the normal hard state of Cyg X-1, i.e. it was always in the soft state. There is no signature of the high-energy cutoff in the data. The eqpair model gives a good description of the data (see Fig. 3) indicating that about 10-20% of the total power goes to accelerate non-thermal electrons. The C/$`\chi `$-state  differs, however, from the B/$`\gamma `$-state in that the 20-200 keV tail is produced by thermal Comptonization in the former and by non-thermal Compton scattering in the later.

## 3 Old and New Alternatives

### 3.1 The “standard” model

The black body looking soft component is associated with the optically thick accretion disk by most of the researchers. The broad-band spectra are often fitted by the “standard” model consisting of a disk-blackbody (soft component) and a power-law (hard). There are numerous problems with such modeling. First, a black body is a bad representation of the spectrum expected from the accretion disk (e.g. ). Real data also show that the soft bumps in the Cyg X-1 soft state  and GRS 1915+105  cannot be fitted by a black body (or multicolor disk). Thermal Comptonization of a blackbody is a much better description of these spectra. Second, a power-law, even exponentially folded, is a very bad representation of the Comptonization spectra. At the lower end, Comptonization spectrum cuts off below the seed photon energy while a power-law has no break there. The normalization of the blackbody thus can be underestimated by a large factor (see e.g. ). The conclusions (e.g. variations of the inner disk radius) resulting from fitting the data with this “standard” model thus should be taken with a grain of salt.

At the higher end, the (thermal) Comptonization spectrum has a much sharper cutoff than an exponentially folded power-law. This difference in the spectral shape is important when we model the broad-band spectra from accreting black holes with the later model adding a Compton reflection component (model pexrav  from XSPEC) since the amplitude of the reflection component strongly depends on the assumed shape of the underlying continuum. Thus, we would advise not to use pexrav when modeling Comptonization spectra close to the black body or to the high energy cutoff.

### 3.2 X-rays from the jet?

A very interesting correlation between radio and X-ray fluxes has been discovered recently . There are two possible origins of this correlation. One is that the level of X-ray emission is related to the rate of ejection of radio-emitting clouds, forming a compact jet (e.g., ). Another is that the X-ray emission of black hole binaries is dominated by the synchrotron emission of the jet .

We note here there are many strong arguments against the second interpretation. The broad-band X/$`\gamma `$-ray spectra of black hole binaries in the hard state are very well modeled by thermal Comptonization and Compton reflection (e.g., ). The presence of reflection implies that the X-ray emission is not strongly beamed away from the disk. The thermal-Compton origin of the primary X-ray emission is strongly supported by a remarkable uniformity of the both energy and shape of the high-energy cutoffs of black hole binaries in the hard state observed by OSSE . This cutoff is naturally accounted for by thermostatic properties of thermal Comptonization as well as $`e^\pm `$ pair production (e.g., ), as it corresponds to the transition to relativistic temperatures. At higher temperatures, cooling becomes extremely efficient and copious pair production starts. This reduces the energy available per particle leading to the temperature decrease. On the other hand, $`m_\mathrm{e}c^2`$ plays no particular role for non-thermal synchrotron emission (cf. variable cutoff energy during flares in blazars). Thus, accounting for the observed cutoff energies requires fine-tuning of product of the square of the maximum electron energy and the magnetic field strength in the non-thermal synchrotron models. In addition, the jet model has problems reproducing the actual shape of the cutoff. For example, the synchrotron model of the hard state of Cyg X-1 (see fig. 3a in ) when matched to the 100 keV flux overestimates the 1 MeV flux (see ) by a factor of 8.

An additional evidence against a substantial part of X-rays being non-thermal synchrotron is provided by spectral variability. In the case of Cyg X-1, the ASM/BATSE data show spectral pivoting around $`40`$–50 keV (see  and Fig. 4). The characteristic variations of the power-law slope $`\mathrm{\Delta }\mathrm{\Gamma }`$ from those data is $`0.2`$–0.3. This power-law spectral variability extended to 15 GHz would imply a huge variability of the radio flux, by several orders of magnitude. However, the range of the variability of the 15 GHz flux correlated with the ASM flux is by a factor of several, basically the same as the range of the variability of the ASM flux itself . If both radio and X-rays were due to non-thermal synchrotron emission, their observed variability pattern should yield the rms in X-rays virtually independent of energy. This is clearly in strong disagreement with the data shown in Fig. 4, and, in particular, with the ASM 1.5–3 keV flux being strongly anticorrelated with the 100–300 keV flux from BATSE .

The amplitude of Compton reflection and the iron line imply that dense and rather cold material occupies a solid angle $`\pi `$ as viewed from the X-ray source. Smearing of these components and the observed correlations with the spectral slope  clearly identifies the reflector with the accretion disk and implies the origin of the continuum emission within 30–100 gravitational radii from the black hole which would also be consistent with rapid X-ray variability. The above arguments also rule out the dominant contribution of the nonthermal Compton emission  from the jet to the observed hard state spectra of black holes.

All these arguments strongly support the interpretation of the correlated radio emission as being due to ejection of clouds from the X-ray source which could be similar to coronal mass ejections (CME) observed at the Sun. The recently discovered strong X/$`\gamma `$-rays flares from Cyg X-1  could be the extremes of such an activity. The right panel in Fig. 4 shows the flaring activity of Cyg X-1 in April 1999. The episode D–E shows strong flare in the 20–100 keV band with a weaker activity above 100 keV and no detectable signal above 300 keV. This indicates that there exist at least two independent spectral components (see ): one could be related to the inner hot disk or the magnetized corona, while another to the base of the jet.

### 3.3 Bulk motion Comptonization

The power-law like spectra of black holes in the soft state were interpreted as resulting from bulk motion Comptonization in the converging flow . The specific feature of that model is a cutoff at $`100`$–200 keV. (We note here that the XSPEC version of the model bmc has no cutoff built in.) The data show no signatures of the cutoff at least up to 500 keV in GRS 1915+105  and up to 10 MeV in Cyg X-1 . This supports their non-thermal origin and strongly rules out a significant contribution from the bulk motion Comptonization. (See further discussion in .)

## Acknowledgments

This work was partly supported by the Academy of Finland and grants from the Polish Committee for Scientific Research (5P03D00821, 2P03C00619p1,2) and the Foundation for Polish Science.


text:
# Particle Correlations at RHIC – Scrutiny of a Puzzle

## Abstract

We present calculations of two-pion and two-kaon correlation functions in relativistic heavy ion collisions from a relativistic transport model that includes explicitly a first-order phase transition from a thermalized quark-gluon plasma to a hadron gas. We compare the obtained correlation radii with recent data from RHIC. The predicted $`R_{\mathrm{side}}`$ radii agree with data while the $`R_{\mathrm{out}}`$ and $`R_{\mathrm{long}}`$ radii are overestimated. We also address the impact of in-medium modifications, for example, a broadening of the $`\rho `$-meson, on the correlation radii. In particular, the longitudinal correlation radius $`R_{\mathrm{long}}`$ is reduced, improving the comparison to data.

Two-particle correlations at small relative momenta have been predicted to be particularly sensitive to a phase transition from quark-gluon matter to hadronic matter . For a first-order phase transition, larger hadronization times were expected to lead to considerably enhanced correlation radii, characterizing the space-time extension of the particle-emitting source, compared to, for example, a purely hadronic scenario. The radii should also depend on the critical temperature $`T_c`$, the latent heat, the initial specific entropy density or the initial thermalization time of the quark-gluon phase.

Here, we discuss relativistic transport calculations at RHIC energies that describe the initial dense stage by hydrodynamics  and the later more dilute stages by microscopic transport  of the particles. The two models are matched at the hadronization hypersurface . In the hadronic phase the particles are allowed to rescatter and to excite resonances based on cross sections as measured in vacuum. For the initial dense (hydrodynamical) phase a bag model equation of state exhibiting a first-order phase transition is employed. Hence, a phase transition in local equilibrium that proceeds through the formation of a mixed phase, is considered. The details of this relativistic hybrid transport model can be found elsewhere .

We first briefly summarize the main conclusions obtained in previous work . Then, we show in detail the results for pions when calculating explicitly the correlation functions from the source function of the transport model. These results are subsequently compared to calculations that take in-medium modifications into account and to experimental data for central Au+Au collisions at $`\sqrt{s}_{NN}=130`$GeV .

Model calculations have demonstrated that the freeze-out hypersurfaces of pions extend to rather large radii and times compared to the size of the mixed phase . In this hadronic phase many soft collisions take place which hardly modify the single-particle spectra but have a strong impact on the correlation functions that measure the final freeze-out state. Studying this in detail leads to the following conclusions : (i) The dissipative hadronic phase leads to a rather large duration of emission. (ii) The $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ ratio, thought to be a characteristic measure of this emission duration, increases with transverse momentum. (iii) The specific dependencies of the interferometry radii on the QGP properties are rather weak due to the dominance of the hadronic phase. This even results in qualitative differences if calculations with and without this subsequent hadronic phase are compared (dependence on the critical temperature) .

The correlations of kaons provide a severe test of the pion data and have several advantages . In particular, the kaon density is much lower than the pion density . Hence, multiparticle correlations that might play a role for the pions are of minor importance for the kaons. Also, the contributions from long-lived resonances are under better control for kaons. The $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ ratio for kaons is shown in Fig. 1. The sensitivity to $`T_c`$ and to the specific entropy density (SPS vs. RHIC) is enhanced at larger transverse momenta ($`K_T1`$GeV/c). This enhanced sensitivity is also driven by a strong increase of the direct emission component (from the phase boundary) at high $`K_T`$ as shown in Fig. 2. More and more kaons (up to $`30\%`$) escape the initial stages (unperturbed by the hadronic phase).

The pion correlation parameters $`R_\mathrm{o},R_\mathrm{s},R_\mathrm{l}`$and $`\lambda `$ are obtained from the explicit calculation of the correlation functions  in the respective transverse momentum bins and subsequent fitting of these correlation functions to a Gaussian form of the correlator $`C_2=1+\lambda \mathrm{exp}(R_\mathrm{o}^2q_\mathrm{o}^2R_\mathrm{s}^2q_\mathrm{s}^2R_\mathrm{l}^2q_\mathrm{l}^2)`$. The results of fitting the 3-dimensional pion correlation functions (as obtained from the transport calculations (RHIC initial conditions, $`T_c160`$MeV) + correlation after burner (by Pratt) ) are shown, together with the experimental STAR and PHENIX data , in Figure 3. Even for a first-order phase transition scenario, the calculated interferometry radii (T) are not unusually large. Finite momentum resolution (f.m.r.) reduces the true correlation radii and the $`\lambda `$ intercept parameter, in particular, at higher $`K_T`$ . This is corrected for in the experimental analysis.

While the $`R_{\mathrm{side}}`$ radii show agreement the $`R_{\mathrm{out}}`$ radii are too large compared to the data. However, these pion radii are considerably smaller than the corresponding radii obtained from the coordinate-space points and using expressions for the Gaussian radius parameters based on a saddle-point integration over the source function . Only the calculations as shown in Figure 3 which are obtained from the complete calculation and the performed fits should be compared to data. The $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ ratio is also larger than unity for the fitted values and confirms the so-called HBT-puzzle, i.e., the RHIC data from STAR and PHENIX  indicate a decreasing ratio with $`K_T`$ (even below 1) and the calculations show only ratios larger 1. We remark that at the CERN-SPS (see Refs. in ) the $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ ratio seems to increase with $`K_T`$, being larger than 1, implying a real qualitative change of the reaction dynamics from SPS to RHIC energies. A discussion of theoretical assumptions and uncertainties, experimental corrections, and possible solutions of this problem (including strongly opaque sources or scenarios with strong supercooling leading to spinodal instabilities) is provided in .

In-medium modifications may alter the correlation radii due to their impact on the mean free path via modified cross sections at finite temparatures or densities ($`T_{\mathrm{f}.\mathrm{o}.}100120`$MeV). Increasing the width of the $`\rho `$ from the vacuum value of $`\mathrm{\Gamma }_\rho =150`$MeV to $`\mathrm{\Gamma }_\rho ^{\mathrm{i}.\mathrm{m}.}=350`$MeV enhances the opacity of the system (more $`\pi \pi \rho `$ processes), with a particular reduction of the radius $`R_{\mathrm{long}}`$, improving the comparison to data (Fig. 3).

We have demonstrated that the kaon interferometry measurements, in particular at high $`K_T`$, will provide an exellent probe of the space-time dynamics (close to the phase boundary). In addition they represent a severe test of the pion correlations and may help to better understand the HBT-puzzle , i.e., the difference in the $`K_T`$ dependence of the $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ ratio between the model predictions and the experimental RHIC data. The differences are due to the $`R_{\mathrm{out}}`$ radii (which are larger in the model calculations) while the $`R_{\mathrm{side}}`$ radii seem to be described reasonably. An in-medium broadening of the $`\rho `$-meson reduces in particular $`R_{\mathrm{long}}`$, improving the comparison to data.

Acknowledgments:  SS is supported by the Alexander von Humboldt Foundation through a Feodor Lynen Fellowship and DOE Grant No. DE-AC03-76SF00098.


text:
# Generic slow relaxation in a stochastic sandpile

## Abstract

Simulations of a stochastic fixed-energy sandpile in one and two dimensions reveal slow relaxation of the order parameter, even far from the critical point. The decay of the activity is best described by a stretched-exponential form. The persistence probability (for a site not to have toppled up to time $`t`$), also exhibits stretched-exponential relaxation. The results suggest a connection between sandpile models and structural glasses.

PACS: 05.40.-a, 05.65.+b, 64.70.Pf, 05.50.+q

Email address: dickman@fisica.ufmg.br

Sandpile models have attracted much interest in recent years, as paradigms of scale-invariance in the apparent absence of tuning parameters , and as intriguing examples of absorbing-state phase transitions . Most studies of sandpiles have focused on the scale-invariant stationary state under slow driving , or on scaling properties in the vicinity of the absorbing-state transition . Relatively little attention has been given to the dynamic properties of sandpiles away from the critical point.

A central feature of sandpile models is the presence of a conserved field, the density of particles. This field couples to the activity density, which is the order parameter. When, as in the case of FES, the conserved field is frozen in the absence of activity, the critical behavior is expected to fall in a universality class distinct from that of directed percolation (DP) . (In the absence of such a conservation law, DP is generic for absorbing-state phase transitions . Relaxation of the order parameter to its quasi-stationary value in a sandpile at the critical density is best characterized as a power-law , but with certain anomalies in the one-dimensional case . Besides having an effect on the critical behavior, it is reasonable to expect the conservation law to modify relaxation properties away from the critical point as well.

In simple models exhibiting an absorbing-state phase transition, such as the contact process , relaxation to the stationary state is expected to be monotonic and exponential, away from the critical point . The simulation results reported below show that in sandpiles, off-critical relaxation is considerably slower, following a stretched-exponential or, in certain cases, algebriac decay. Stretched-exponential functions have been reported for avalanche size distributions in experiments on sand and rice piles , and in granular avalanche models , but not, to my knowledge, in the context of sandpile relaxation dynamics.

I study a variant of Manna’s sandpile , defined on a lattice of $`L^d`$ sites (in $`d`$ dimensions), with periodic boundaries. The configuration is specified by the number of particles $`z_i=0,1,2,\mathrm{}`$ at each site; sites with $`z_i2`$ are said to be active. A Markovian dynamics is defined by the toppling rate, which is unity for all active sites, and zero for sites with $`z_i<2`$. When a site $`i`$ topples, it sends two particles to adjacent sites ($`z_iz_i2`$); the particles move independently to randomly chosen nearest neighbors $`j`$ and $`j^{}`$ ($`j`$, $`j^{}\{i+1,i1\}`$ in the one-dimensional case). Thus $`j=j^{}`$ with probability $`1/2d`$. The dynamics conserves the number of particles, $`N`$, which is fixed by the initial configuration.

The system evolves via a sequential dynamics: the next site to topple is chosen at random from a list of active sites; the time increment associated with a toppling is $`\mathrm{\Delta }t=1/N_A`$, where $`N_A`$ is the number of active sites just prior to the event. Initial configurations are generated by distributing $`\zeta L^d`$ particles randomly over the lattice, yielding an initial distribution that is spatially homogeneous and uncorrelated. Previous studies confirm the existence of a continuous phase transition from an absorbing to an active phase at a particle density $`\zeta _c=0.94885`$ in one dimension , and 0.71695 in 2-$`d`$ . For $`\zeta <\zeta _c`$, the stationary value of the order parameter (the density of active sites) is zero.

I begin with the simplest case, relaxation of the order parameter $`\rho `$ in one dimension, for $`\zeta <\zeta _c`$. For $`\zeta =0.5`$, far below the critical value, the pattern of relaxation is essentially the same for systems of 1000 or more sites. Figure 1 shows that $`\rho (t)`$ (obtained from averages over $`5\times 10^5`$ realizations of a system of 5000 sites), is well described by a stretched exponential, $`\rho (t)\mathrm{exp}[(t/t_0)^\beta ]`$ with $`\beta =0.45`$ ($`t_0`$ represents a characteristic timescale for relaxation). In this and subsequent analyses, the exponent $`\beta `$ is determined using the criterion of zero curvature, in the asymptotic region of the graph of $`\mathrm{ln}[\rho (t)\overline{\rho }]`$ versus $`t^\beta `$ ($`\overline{\rho }`$ is the asymptotic activity density, which is of course zero for $`\zeta <\zeta _c`$). A rough estimate of the uncertainty in $`\beta `$ values is $`\pm 0.02`$.

Next I examine the relaxation for $`\zeta =0.9`$, much nearer the critical point. There are now considerable finite-size effects and much larger systems (up to $`10^5`$ sites) are required to observe the asymptotic behavior. For small systems $`\rho (t)`$ appears to decay faster than a stretched exponential, but as $`L`$ is increased, the long-time slow relaxation develops. For $`L=10^5`$, $`\rho (t)`$ again follows a stretched exponential, with $`\beta =0.39`$. It should be emphasized that while the simulation data cannot be interpreted as proving stretched-exponential relaxation, they do serve to rule out definitively both exponential and power-law relaxation of $`\rho (t)`$ in the subcritical regime.

In the supercritical regime, the relaxation of the order parameter to its stationary value $`\overline{\rho }`$ is nonmonotonic. $`\rho (t)`$ decays rapidly at first, and then slowly approaches the stationary value from below. For $`\zeta =1`$, I find $`\overline{\rho }=0.118222(8)`$, with $`\mathrm{\Delta }\rho \overline{\rho }\rho (t)`$ decaying asymptotically as a power-law, $`t^{0.54}`$ (see Fig. 2). The initial decay is again well described by a stretched exponential, with an exponent $`\beta =0.28`$. For somewhat higher densities ($`\zeta =1.2`$ and 1.3), the asymptotic approach to $`\overline{\rho }`$ (from below), is again via a power law, with an exponent of 0.5 - 0.52.

I also studied a variant of the stochastic sandpile introduced in Ref. , in which the toppling rate at a site with $`z`$ particles is $`z(z1)`$. (In this case $`\zeta _c=0.9493(2)`$.) For $`p=1/2`$, the decay to $`\rho =0`$ is described with high precision by a stretched exponential with $`\beta =0.475(25)`$. In the supercritical regime the approach to the stationary value is nonmonotonic, similar to that observed above. Unlike the constant-rate model, however, the relaxation is well described by an expression consisting of two stretched exponentials:

$$\rho (t)=\overline{\rho }+A_1\mathrm{exp}[\lambda _1t^{1/2}]A_2\mathrm{exp}[\lambda _2t^{1/2}],$$
(1)
with $`A_1`$, $`A_2`$, $`\lambda _1`$, and $`\lambda _2`$ all positive constants ($`A_1A_2`$, and $`\lambda _1\lambda _2`$).

Slow relaxation appears to be robust under changes in the toppling rate. Of equal interest is the nature of relaxation in two or more dimensions. At a density of $`\zeta =0.5`$ (well below $`\zeta _c=0.71695(5)`$), studies using lattices of $`L^2`$ sites, with $`L=40`$, 80, and 160, reveal that the order parameter decays as a stretched exponential, but with an exponent $`\beta =0.81`$. Nearer the critical point ($`\zeta =0.7`$), finite-size effects are prominent, as was found in the one-dimensional case. Studies using system sizes $`L`$ of up to 1280 again confirm stretched-exponential decay with $`\beta =0.81`$, as shown in Fig. 3.

The relaxation above the critical density (in a system of size $`L=160`$ at $`\zeta =0.75`$; $`\overline{\rho }=0.04995`$), is nonmonotonic, as in one dimension. The excess density is again well described by the difference between two stretched exponentials, as in Eq. (1), but in this case the exponent $`\beta `$ associated with the early decay is about 0.4, while the long-time approach to $`\overline{\rho }`$ (from below) is characterized by an exponent of 0.8 (see Fig. 3, inset).

As we have seen, relaxation of the order parameter in the stochastic sandpile is characterized by stretched-exponential (or in some cases, power-law) functions. Recently, O’Donoghue and Bray  demonstrated stretched-exponential decay of the persistence probability, i.e., that a given site has never been visited by a diffusing particle, in certain one-dimensional reaction/diffusion processes. This suggests a study of persistence in the sandpile. In fact, since a site that is initially below threshold for toppling must be visited by another diffusing particle or particles before it can topple, persistence seems particularly relevant to sandpile relaxation dynamics. (Studies of dynamic critical exponents for persistence in conserved lattices gases, which are closely related to sandpiles, were reported by Lübeck ; here, however, we focus on the dynamics away from the critical point.)

In studying persistence in a sandpile, it appears useful to group sites according to their initial occupation number $`z`$. Let $`p(t;z)`$ denote the persistence probability of sites whose initial height is $`z`$. For $`z=0`$ or 1, I define $`p(t;z)`$ as the probability that a site has not toppled up to time $`t`$. For $`z>2`$, however, such a definition is not very interesting, since (in the constant-rate model), the probability not to have toppled up to time $`t`$ is simply $`e^t`$. For $`z2`$, therefore, $`p(t;z)`$ is defined as the probability not to have toppled a second time. I study the one-dimensional fixed-rate model with particle densities $`\zeta =0.5`$, 0.9, and 1, as above.

For $`\zeta <\zeta _c`$, all activity ceases after a finite time, so that the persistence probability approaches a nonzero value $`\overline{p}(z)`$ at long times. For $`\zeta =0.5`$, the asymptotic relaxation of $`p(t;z)`$ to $`\overline{p}(z)`$ again follows a stretched exponential; the best estimates for the exponent $`\beta `$ are 0.42 for $`z=0`$ and $`\beta =0.39`$ for $`z=1`$, 2 and 3. For $`\zeta =0.9`$, a study using $`L=10^4`$ yielded $`\beta 0.44`$, 0.47, 0.49, and 0.41 for $`z=0`$, 1, 2, and 3, respectively. Finally, for $`\zeta =1`$, the asymptotic decay of the persistence probability (studied on a ring of $`2\times 10^4`$ sites) follows a stretched exponential with $`\beta `$ values of 0.42, 0.41, 0.38, and 0.45 for the various $`z`$ values. The persistence probabilitites for $`\zeta =0.9`$ and $`\zeta =1`$ are shown in Fig. 4. (In all of these studies, the exponents for $`z=2`$ and 3 are less certain, due to poorer statistics, and to the fact that stretched-exponential behavior sets in later than for $`z=0`$ or 1. I find no clear trend in the estimated exponent value, either with $`z`$ or with $`\zeta `$.) The foregoing results can be summarized as indicating that in one dimension, the asympotic relaxation of the persistence probability $`p(t;z)`$ follows a stretched exponential with an exponent $`\beta `$ of about 0.42.

The simulation results indicate that the relaxation of the order parameter and of the persistence probability, in a stochastic sandpile away from its critical point, generically exhibits stretched-exponential scaling. The exponent $`\beta `$ generally takes values in the neighborhood of 0.4 in the one-dimensional case, while in two dimensions (for which only order-parameter relaxation was studied) the value is about 0.8. Stretched-exponential relaxation may be understood as a consequence of the depletion of active elements, leading to a decreasing relaxation rate, $`(1/\rho )|d\rho /dt|`$, as time goes on. The larger exponent in two dimensions is plausible, given the larger number of paths in configuration space. Detailed explanations of the streteched-exponential form, and of specific exponent values, remain as open challenges.

The observation of stretched-exponential relaxation suggests a connection with glassy dynamics. Superficially, the sandpile model seems to have little connection with dynamics of a dense fluid, but the two problems are related at a more abstract level. In dense fluids, not all the empty space in the system is available for particle movement, leading to a highly cooperative dynamics . Some larger grouping of “voids” is required for relaxation on scales beyond that of the local cage. (For simplicity, I frame the analogy in terms of a hard-sphere fluid, in which the relevant variable is density not temperature.) Similarly, in a sandpile not all particles are available for movement (and relaxation): only those with companions in the same cell are mobile. These observations suggest a parallel between particles in the sandpile and parcels of unoccupied space in a dense fluid. The dynamic arrest in the fluid, when the density approaches an (apparent) transition value, corresponds to the absorbing-state phase transition in the sandpile. The associated order parameter - mobile free volume in the fluid - vanishes in the absorbing phase. (The total unoccupied volume is of course conserved, but the fraction that is mobile or “active” is not.) From this vantage, the liquid-glass transition corresponds to a dynamic transition, akin to directed percolation, of free volume, as opposed to the static percolation transition suggested some time ago by Cohen and Grest . While the sandpile models studied to date seem too simplistic to capture the dynamics of a dense fluid near the glassy state, the results of the present study suggest that it would be worthwhile pursuing this analogy by developing somewhat more realistic particle models, and studying diffusion and the response to external driving, in the hope of gaining further insight into real glasses.

Acknowledgments

I thank Kim Christensen and Francisco F. Araujo Jr. for helpful discussions. This work was supported by CNPq, CAPES, and FAPEMIG, Brazil.

FIGURE CAPTIONS

FIG. 1. Relaxation of the order parameter in the one-dimensional sandpile, $`\zeta =0.5`$, $`L=5000`$. Inset: semi-log plot of the same data.

FIG. 2. Relaxation of the order parameter in the one-dimensional sandpile, $`\zeta =1`$, $`L=2\times 10^4`$. Inset: log-log plot of the excess $`\mathrm{\Delta }\rho =\rho (t)\overline{\rho }`$.

FIG. 3. Relaxation of the order parameter in the two-dimensional sandpile, $`\zeta =0.7`$, $`L=1280`$. Inset: $`\mathrm{ln}|\mathrm{\Delta }\rho |`$ versus $`\mathrm{ln}t`$ for $`\zeta =0.75`$; points: simulation; solid curve: fit using two stretched exponentials as described in text.

FIG. 4. Persistence probabilities $`p(t;z)`$ for $`z=0`$, 1, 2, and 3 (upper to lower) in the one-dimensional model with $`\zeta =1`$. Inset: a similar plot for $`\zeta =0.9`$.


text:
# Critical load and congestion instabilities in scale-free networks

## Abstract

We study the tolerance to congestion failures in communication networks with scale-free topology. The traffic load carried by each damaged element in the network must be partly or totally redistributed among the remaining elements. Overloaded elements might fail on their turn, triggering the occurrence of failure cascades able to isolate large parts of the network. We find a critical traffic load above which the probability of massive traffic congestions destroying the network communication capabilities is finite.

Complex heterogeneous connectivity patterns have been recently identified in several natural and technological networks strogatz ; bara02 ; doro . The Internet and the World-Wide-Web (WWW) networks, where nodes represent routers or web pages and edges physical connections or hyper-links, appear to have a topology characterized by the presence of “hubs” with many connections to peripherical nodes. Empirical evidence recently collected shows that this distinctive feature finds its statistical characterization in the presence of heavy-tailed degree distributions bara00 ; fal ; broder00 ; gov00 ; gcalda ; gov01 ; caida ; romu01 . In the Internet, for instance, the statistical analysis reveals that the degree distribution $`P(k)`$, defined as the probability that any node has $`k`$ links to other nodes, is well approximated by a power-law behavior $`P(k)k^\gamma `$, with $`\gamma 2.2`$ fal ; gov00 ; romu01 . This makes the Internet a capital example of the recently identified class of scale-free (SF) networks bara02 . The statistical physics approach has been proved to be a valuable tool for the study of complex networks, and several interesting results concerning dynamical processes taking place on complex networks have been recently reported. In particular, the absence of the percolation newman00 ; havlin01  and epidemic pv01a ; virusreview  thresholds in SF networks has a large impact because of its potential practical implications. The absence of the percolating threshold, indeed, prompts to an exceptional tolerance to random damages barabasi00 . This is a property that assumes a great importance in communication networks, guaranteeing the connectivity capabilities of the system.

Percolation properties of SF networks refer only to the static topological connectivity properties newman00 ; havlin01 . On the other hand, in the Internet and other communications networks, many instabilities are due to traffic load congestions labo1 ; labo2 . The traffic load carried on the failing nodes or connections is automatically diverted to alternative paths on the networks and instabilities can spread from node to node by an avalanche of traffic congestions and overloads. For instance, route flaps have led to the transient loss of connectivity for large portions of the Internet. These instabilities are thus of a dynamical nature and depend on how information is routed and distributed in the network. The models proposed so far, however, deal with regular structures jap1 ; jap2 ; watts02  and do not take into account the complex topology of SF networks.

In this paper, we propose a simple model aimed at the study of failure cascades generated by the redistribution of traffic load by congested links or nodes in SF networks. We find that the system behavior depends on the average traffic load imposed to the network. Above a critical value of the average traffic load, a single failure has a finite probability of triggering a congestion avalanche affecting a macroscopic part of the network. The present analysis thus reveals the existence of a transition from a free phase to a congested one as a function of the amount of traffic carried by the network. Contrary to what happens for the static percolation transition havlin01 ; newman00 , loaded SF networks exhibit a finite threshold above which the system can develop macroscopic instabilities with respect to small damages if we consider the dynamics of the traffic carried on top of them. The results provided here represent a first step towards a more complete modeling of traffic instabilities in real communication networks.

In order to include the degree fluctuations of SF networks we shall use in the following the Barabási-Albert model bar99 . This is a stochastic growth model in which one starts from a small number $`m_0`$ of nodes and at each time step a new node is introduced. The new node is connected preferentially to $`m`$ old ones (for the simulations we used $`m=3`$) with a probability $`\mathrm{\Pi }(k_i)=k_i/_jk_j`$ proportional to the degrees $`k_i`$ of the nodes. The repeated iteration of this scheme gives as a result a complex network with a topological structure characterized by a power-law degree distribution $`P(k)=2m^2k^\gamma `$ with $`\gamma =3`$ and average degree $`k=2m`$. In principle, one might also consider a more general class of complex networks with variable power-law degree distributions bara02 ; doro .

To simulate the flow of data packets on SF networks, taking into account the load redistribution in case of damages, we need to specify the initial state of the network; i.e. the load of traffic flowing through each link. An estimate of such load, assuming that the routing takes place following the minimum path, is given by the total number of shortest paths between any two nodes in the network that pass through the node $`i`$. This magnitude is called betweenness or load new01 ; goh01  and has been recently studied in SF networks. This property of the network can also be defined in terms of links. In real systems, however, the amount of traffic carried by each link is a fluctuating quantity that depends on many variables such as number of users, routing agreements, and available bandwidth. For this reason, we associate to each link connecting the nodes $`i`$ and $`j`$ of the network a load $`\mathrm{}_{i,j}`$ drawn from a probability distribution that specifies the initial traffic load of the system. For simplicity, we have considered a uniform distribution $`U(\mathrm{})`$ for $`0<\mathrm{}<1`$, taking the form

$$U(\mathrm{})=\{\begin{array}{ccc}\frac{1}{2\mathrm{}},& \mathrm{}[0,2\mathrm{}]\hfill & \mathrm{if}\mathrm{}0.5\hfill \\ & & \\ \frac{1}{2(1\mathrm{})},& \mathrm{}[2\mathrm{}1,1]\hfill & \mathrm{if}\mathrm{}0.5\hfill \end{array}.$$

This uniform distribution implies that the minimum initial load carried by a link is bounded by a nonzero value for an average load $`l>0.5`$, which means that there will be no links with load smaller than this lower bound. The results reported in this paper were obtained with the initial distribution $`U(\mathrm{})`$. In order to test the universality of the critical behavior, we have also considered the distribution $`F(\mathrm{})=(\mathrm{}^11)(1\mathrm{})^{(\mathrm{}^12)}`$, $`\mathrm{}[0,1]`$, which allows the existence of links with a very small load, irrespective of the average load $`\mathrm{}`$ flowing through the system. Both distributions $`U(\mathrm{})`$ and $`F(\mathrm{})`$ yield the same qualitative results. Along with the load, we associate to each link the same capacity $`C`$ that, without loss of generality, we fix equal to one. This choice can be considered as a first approximation since the actual difference between line bandwidths in communication networks can be large. In this perspective, we consider as the most important source of heterogeneity the flow of different amounts of load through the network.

The dynamics of the model is defined by a simple threshold process. A link is selected at random and overloaded by raising its traffic. When the load carried by a link is $`\mathrm{}_{i,j}>C`$, i.e. when it exceeds the link’s capacity, the link is considered congested and the load it carries is diverted among its (not overloaded) neighboring links. This amounts to consider that the time scale of the local congestion is greater than the time scale characterizing the reorganization of the routing procedure. The redistribution of the load on its turn might provoke that other links become overloaded, thus triggering a cascade of failures. We have explored two physically different settings of the load redistribution rule. The first consists on equally distributing the load of a congested link among the non-congested neighboring links. We refer to it as the deterministic redistribution rule, respectively. The second case will be called random redistribution because when a link is overloaded, a random amount of load is redistributed to each of the remaining working links in its neighborhood. Finally we note that in the rare event in which the congested link has no active neighbors, its load can be equally shared among all the remaining working lines of the network or just be considered as lost from the network. This amounts to a conserved or dissipative redistribution rule. Many physical systems display criticality only when energy is conserved jen98 ; marro . In distributed networks such as the Internet, however, it is common to discard packets if there is not a route available at the moment. As we shall see in the following, the results do not depend qualitatively on the conserved nature of the traffic load.

We have performed large-scale numerical simulations by applying repeatedly the rules stated above on BA networks. The sizes of the networks used in the simulations range from $`N=5\times 10^3`$ nodes ($`15\times 10^3`$ links) to $`N=10^5`$ nodes ($`3\times 10^5`$ links). All numerical results have been obtained by averaging over 10 different networks and, at least, 100 different realizations of the initial load distribution.

In order to inspect the occurrence of dynamic instabilities, we construct the phase diagram of the system. The order parameter can be identified as the probability $`P_𝒢`$ of having a giant component $`𝒢`$ of connected nodes with size of the order of the network size. The giant component is defined as the largest component of the network made by nodes connected by active links, after the system has reached a stable state (when $`\mathrm{}_{i,j}<C`$ for all $`i`$ and $`j`$). The existence of a giant component implies that a macroscopic part of the network is still functional. If the giant component of the network is zero, the communication capabilities of the network are destroyed and a congestion of the order of the system size builds up. It is worth noticing that although we have determined the giant component size in terms of nodes and the dynamical rules of the model are expressed in terms of links, the results are completely equivalent since a connected node is defined as a node with at least one active link.

In Fig. 1 we plot the order parameter $`P_𝒢`$ as a function of the average load $`\mathrm{}`$. At low values of the average load, the network always reaches a stable state in which the number of isolated nodes is very small and with probability $`P_𝒢=1`$ the network has a giant component of connected nodes of the order of the system size. When increasing the load imposed on the network, the system starts to develop instabilities. In particular, above a critical load $`\mathrm{}_c^I`$, whose value depends on the model considered, with a finite probability the system evolves to a congested state without giant component of connected nodes; i.e. the largest set of connected active nodes has a density of order $`N^1`$. This implies a probability of having a giant component $`P_𝒢<1`$, which is decreasing as the load is progressively increased. At an average load $`\mathrm{}_c^{II}0.82`$ we get that $`P_𝒢=0`$, signalling that, with probability one, any instability will propagate until the complete fragmentation of the network. It is worth remarking that this scenario is rather different from the percolation one in which the probability of having a giant component is abruptly dropping from one to zero at the transition point. Here, the probability decays continuously to zero and we have a wide region of $`\mathrm{}`$ where the initial instability can trigger a destructive congestion with probability $`1P_𝒢`$. Fig. 2 illustrates the probability $`p(S)`$ that the isolated network has a size $`S`$ in the case that no giant component of connected nodes has survived. The distribution is rather peaked also at relatively small values of $`\mathrm{}`$, almost affecting the totality of the network.

The phase diagram obtained in Fig. 1 points out that the value of the average load at which $`P_𝒢=0`$ is relatively high. On the other hand, the value at which $`P_𝒢`$ is appreciably smaller than one is well below the theoretical capacity of the network measured as the capacity $`C=1`$ of the individual links ($`\mathrm{}_c^I0.150.3`$, see Fig. 1). This evidence defines a wide region of load values in which there is a small but finite probability that a small instability propagates through the system and could explain why in real communication networks congestions of diverse degree can be encountered from time to time. In Fig. 1 we also report curves obtained for the different dynamical rules defined previously and it is interesting to remark that the highest level of stability has been obtained for the deterministic and conservative version of the model ($`\mathrm{}_c^I0.3`$). Moreover, the dynamics of the transition does not depend on the distribution of initial loads and the network size $`N`$.

Another way to shed light on the congestion dynamics is to inspect the process of generation of the congestion instability. We define the size $`s`$ of a congestion burst or avalanche as the total number of simultaneously overloaded links. The cumulative distribution $`P(s)`$ of avalanches of size larger than $`s`$ for several values of the average load imposed on the network and four different system sizes have been plotted in Fig. 3 for the random and dissipative definition of the present model. The main plot is in log-linear scale, so that a straight line corresponds to a power law of the form $`p(s)s^1`$ for the probability of observing an avalanche of size $`s`$. Power-laws with exponent $`1`$ have been found for several characteristic features of Internet traffic such as latency times, queue lengths, and congestion lengths jap1 ; mag ; will . In the figure we focus on region close to the stable region $`\mathrm{}0.20`$, which means that the power-law behavior extends to values far from the instability transition. This fact confirms that it is not necessary that the network operates very close to a critical point in order to observe power-laws in the distribution of several quantities. The inset in Fig. 3 shows that the cumulative size of overloaded links also scales with the system size, the scaling dynamics, however, remaining the same. This may help understand why power-law distributions observed in real communication networks have been measured for different network sizes, i.e. both for local networks and for networks that extend to a very large scale.

In summary, we have introduced a simple threshold model aimed at the description of instabilities due to load congestion that takes into account the topological properties of SF networks. The results obtained point out that the network can freely handle traffic up to some critical average load $`\mathrm{}_c^I`$. Above this level the network faces partial congestions that start to build up local bottlenecks in various places and small instabilities might trigger macroscopic outages with a finite probability. Above a critical load value $`\mathrm{}_c^{II}0.82`$ any small instability leads to the whole network collapse. In the intermediate region of network load, the number of simultaneous line casualties follows a power-law resembling what has been observed in experimental studies of the Internet. We hope that our work will provide hints for accurate modeling of the Internet and the WWW large-scale traffic behavior.

###### Acknowledgements.

This work has been partially supported by the European Commission - Fet Open project COSIN IST-2001-33555. R.P.-S. acknowledges financial support from the Ministerio de Ciencia y Tecnología (Spain).


text:
# Lasing from a single quantum wire

## Abstract

A laser with an active volume consisting of only a single quantum wire in the 1-dimensional (1-D) ground state is demonstrated. The single wire is formed quantum-mechanically at the T-intersection of a 14 nm Al0.07Ga0.93As quantum well and a 6 nm GaAs quantum well, and is embedded in a 1-D single-mode optical waveguide. We observe single-mode lasing from the quantum wire ground state by optical pumping. The laser operates from 5 to 60 K, and has a low threshold pumping power of 5 mW at 5 K.

A quantum wire laser is a novel semiconductor laser based on a 1-dimensional (1-D) active region such that the electron and hole carriers are allowed to move only in that one direction. In contrast, more familiar semiconductor lasers are 3-D double-hetero-structure lasers and 2-D quantum well lasers. Because the 1-D density of states (DOS) becomes more enhanced at the bottom of its band edge than the 2-D or 3-D DOS, a quantum wire laser is expected to show improvement in lasing performance wei ; arakawa ; asada ; arakawa2 .

However, fabrication of a single quantum wire that lases from the lowest quantum state of the wire is difficult. Indeed all quantum wire lasers so far reported are either multiple-wire lasers weg ; rubio ; sirigu ; kim ; aki1 ; aki2 , which have many quantum wires in the active region, or they are wide single-wire lasers which lase only in their excited-subband states kapon ; wata ; pies . The excited-state lasing is typically characterized by a lasing energy above the energy of ground-state spontaneous emission at low pump levels. Since excited subbands allow motion of carriers in directions perpendicular to the axis of the wire, such excited-state wire lasing is not expected to have 1-D characteristics.

In this paper, we demonstrate a single quantum wire laser with only one 1-D subband, and observe stable single-mode lasing from this ground state subband of our quantum wire.

The single quantum wire laser is fabricated by an advanced crystal growth method called cleaved-edge overgrowth with molecular beam epitaxy (MBE) pfei , in which two MBE growth steps are separated by an in situ wafer cleave. In the first MBE growth, we successively grew at 600 oC on a non-doped (001) GaAs substrate, a 500 nm GaAs buffer layer, a 1.5 $`\mu `$m Al0.5Ga0.5As cladding layer, a 250 nm Al0.35Ga0.65As barrier, a 14 nm Al0.07Ga0.93As quantum well (stem well), a 250 nm Al0.35Ga0.65As barrier, 1.5 $`\mu `$m Al0.5Ga0.5As cladding layer, and a 6.5 $`\mu `$m GaAs cap layer. Then, we cleaved the wafer in the MBE chamber, and grew at 490 oC on the newly exposed (110) edge, a 6 nm GaAs quantum well (arm well), a 10 nm Al0.5Ga0.5As barrier, a 111 nm Al0.1Ga0.1As core layer, a 960 nm Al0.5Ga0.5As cladding layer, and a 10 nm GaAs cap layer. After the growth of the arm well, we interrupted growth and annealed the GaAs surface for 10 minutes at 600 oC yoshi .

Figure 1 shows a schematic cross-sectional view of the single quantum wire laser structure. At a T-shaped intersection of a 14 nm-thick stem well and a 6 nm-thick arm well, quantum-mechanical confinement of electrons forms a quantum wire (T-wire). Contour curves in the blow up of the T-wire in Fig. 1 show the wave function of 1-D electrons. Previous experimental and theoretical studies aki1 ; szy  show that such a T-wire has no confined higher 1-D electron subband, and is in the 1-D ground state quantum limit. A photoluminescence measurement shows that energies of the ground-state exciton in the T-wire, arm well and stem well are respectively 1.582, 1.602 and 1.636 eV at 5K. The single quantum wire is embedded in a core of T-shaped optical waveguide (T-waveguide). Contour curves at the top part of Fig. 1 show optical field intensity, or optical mode (transverse mode), for lasing. The optical confinement factor $`\mathrm{\Gamma }`$, defined as the portion of the optical intensity overlapped with the 14 nm $`\times `$ 6 nm region of the quantum wire, is optimized by choosing the thickness of the core layers in the T-waveguide on the basis of a finite-element-method calculation. In the present design, the value of $`\mathrm{\Gamma }`$ is 5$`\times `$10-4.

A series of lasers each with 500 $`\mu `$m long optical cavities were cut from the wafer by (1$`\overline{1}`$0) cleavage, and the cavity-mirror surfaces were coated by 120 nm- and 300 nm-thick gold layers with estimated reflectivities of 97 %. Each laser was pumped optically with cw titanium-sapphire-laser light, mechanically chopped into a 1 % duty ratio to minimize sample heating. The excitation light energy was 1.6455 eV. The pump laser light in a filament shape of about 1 $`\mu `$m width focused by a cylindrical lens and a 0.5 numerical aperture objective lens was incident on the top (110) surface of the sample illuminating the 500 $`\mu `$m long T-wire uniformly. Laser emission was detected via the 120 nm-thick gold mirror which has an estimated transmissivity less than 0.1 %.

Figure 2 shows typical laser emission spectra from a device at 5K for various excitation input powers $`P_{in}`$. At an input power $`P_{in}`$=8.3 mW, multi-mode laser emission of the T-wire is observed at 1.578 eV. It changes to single-mode at $`P_{in}`$=17 mW, and shows a slight red-shift with mode hopping as the input power increases. The energy stability of our quantum wire laser with increasing pumping power is remarkable. The total red-shift from laser threshold to pumping at 260 mW or $`\times `$ 50 threshold is only 1.5 meV. The details of this small shift with increasing pumping power are seen most clearly in Fig. 3. Note that for all pumping powers the lasing energies lie just below the ground-state exciton (1.582 eV) in the T-wire, proving that the lasing indeed occurs in the T-wire ground state.

At higher energies, laser emission from the arm well and the stem well are also observed, as denoted in Fig. 2. Lasing by the arm well, shown with a magnified scale by a factor of 20, occurs only at very high pumping power $`P_{in}`$=260 mW. It cannot occur for lower powers, since photo-excited carriers in the arm well readily empty into the T-wire. At $`P_{in}`$=260 mW, however, the T-wire state saturates allowing the arm well to fill sufficiently to begin separately lasing.

Lasing by the stem well is also observed for input powers above $`P_{in}`$=17 mW. Not surprisingly the emission is stronger than that of the T-wire, and is shown with a reduced scale by a factor of 1/150. In contrast to the case of the arm well, carrier flow from the stem well to the T-wire causes only a minor effect on carrier density in the stem well, since the excitation light generates carriers in the stem well over a large region extending many microns away from the T-wire. The stem well laser emission has multiple longitudinal modes. It shows a larger red-shift of 5 meV at $`P_{in}`$=260 mW. The stem well lasing can be regarded as a prototype of a gain-guided 2-D quantum well laser.

It is interesting to compare the lasing of the T-wire with that of the stem well and the arm well, since it highlights in a single device intriguing differences in performance of a 1-D laser in comparison with 2-D lasers. First, the T-wire tends to lase in a single mode in contrast to multi-mode lasing of the arm and stem wells. This should be partly attributed to a 1-D optical waveguide, and to narrower gain spectra of the T-wire than those of the arm and stem wells. Second, the T-wire shows comparatively small red-shifts of 1.5 meV for increased input power $`P_{in}`$=260 mW, while the stem well shows larger red-shift of 5 meV. This reflects difference in many-body Coulomb interaction effects in a 1-D wire and a 2-D well aki3 . Third, the T-wire lasing has a lower threshold than the others. This last point is all the more remarkable since the T-wire laser has a small optical confinement factor $`\mathrm{\Gamma }`$ = 5$`\times `$10-4 in contrast to large $`\mathrm{\Gamma }`$ = 0.02 in the stem well laser and $`\mathrm{\Gamma }`$ = 0.01 in the arm well laser. Clearly the single T-wire laser has a high enough gain to compensate for the very low $`\mathrm{\Gamma }`$ factor. These three features of the T-wire are also observed in single wire lasers we fabricated with 800 $`\mu `$m-long cavities and also in multiple wire lasers with 20 wires aki2 .

Figure 3 (a) shows lasing spectra of the T-wire more in detail for the various input powers up to $`P_{in}`$=260 mW. Though the width of each emission line is limited by the spectral resolution 0.2 meV of the spectrometer, the results reveal interesting details of longitudinal modes in the T-wire lasing, where modes a \- e denote respective longitudinal lasing modes. Multi-mode lasing with several longitudinal modes centered at a starts at input power of $`P_{in}`$= 5 mW. Transition from multi-mode lasing to single-mode lasing at b occurs at low input power of $`P_{in}`$=10 mW. Single mode lasing at such low input power is one of the remarkable characteristics of the T-wire laser. For increased input power, we observe mode hopping from b to c, to d, and to e with accompanying red-shifts. Below input power $`P_{in}`$=130 mW, lasing-mode separation, or mode hopping distance, is $`\mathrm{\Delta }E_1`$=0.30 meV. This $`\mathrm{\Delta }E_1`$ corresponds to the mode separation of a Fabry-Perot resonator with 500 $`\mu `$m cavity length, when the effective refractive index is 4.1. Above $`P_{in}`$=130 mW, however, it becomes $`\mathrm{\Delta }E_2`$=0.60 meV. The reason of the doubled mode separation at high pumping is not understood.

Figure 3 (b) plots laser emission intensity of the T-wire at 5K as a function of input power. The inset shows magnified plots near threshold, which gives the value of the threshold input power at 5 mW.

For higher input powers above $`P_{in}`$=20 mW, laser emission intensities of the T-wire plotted by black dots and denoted as Total in Fig. 3 (b) increases showing irregular oscillations and dips. We also plot intensities of respective longitudinal modes a \- e in T-wire lasing. It turns out that the dips in the total intensity occur exactly at those input powers where the laser switches from one lasing mode to the next. The reason for this interesting effect is also not at present understood.

Figure 4 (a) shows the laser emission spectra of the T-wire at 5, 20, 40 and 60K for input powers of 42 mW and 130 mW. Single-mode lasing is observed up to 40K, while only multi-mode lasing occurs up to 60K. With the increase in temperature, the spectral peaks are shifted toward lower energies as expected due to the shrinkage in the GaAs band-gap energy.

Figure 4 (b) plots laser emission intensity of T-wire at 5, 20, 40 and 60 K against input powers up to 100 mW. Notice that the lasing intensities at 20 K are larger than those at 5K. This is caused by enhanced carrier diffusion in the stem well and the arm well which provides more carriers to the T-wire at 20K than at 5K. At 40 K and 60 K, however, lasing intensities become smaller and the lasing threshold increases. This is probably because carriers in the T-wire are thermally activated to escape to distant regions of the arm well and the stem well.

In summary, we fabricated a single quantum wire laser in the 1-D quantum limit. We observed lasing at the ground state of a 14 nm $`\times `$ 6 nm wire at 5 - 60 K via optical pumping. It shows a low threshold input power of 5 mW at 5 K.

One of the authors (H. A.) acknowledges the financial support from the Ministry of Education, Culture, Sports, Science and Technology, Japan.


text:
# Spherical collapse with heat flow and without horizon

## Abstract

We present a class of solutions for a heat conducting fluid sphere, which radiates energy during collapse without the appearance of horizon at the boundary at any stage of the collapse. A simple model shows that there is no accumulation of energy due to collapse since it radiates out at the same rate as it is being generated.

One of the interesting problems related to gravitational collapse of stars is to study the solutions of various spherically symmetric fluid distributions with radial heat flux in the interior satisfying all the energy conditions and having reasonable physical behaviour. In this case the exterior spacetime is described by the radiating Vaidya metric . The junction conditions at the boundary of the isotropic fluid sphere with dissipation in the form of the radial heat flux matching with the radiating Vaidya metric were previously studied by Santos . In particular special attention was given to the model which from an initial static configuration starts collapsing when dissipation takes place in the form of heat flux . The initial static configuration evolves until the horizon is formed.

In this letter we present an interesting solution of a collapsing fluid sphere with radial heat flux. The interesting feature of it is that collapse never encounters horizon and it continues until a (naked) singularity is reached. The matter content of this model has reasonable properties and the non-occurrence of the horizon is due to the fact that the rate of mass loss is balanced by the fall of boundary radius. In fact, one gets a class of such solutions although the explicit ex

ample is given for a very simple heat conducting fluid sphere, which satisfies all the energy conditions and has very reasonable physical properties.

The metric in the interior of a shearfree spherically symmetric fluid distribution is given by

$$ds^2=A^2(r,t)dt^2+B^2(r,t)[dr^2+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2]$$
(1)
The stress energy tensor of a nonviscous heat conducting fluid reads

$$T^{\mu \nu }=(\rho +p)v^\mu v^\nu +pg^{\mu \nu }+q^\mu v^\nu +q^\nu v^\mu .$$
(2)
The heat flow vector $`q^\mu `$ is orthogonal to the velocity vector so that $`q^\mu v_\mu =0`$. Assuming comoving coordinates $`v^\mu =(g_{00})^{1/2}`$ and $`q^\mu =(0,q^1,0,0)`$ nontrivial Einstein equations are

$`\rho `$ $`=`$ $`{\displaystyle \frac{4B^{}}{rB^3}}+{\displaystyle \frac{3\dot{B}^2}{A^2B^2}}+{\displaystyle \frac{B^2}{B^4}}{\displaystyle \frac{2B^{\prime \prime }}{B^3}}`$ (3)
$`p`$ $`=`$ $`{\displaystyle \frac{2A^{}}{rAB^2}}+{\displaystyle \frac{2B^{}}{rB^3}}+{\displaystyle \frac{2\dot{A}\dot{B}}{A^3B}}{\displaystyle \frac{\dot{B}^2}{A^2B^2}}{\displaystyle \frac{2\ddot{B}}{BA^2}}+{\displaystyle \frac{2A^{}B^{}}{AB^3}}+{\displaystyle \frac{B_{}^{}{}_{}{}^{2}}{B^4}}`$ (4)
$`p`$ $`=`$ $`{\displaystyle \frac{A^{}}{rAB^2}}+{\displaystyle \frac{B^{}}{rB^3}}+{\displaystyle \frac{2\dot{A}\dot{B}}{A^3B}}{\displaystyle \frac{\dot{B}^2}{A^2B^2}}{\displaystyle \frac{2\ddot{B}}{BA^2}}+{\displaystyle \frac{A^{\prime \prime }}{AB}}{\displaystyle \frac{B_{}^{}{}_{}{}^{2}}{B^4}}+{\displaystyle \frac{B^{\prime \prime }}{B^3}}`$ (5)
$`q^1`$ $`=`$ $`{\displaystyle \frac{2A^{}\dot{B}}{A^2B^3}}{\displaystyle \frac{2B^{}\dot{B}}{AB^4}}+{\displaystyle \frac{2\dot{B}^{}}{AB^3}}`$ (6)
The exterior Vaidya metric already mentioned is given explicitly in the form

$$ds^2=\left(1\frac{2M(v)}{\overline{r}}\right)dv^22d\overline{r}dv+\overline{r}^2(d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2)$$
(7)
where v is the retarded time and M(v) is the exterior Vaidya mass. The junction conditions then yield the following system of equations valid on the boundary $`r=r_\mathrm{\Sigma }`$ ,

$`(rB)_\mathrm{\Sigma }`$ $`=`$ $`\overline{r}_\mathrm{\Sigma }`$ (8)
$`p_\mathrm{\Sigma }`$ $`=`$ $`(q^1B)_\mathrm{\Sigma }`$ (9)
and

$$m_\mathrm{\Sigma }=\left[\frac{r^3B\dot{B}^2}{2A^2}r^2B^{}\frac{r^3B^2}{2B}\right]_\mathrm{\Sigma }$$
(10)
where $`m_\mathrm{\Sigma }`$ is the mass function calculated in the interior at $`r=r_\mathrm{\Sigma }`$ . The junction conditions yield $`m_\mathrm{\Sigma }=M(v)`$.

At this stage we consider a particular form of the metric coefficients given in (1) and choose them separable in r and t co-ordinates as

$`A`$ $`=`$ $`a(r)`$ (11)
$`B`$ $`=`$ $`b(r)R(t)`$ (12)
so that $`\rho `$ , p and $`q^1`$ can now be expressed in the form

$`\rho `$ $`=`$ $`{\displaystyle \frac{1}{R^2}}\left[{\displaystyle \frac{3}{a^2}}\dot{R}^2{\displaystyle \frac{1}{b^2}}\left({\displaystyle \frac{2b^{\prime \prime }}{b}}{\displaystyle \frac{b^2}{b^2}}+{\displaystyle \frac{4b^{}}{rb}}\right)\right]`$ (13)
$`p`$ $`=`$ $`{\displaystyle \frac{1}{R^2}}\left[{\displaystyle \frac{1}{a^2}}(2R\ddot{R}+\dot{R}^2)+{\displaystyle \frac{1}{b^2}}\left({\displaystyle \frac{b^2}{b^2}}+{\displaystyle \frac{2a^{}b^{}}{ab}}+{\displaystyle \frac{2}{r}}\left({\displaystyle \frac{a^{}}{a}}+{\displaystyle \frac{b^{}}{b}}\right)\right)\right]`$ (14)
$`q^1`$ $`=`$ $`{\displaystyle \frac{2a^{}\dot{R}}{R^3a^2b^2}}`$ (15)
where ‘dot’ and ‘dash’ indicate derivatives with respect to time and radial coordinate. The isotropy of pressure would give the equation,

$$\frac{a^{\prime \prime }}{a}+\frac{b^{\prime \prime }}{b}2\frac{b_{}^{}{}_{}{}^{2}}{b^2}2\frac{a^{}b^{}}{ab}\frac{a^{}}{ra}\frac{b^{}}{rb}=0.$$
(16)
The boundary condition (9) now yields at $`r=r_\mathrm{\Sigma }=r_o`$ in view of (14) and (15)

$$2R\ddot{R}+\dot{R}^2+m\dot{R}=n$$
(17)
where m and n are constants. Such an equation was first given in  in the context of a very particular case of Maiti’s solution . The general solution of (16) in closed form is not available and a very simple solution is

$$R(t)=Ct.$$
(18)
An interesting feature of the above solution is that at the boundary, $`m_\mathrm{\Sigma }/\overline{r}\mathrm{\Sigma }`$ turns out to be independent of time. It can be demonstrated very easily from the expression (10). We then get using the solution (18)

$$\frac{2m_\mathrm{\Sigma }}{\overline{r}_\mathrm{\Sigma }}=\frac{2m_\mathrm{\Sigma }}{(rB)_\mathrm{\Sigma }}=2\left[\frac{C^2r_0^2b_0^2}{2a_0^2}\frac{r_0b_0^{}}{b_0}\frac{r_0^2b_0^2}{2b_0^2}\right]$$
(19)
where $`r_0`$ is written for the comoving radial coordinate at the boundary and $`b(r_0)=b_0,b^{}(r=r_0)=b_0^{}`$. It is quite possible to adjust the aparameters in the above equation so as to keep $`2m_\mathrm{\Sigma }/\overline{r}_\mathrm{\Sigma }<1`$ in order to avoid the appearance of horizon at the boundary.

A simple collapsing model with heat flow but without horizon

We set $`b(r)=1`$ in the pressure isotropy equation (16) and consider the special solution ,

$$A=a(r)=(1+\xi _0r^2).$$
(20)
and then in view of equations (17-20) we get

$$m=4\xi _0r_0,n=4\xi _0(1+\xi _0r_0^2),C=\frac{1}{2}\left(|m|+(m^2+4n)^{1/2}\right).$$

In the solution $`R(t)=Ct`$, the constant C is chosen to be positive so that collapsing phase corresponds to $`\mathrm{}<t<0`$. We further have

$$C^2<4\xi _0(1+\xi _0r_0^2).$$
(21)
The explicit expressions for the density, pressure and heat flow vector are now

$`\rho `$ $`=`$ $`{\displaystyle \frac{3}{t^2(1+\xi _0r^2)^2}}`$ (22)
$`p`$ $`=`$ $`{\displaystyle \frac{1}{t^2(1+\xi _0r^2)^2}}\left[{\displaystyle \frac{4\xi _0}{C^2}}(1+\xi _0r^2)1\right]`$ (23)
$`q^1`$ $`=`$ $`{\displaystyle \frac{4\xi _0r}{(1+\xi _0r^2)^2}}{\displaystyle \frac{1}{C^2t^3}}`$ (24)
All the above quantities diverge at $`t0`$. The above expressions show that $`\rho >0,p>0`$ and $`\rho ^{}<0`$, and $`p^{}<0`$ would require $`C^2<2\xi _0(1+\xi _0r^2)`$. It would be satisfied throughout the interior if $`C^2<2\xi _0`$. Further $`(\rho p)>0`$ everywhere within the interior implies $`C^2>\xi _0(1+\xi _0r_0^2)`$. All the above physically reasonable conditions will be satisfied provided the following restrictions are imposed on the magnitude of $`C^2`$

$$2\xi _0>C^2>\xi _0(1+\xi _0r_0^2).$$
(25)
Since there is heat flow in the radial direction the fluid must satisfy another condition in order to be consistent with all the energy conditions and it is $`(\rho +p)>2|q|`$, where $`|q|(g_{\mu \nu }q^\mu q^\nu )^{1/2}`$. This would require,

$$1+\frac{2\xi _0}{C^2}(1+\xi _0r^2)>4\xi _0r/C$$

which would always be true as it could be written as

$$\left[1\frac{2\xi _0r}{C}\right]^2>\frac{2\xi _0}{C^2}(1\xi _0r^2).$$
(26)
Now in view of eq. (19), since b = 1 it follows immediately that

$$1\frac{2m_\mathrm{\Sigma }}{\overline{r}_\mathrm{\Sigma }}=\left[1\frac{C^2r_0^2}{(1+\xi _0r_0^2)^2}\right]$$
(27)
So once we have $`C^2<1/r_0^2+\xi _0^2r_0^2+2\xi _0`$, it is clear that the boundary surface can never reach the horizon. Here in fact both the mass function and the area radius of the radiating sphere vary linearly with t and hence the ratio is independent of time. It may be further noted that though the density, pressure and curvatures would diverge as $`1/t^2`$ when $`t0`$, yet the mass function which would go as t , would remain finite and go to zero in this limit. Thus naked singularity would be a weak curvature singalarity. For strong curvature singularity the curvature should diverge at least as $`t^3`$. This simple example shows that there is no accumulation of energy due to collapse as it is being radiated out at the same rate as it is being generated.

Acknowledgement: AB would like to thank IUCAA for hospitality.