text:
# Partially ionizing the universe by decaying particles

## I Introduction

It was known for many years that the universe, which once became almost neutral after recombination, reionized before $`z5`$ by using Gunn-Peterson (GP) test of high redshift quasars. It was not until recently, however, that direct evidences to determine the reionization epoch were found.

First evidence of the reionization epoch was obtained by Sloan Digital Sky Survey (SDSS) SDSS . The spectra of quasars observed by SDSS team revealed that there existed neutral hydrogen at the redshift $`z6`$ by GP test, which implies that we have started to see the completion of reionization, although the amount of neutral hydrogen does not have to be very large to explain the observed GP trough.

Recent Wilkinson Microwave Anisotropy Probe (WMAP) observation provided second important evidence of reionization WMAP . WMAP found that the optical depth due to the reionization is $`\tau _{op}=0.17\pm 0.04`$, which may correspond to the reionization epoch $`z20`$.

Therefore, it is observationally required that the reionization process was very efficient at around $`z20`$, but the process did not complete until $`z6`$. The conventional view of the reionization process is that the universe is reionized by UV fluxes from early formed stars and quasars. Many calculations have shown that the universe can be reionized at $`z610`$ Loeb . Thus, it is possible to explain the full reionization at $`z=6`$ observed by SDSS but not what WMAP observed. Recent studies Ciardi ; FK03  showed that the reionization epoch can be earlier assuming the non-standard initial mass function (IMF) or unrealistically large UV fluxes from galaxies. However, the UV fluxes from stars and quasars increase as time and soon exceed the critical value to fully reionize the universe. Therefore it is expected that the universe achieved the full reionization much earlier than $`z6`$ in this scenario. Therefore we need to consider really extreme and perhaps unrealistic assumptions about IMF, escaping UV fraction from galaxies and so on, to have the reionization process with two stages; for example, see Astro .

Alternative to along this line, one may also consider other sources for ionizing hydrogens and heliums. It could be decaying particles, but it has been considered that enduring photon emission from particle decay disables the picture of GP trough confirmed by SDSS. In this article, however, we will show that the optical depth $`0.17`$ can be explained naturally by decaying particles with rather large parameter space of mass and lifetime; in some cases, the lifetime could be longer than the age of the universe. We follow the evolution of hydrogen, helium, and electron temperature including UV sources of quasars, stars, and decaying particles. Here we simply assume the particles decay into two photons with energy half of the particle mass. Notice that Hansen and Haiman HH  tried to explain the optical depth of $`\tau _{op}0.17`$ by a particle of mass $`200`$ MeV and the lifetime $`4\times 10^{15}`$ sec, which is $`z20`$ reported by WMAP team WMAP . They considered the decaying particle as some sterile neutrino which decay into electrons (and positrons) in order to achieve half partially reionized era before $`z6`$. Here we will see that photons emitted through particle decay keeps much smaller partially ionization before $`z6`$, and explains rather large optical depth.

## II UV sources for ionization

We assume that the particle $`\varphi `$ emits two photons with monochromatic energy of half mass of that particle, i.e., $`E_\gamma =m_\varphi /2`$. The number density of $`\varphi `$-particle is written as

$$n_\varphi =n_\varphi (0)(1+z)^3e^{\frac{t}{\tau _\varphi }},$$
(1)
where $`\tau _\varphi `$ is the lifetime of $`\varphi `$-particle. If $`E_\gamma >13.6`$ eV, those photons emitted can ionize hydrogen atoms. Then one can write a source term for the decaying particle as

$$\left(\frac{dn_\gamma }{dt}\right)_{dp}=\frac{n_\varphi (t)}{\tau _\varphi }.$$
(2)
In order to calculate how many photons are emitted, the abundance, mass, and the lifetime of the particle should be fixed. We parameterize them as free parameters by $`\mathrm{\Omega }_\varphi `$, $`E_\gamma `$, and $`\tau _\varphi `$. We will see below the allowed region in these parameter space. Notice that we restrict ourselves only for $`E_\gamma 10^5`$ eV. This is because $`E_\gamma `$ MeV case involves e+e- pair production, and it is too complicated. Moreover, $`E_\gamma O(10)`$ MeV may destroy light elements such as D, T, 3He, and 4He, spoiling successful big bang nucleosynthesis.

Other than decaying particles, we also include usually considered contributions to reionizing photons from quasars and stars. UV photons emitted from early formed stars are indeed responsible for achieving full ionization at $`z6`$. UV photons from quasars may keep the ionization fraction unity, but this contribution is very small, and may not be needed for that purpose (Of course, quasars exist and there are contributions from them). We follow the argument of Fukugita and Kawasaki FK94 ; FK03  to follow the thermal history from $`z>10^3`$ including reionization epoch, calculating the ionization fraction of hydrogen and helium, the photon spectrum, and the electron temperature (see also FK90 ; FK93  for treatment of high energy photons from decaying particle).

Some of the difference from Ref. FK94 ; FK03  are adoption of decreasing star formation rate for smaller collapsed objects. In Refs.FK94 ; FK03  it was assumed that the collapsed objects which satisfy the cooling condition become star-forming galaxies and a constant fraction $`f`$ of baryons go into stars. In the present paper, we introduce the star-forming efficiency $`f_{star}`$ which depends on the mass of the objects as

$$f_{star}=\{\begin{array}{cc}f\left(\frac{M}{3\times 10^{12}M_{}}\right)^{0.8}& (M<3\times 10^{12}M_{}),\hfill \\ f& (M3\times 10^{12}M_{}),\hfill \end{array}$$
(3)
where $`M`$ is the mass of the collapsed object. With $`f_{star}`$ the mass function of galaxies obeys the form of an empirical Schecheter function at $`z=0`$. Since we expect that the UV fluxes from stars reionize the universe at relatively low redshift ($`z6`$), we do not have to take extreme values for model parameters such as the UV escape fraction $`F_{UV}`$ and $`f`$. We take $`F_{UV}=0.04`$ We take $`F_{UV}=0.0050.01`$ for longer lifetime than the age of the universe for $`E_\gamma =15`$ eV. and $`f=1`$ which give the metalicity $`Z0.01`$.

In order to estimate the UV contribution from quasars we use the luminosity function in Refs.Boyle  and Fan , from which we obtain the UV emissivity ( 1/s/cm3) at 4400Åas

$$\left(\frac{dn_\gamma }{dt}\right)_{QSO}=\{\begin{array}{cc}5\times 10^{23}10^{0.47(z2)}h\hfill & (z>2),\hfill \\ 5\times 10^{23}h\left(\frac{1+z}{3}\right)^{3.5}\hfill & (z2),\hfill \end{array}$$
(4)
where $`h`$ is the Hubble constant in units of 100km/s/Mpc. We also assume a power low spectrum for quasars,

$$\left(\frac{d^2n_\gamma }{dtdE_\gamma }\right)_{QSO}\{\begin{array}{cc}E_\gamma ^{0.7}\hfill & (E_\gamma <10.2\mathrm{eV}),\hfill \\ E_\gamma ^{1.4}\hfill & (E_\gamma 10.2\mathrm{eV}).\hfill \end{array}$$
(5)
We set the cosmological parameters following WMAP-only WMAP : $`H_0=72`$ km/s/Mpc, $`\mathrm{\Omega }_m=0.29`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=1\mathrm{\Omega }_m`$, $`\mathrm{\Omega }_b=0.047`$, and $`\sigma _8=0.9`$. Apparently, we are dealing with $`\mathrm{\Lambda }`$CDM universe, and the decaying particle is assume to be a part of cold dark matter, Particles with very small mass such as $`10^5`$ eV can be regarded as cold dark matter in such a way that a scalar condensate is corehently oscillating. but could be hot dark matter with negligible contribution to the total density of the universe.

## III Ionization history

We follow the evolution of the fraction of ionized hydrogen and helium, in addition to electron temperature. Figure 1 shows ionization history of hydrogen for $`E_\gamma =15`$ and $`10^5`$ eV for various lifetime $`\tau _\varphi `$. In each case, the abundance of $`\varphi `$-particles is determined in such a way that the optical depth reads $`\tau _{op}0.17`$, which we define by

$$\tau _{op}=_0^{\mathrm{}}𝑑z\sigma _T\left(\frac{dt}{dz}\right)[n_{HII}n_{HII}|_{sr}],$$
(6)
where $`\sigma _T`$ is the Thomson cross section and $`n_{HII}|_{sr}`$ is the number density for standard recombination. We subtract this term in order to estimate only the effect of reionization.

In cases for $`E_\gamma 30`$ eV, the lifetime of $`\varphi `$-particle is allowed for many orders of magnitude. For longer lifetime than the age of the universe, it is the value of $`n_\varphi /\tau _\varphi `$ that concerns with the amount of photon emission, so ionization histories of these cases are identical to that of $`\tau =10^{18}`$ sec shown in Fig.1. Notice that the upper limit of the lifetime is around $`4\times 10^{22}`$ sec, when $`\varphi `$-particles is (cold) dark matter, i.e., $`\mathrm{\Omega }_\varphi \mathrm{\Omega }_m`$. For shorter lifetime, the abundance should be much less than that contributes to the energy density of the universe. See Fig. 2.

On the other hand, if larger energy of photons are emitted, we have to be careful about cosmic X-ray backgrounds. We adopt the following simple formulas to compare with the calculated photon flux KY :

$`F_{\gamma ,obs}`$ $``$ $`8\left({\displaystyle \frac{E_\gamma }{\mathrm{keV}}}\right)^{0.4}(0.2\mathrm{keV}E_\gamma 25\mathrm{k}\mathrm{e}\mathrm{V}),`$ (7)
$``$ $`380\left({\displaystyle \frac{E_\gamma }{\mathrm{keV}}}\right)^{1.6}(25\mathrm{k}\mathrm{e}\mathrm{V}E_\gamma 350\mathrm{k}\mathrm{e}\mathrm{V}),`$ (8)
$``$ $`2\left({\displaystyle \frac{E_\gamma }{\mathrm{keV}}}\right)^{0.7}(350\mathrm{k}\mathrm{e}\mathrm{V}E_\gamma 2\mathrm{M}\mathrm{e}\mathrm{V}),`$ (9)
where $`F_{\gamma ,obs}`$ is measured in units of (cm2 sr sec)-1.

We find that long lifetime is not allowed, because the emitted photon flux exceeds the observed diffuse photon background. Thus, the lifetime is rather restricted in the narrow range, and should typically be $`\tau _\varphi 10^{16}`$ sec. This constraint is taken into account in Fig. 2, thus the parameter space is restricted for $`E_\gamma 10^2`$ eV. Hydrogen ionizing histories for $`E_\gamma =10^5`$ eV are shown in Fig. 1 as an example. They look very similar to those for $`E_\gamma =1530`$ eV. In this case, even $`\tau _\varphi =10^{16}`$ sec seems to be excluded as is shown in Fig. 2.

Evolutions of ionization fraction of hydrogen look more or less similar in any emitted photon energy for the same lifetime. In order to distinguish them, we need to look at the ionization history of helium. Since larger photon energy can excite helium atoms, one can see the qualitative difference from low energy emitted photons. In Figs. 3 and 4, ionization histories of HeII and HeIII are shown, respectively. We only plotted for $`E_\gamma =15`$ and $`10^5`$ eV cases as in Fig. 1.

It is trivial to say that rather high ionization fraction begins at high redshift in large $`E_\gamma `$ cases. Therefore, although it will be very difficult to observe ionization fraction of helium at $`z20`$, it might be one of the tools for (dis)proving the particle decay scenario.

## IV Cosmic microwave background spectrum

The evolutions of ionization fraction of hydrogen and helium look very different from those instantaneous reionization history usually considered in the literatures, although we obtain the optical depth of $`0.17`$. So it is interesting to know what we can see, or more precisely, whether we can distinguish them, in the cosmic microwave background (CMB) power spectrum. Figure 5 shows the cross correlation spectra between temperature and E-mode polarization. Our models seem to be less power in the small $`\mathrm{}`$ region than the model with instantaneous reionization. This is because particle decay takes place for longer period, so the power spreads to larger $`\mathrm{}`$ for fixed value of the optical depth. From the TE spectrum, it is not apparent to distinguish between these models within these error bars. Even difference between the instantaneous reionization and ionization by decaying particles can hardly seen. Moreover, notice that larger optical depth $`\tau _{op}0.43`$ may account for the observational results of the TE spectrum, too. Of course, it is known that larger optical depth makes the first acoustic peak of the TT spectrum lower. However it can be adjusted by introducing larger spectral index, so-called $`\tau _{op}n`$ degeneracy.

Unlike TE spectra, however, EE spectra seem to serve us a much powerful tool for discriminating among different ionization histories,§§§ Even in more conventional cases, differences of EE spectra seem to have ability to distinguish among different ionization histories Holder ; HuHolder . as shown in Fig. 6. Peculiar feature of enduring photon emission from particle decays is larger power around $`\mathrm{}=`$ (a few)$`\times 10`$. In addition, there is different feature between long and short lifetime: for long lifetime, the spectrum is rather flat, while there is a dip at $`\mathrm{}10`$ for short lifetime. Thus, EE spectrum hoped to be seen by Planck satellite or even WMAP should (dis)prove the photons from decaying particles in the near future.

To end this section, we comment on the Compton $`y_c`$-parameter. We also follow its evolution, and obtained that the final value is (several)$`\times 10^7`$ in all the parameters we choose, which is well below the current observational limit COBE .

## V Conclusion and discussion

We have studied comprehensively the ionization history of the universe by photons emitted from decaying particles, together with the effects of stars, quasars, and shown that the decaying particle can explain rather large optical depth observed by WMAP without conflicting other astrophysical and cosmological constraints. One of the amazing fact is that the lifetime of order $`10^{15}`$ sec (in other words, $`z20`$) is not necessarily required for successful scenario, and larger parameter space is allowed. Since the early star formation ejecting UV photons may be difficult to lead large optical depth without adjusting parameters beyond their natural values, this alternative is at least as good as them.

For $`E_\gamma =10^210^5`$ eV, larger lifetime is prohibited, because the emitted flux exceeds the observed diffuse photon background, and $`\tau _\varphi 10^{14}10^{16}`$ sec is only allowed in these cases. On the other hand, if some bump is observed in the photon spectrum, it could be a great clue for decaying particle scenario.

Longer lifetime is also allowed for $`E_\gamma 30`$ eV. In the limiting case, the decaying particle can be the dark matter, provided that its lifetime is around $`4\times 10^{22}`$ sec. Notice that our situation is different from Ref. Sciama , in which they investigated for decaying neutrino dark matter and thus using hot dark matter scenario. In this case, ionization fraction becomes unity at $`z2030`$, which conflicts with SDSS results. In addition, EURD observation may falsify decaying particle scenario with $`E_\gamma 15`$ eV, but it needs careful look, since their conclusion will change if the abundance of galactic dark matter near the Sun is one third of their applied value EURD .

Observation of CMB spectrum is very powerful. Although instantaneous reionization suggests optical depth $`\tau _{op}0.17`$, higher optical depth may be favored for the observed TE spectrum by WMAP if the decaying particle scenario is correct, where slow reionization takes place. It is somewhat difficult to discriminate from each other only by TE spectrum. However, the differences in the EE spectrum is tremendous so we hope that even the lifetime of the decaying particle will be measured by Planck.

Finally, we comment on what can be a candidate for this decaying particle. In the context of no-scale supergravity, scalar particles remain massless at the tree level. They may acquire their masses at one-loop order. Among them, saxion, the scalar partner of axion in supersymmetric theory, might be the best candidate. If we consider the hadronic axion model, the lifetime is estimated as GY

$$\tau _s10^{15}\left(\frac{F_{PQ}}{10^{15}\mathrm{GeV}}\right)^2\left(\frac{m_s}{100\mathrm{e}\mathrm{V}}\right)^1\mathrm{sec},$$
(10)
where $`F_{PQ}`$ is the axion decay constant. It perfectly fits into the allowed region shown in Fig. 2.

## Acknowledgments

NS is supported by Japanese Grant-in-Aid for Science Research Fund of the Ministry of Education, No.14340290.


format: text

text:
# Multifractal Analysis for the Dynamical Heterogeneity in Strongly Correlated Many-Body Systems

## 1 Introduction

We have numerically studied the dynamical heterogeneity in strongly correlated many-body systems, for example, the critical spin stateSN1 , the supercooled liquid near the glass transitionSN2  and the turbulenceSN3 . These systems are scale-invariant and multifractal. In this paper we briefly review how we quantify the dynamical heterogeneity and discuss in detail the turbulence simulation data recently obtained.

In order to quantify the dynamical heterogeneity we have two measures. One is the non-equilibrium parameter $`q`$ of the Rényi-Tsallis distribution function. The other is the singularity spectrum $`f(\alpha )`$. These two measures are closely related and important ingredients of the multifractal analysis.

In the numerical study of the critical spin stateSN1  we have found that the $`q`$-parameter represents the degree of non-equilibrium. We have actually quantified the deviation of the q-parameter from the equilibrium value according to the spatio-temporal scale of the observation.

In the numerical study of the the supercooled liquid near the glass transitionSN2  we have found that the broken-bond distribution which reflects the nature of the cooperatively rearranging region is well described by the singularity spectrum $`f(\alpha )`$ of multifractal. The width of the spectrum becomes broader as approaching the glass transition. Such a broadening is similar to that observed in the numerical study of the Anderson-localization transition.

## 2 Turbulence Simulation

In the numerical study of the turbulenceSN3  we have adopted the lattice Boltzmann method. Our numerical turbulence on $`200^3`$ cubic lattice points is sustained by a random forcing and the Reynolds number is about 500.

### 2.1 mean-field analysis

Turbulence is one of the typical phenomena with multiscale motions. Each phenomena of a scale strongly couples with all the other scales of turbulent motion. In order to analyze such a system a scale-dependent entropy, the so-called $`ϵ`$-entropy $`h(ϵ)`$, works well. For example, the time series of the velocity field in turbulence leads to a non-trivial scaling relation, $`h(ϵ)ϵ^3`$, expected from Kolmogorov’s scaling. Here $`ϵ`$ is the scale of the observation. The existing experimental data are consistent with this scaling. By our numerical experiment simulating the Navier-Stokes equation we have shown the consistency of the scaling. In contrast to the evaluation of the energy spectrum, which is usually employed for testing Kolmogorov’s scaling and determined by the two-point correlation function, Kolmogorov’s scaling is easily observed in the $`ϵ`$-entropy, since it is determined by the mean of the exit time from the observation window and fluctuations are averaged out. Such a unifractal scaling is a mean-field description and fluctuations can be taken into account in a multifractal analysis as shown in the following. In the above mentioned mean-field description the velocity difference in temporal and spatial directions have the same fractal scaling exponent. Thus we have confirmed Kolmogorov’s scaling and Taylor’s hypothesis at the same time.

### 2.2 fluctuation analysis

In order to analyze the fluctuations we use two measures, the non-equilibrium parameter $`q`$ of the Rényi-Tsallis distribution function and the singularity spectrum $`f(\alpha )`$.

The probability distribution function (PDF) for the velocity or vorticity field is non-Gaussian. The non-Gaussianity is quantified by the $`q`$-parameter. The strong non-Gaussianity is observed at small spatio-temporal scale comparable to that for the coherent vortex in the PDF for the velocity or vorticity difference between two space-time points. By using the wavelet denoising we have clarified that the non-Gaussianity or dynamical heterogeneity results from the existence of the coherent vortex which is a strongly-correlated non-equilibrium region.

The strength of the vortices is shown in Fig. 1 to be intermittent. The spatial distribution of the coherent vortices is quantified by the singularity spectrum of multifractal $`f(\alpha )`$ as shown in Fig. 2.

The remainder of this subsection is devoted to some speculations. As seen in our previous studySN2  the width of the singularity spectrum depends on the degree of intermittency so that we expect some Reynolds-number dependence of the spectrum. It has been discussed by many authors and should be claryfied in future systematic study. In our present study the density of the vortex is dilute. Thus the singularity spectrum in Fig. 2 describes the spatial distribution of relatively free vortices. As the Reynolds number is increased the density increases. In this case the interaction among vortices becomes important and the correlation length of the vorticity fluctuation becomes large. While we can observe the vortex only as an individual elementary excitation in our numerical experiment, some collective excitation is expected to dominate at higher Reynolds number. In the limit of divergently large Reynolds number the correlation length becomes divergently large so that we can expect full scale-invariance. We can find a resemblance to the case of the scaling theory in polymers where an ideal scaling relation is realized for dense solutions where polymers are strongly entangled. In the limit of high Reynolds number each boxes counting the coherent vortex for calculating $`f(\alpha )`$ will be filled by almost equal number of vortices so that intermittency will disappear and unifractal Kolmogorov’s scaling will prevail.

### 2.3 summary for turbulence simulation

Numerical simulation data, in real space and time, for a forced turbulence on the basis of the lattice Boltzmann method have been analyzed by unifractal and multifractal schemes.

Our new findings are summarized into two points. First in the unifractal analysis using the exit-time statistics we have verified Kolmogorov’s scaling and Taylor’s hypothesis at the same time. Second in the analysis using the Rényi-Tsallis PDF and the wavelet denoising we have clarified that the coherent vortices sustain the power-law velocity correlation in the non-equilibrium state.

Finally in the multifractal analysis it is clarified that the intermittent distribution of the coherent vortices in space-time is described as a multifractal.


text:
# Simultaneous visible and infrared spectro-polarimetry of a solar internetwork region

## 1 Introduction

Most of the solar surface appears non-magnetic when it is observed in routine synoptic magnetograms (e.g., those obtained at Kitt Peak111http://www.noao.edu/kpno/ or with the MDI instrument on board of SOHO222http://sohowww.nascom.nasa.gov/). However, magnetic fields are detected almost everywhere, also in internetwork (IN) regions, when the polarimetric sensitivity and the angular resolution exceed a threshold. These magnetic fields of the quiet Sun, elusive but now accessible to many existing spectro-polarimeters, have received an increasing interest during the last years. They produce only weak signals, but they cover much of the solar surface and, therefore, they may be carrying most of the unsigned magnetic flux and energy existing on the solar surface at any given time (e.g., Sánchez Almeida 2003, and references therein). This makes the IN potentially important to understand the global magnetic properties of the Sun.

The observational study of the IN fields began in the seventies (Livingston & Harvey 1975; Smithson 1975), but it is still in an initial phase due to the complexity of the magnetic topology. Opposite polarities seem to coexist in one and the same resolution element of typically 1″(see Sánchez Almeida et al. 1996; Sánchez Almeida & Lites 2000; Sigwarth et al. 1999; Lites 2002; Khomenko et al. 2003). In addition, large systematic differences of intrinsic field strength are deduced depending on the diagnostic technique. Techniques based on visible lines favor kG field strengths (Grossmann-Doerth et al. 1996; Sigwarth et al. 1999; Sánchez Almeida & Lites 2000; Socas-Navarro & Sánchez Almeida 2002; Domínguez Cerdeña et al. 2003a, b), whereas inferences based on infrared (IR) lines indicate sub-kG field strengths (Lin & Rimmele 1999; Khomenko et al. 2003). (We add for completeness that estimates based on the Hanle effect indicate the presence of even weaker fields, e.g., Faurobert-Scholl et al. 1995; Bianda et al. 1999.) Sánchez Almeida & Lites (2000) conjectured that this systematic difference actually reflects the existence of a continuous distribution of field strengths ranging from sub-kG to kG (see also Socas-Navarro & Sánchez Almeida 2002). Different observations are biased towards a particular part of such a distribution, producing the observed systematic difference of field strengths. Even more, Sánchez Almeida & Lites (2000) put forward a specific physical mechanism that biases the IR observations towards sub-kG fields. It is caused by the vertical gradient of field strength existing in any magnetic structure (provided they satisfy some sort of mechanical balance with the mean solar photosphere). The Zeeman splitting of the IR lines is so large that the absorption in various layers, with different field strengths, occurs at substantially different wavelengths. The absorption is thus spread over a wide range of wavelengths, which dilutes the polarization signals. On the other hand, sub-kG fields do not disperse the IR absorption. In this case the IR lines are in the so-called weak field regime so, independently of the field strength, the absorption is always produced at the same wavelengths, and a significant polarization signal builds up. Socas-Navarro & Sánchez Almeida (2003) show that the existence of a horizontal gradient can also explain the observed discrepancy. In this particular case a distribution of field strengths is chosen such that the visible lines reflect kG fields because most of the magnetic flux is in the form of kG fields. In the IR lines, however, the kG signal is much reduced because of the wavelength smearing effect described above.

So far the seemingly contradictory visible and IR observations were taken independently (by different authors at different times). If the conjecture on the coexistence of weak and strong fields were correct, a simultaneous observation of the visible and IR lines would have to give co-spatial signals showing both, signatures of kG (visible) and of sub-kG (IR) fields. The present study was meant to test this specific prediction.

The simultaneous observations and few details of the data reduction are described in § 2. Then we explain the data analysis, the classification of the polarimetric profiles, and the determination of magnetic field strengths (§ 3). We present the main results in § 4, and the conclusions in § 5.

## 2 Observations

We took advantage of the Spanish Observatorio del Teide (Tenerife, Spain) having two leading telescopes with sensitive spectro-polarimeters. The IR observations were carried out with the Tenerife Infrared Polarimeter (TIP, Martínez Pillet et al. 1999) at the German Vacuum Tower Telescope (VTT). We obtained spectra of all four Stokes parameters for the IR lines Fe I $`\lambda `$15648 Å (effective Landé factor, $`g_{eff}=3`$) and Fe I $`\lambda `$15652 Å ($`g_{eff}=1.53`$). The observations in the visible range were gathered with the spectro-polarimetric mode of the THEMIS telescope (López Ariste et al. 2000). We chose 3 different spectral ranges for the visible observations. Only one containing the lines Fe I $`\lambda `$6301.5 Å ($`g_{eff}=1.67`$) and Fe I $`\lambda `$6302.5 Å ($`g_{eff}=2.5`$) is analyzed here. The slit width was set to 0$`\stackrel{}{\mathrm{.}}`$5 in both telescopes, and we scanned the solar surface from east to west with a step size of 0$`\stackrel{}{\mathrm{.}}`$5. Each scan consisted of 60 positions, with an integration time of 30 sec per position.

The dataset for this work was obtained on August 10, 2002, in a quiet Sun region close to disk center ($`\mu =0.94`$). We selected IN regions trying to avoid bright features in the H$`\alpha `$ and Ca II K slit-jaw filtergrams supplied by the VTT. Having a few network patches in our fairly large Field-of-View (FOV) is unavoidable, though. In order to point both telescopes to the same region, we benefited from a video link that provided VTT slit-jaw images at THEMIS. First, we pointed both telescopes to an easily identifiable structure, e.g., a sunspot. Then we moved the VTT pointing to a quiet IN region and displaced THEMIS by the same distance. For the data set analyzed here, this method of differential pointing gave an offset of the co-alignment of the two telescopes of 1″. When the offset is corrected it leads to a time lag between the IR and visible data of only 1 minute. We found this offset by shifting the continuum intensity images obtained from the visible and the IR scans (Fig. 1). The optimum shift is chosen as that yielding the minimum difference between the two images.333Actually, we allow for a shift, a rotation of the FOVs, and a scale factor in the direction along the slit. The rotation is only 1°, whereas the scale factor is consistent with the ratio of effective focal lengths in the focal plane of the two polarimeters. Figure 1 shows the resulting intensities in the IR (left) and the visible (right). Despite the difference in angular resolution, it is clear that both images correspond to the same region. The white boxes point out a set locations on the Sun, and one can easily identify the same structures in the boxes of the two spectral ranges. The spatial resolution for the IR observations is 1″–1$`\stackrel{}{\mathrm{.}}`$2, with a granulation contrast of 1.7%, and it becomes 1$`\stackrel{}{\mathrm{.}}`$5–1$`\stackrel{}{\mathrm{.}}`$7 in the visible, with a contrast of 1.9%. The resolution was estimated as the cutoff of the power spectra of the continuum images shown in Fig. 1. The IR resolution is better due to the use of an active tip-tilt mirror that removes image motion during scanning (Ballesteros et al. 1996). The FOV common to visible and IR observations is of the other of 30″ $`\times `$ 35″ (Fig. 1).

The data reduction for the IR and visible ranges was very similar: flat field correction, removal of fringes, demodulation to get the Stokes parameters from the measurements, and correction of seeing-induced crosstalk by combining the two beams of each polarimeter. The instrumental polarization (IP) of the VTT was cleared away using a model Mueller matrix for the telescope, which was calibrated with large linear polarizers inserted at the telescope entrance (Collados 2003). THEMIS data do not need this step since it is an IP free telescope. The reduction process renders Stokes $`I`$, $`Q`$, $`U`$ and $`V`$ profiles in both spectral ranges for each pixel of our FOV. Here we only use intensity and circular polarization. The noise in the polarization measurements was estimated as $`34\times 10^4I_\mathrm{c}`$ in the IR data, and $`7\times 10^4I_\mathrm{c}`$ in the visible data ($`I_\mathrm{c}`$ is the continuum intensity). They correspond to the rms fluctuations of the polarization in continuum wavelengths.

## 3 Data analysis

We select for analysis those IR and visible Stokes $`V`$ profiles whose peak polarization is larger than three times the noise level. This thresholding separates clear signals from very noisy profiles or pure noise. The pixels with both IR and visible signals above their respective threshold correspond to 30% of the FOV. However, an additional 40% of the pixels present polarization signals only in either the visible or the IR. This study deals with those pixels with signals in both spectral ranges (some 1840).

All selected profiles were classified using a Principal Component Analysis (PCA) algorithm, similar to the one used by Sánchez Almeida & Lites (2000, § 3.2). We classify, simultaneously, the IR and visible Stokes $`V`$ profiles. All vectors are scaled to the maximum of the blue lobe of the Fe I $`\lambda `$6302.5 Å Stokes $`V`$ profile. The PCA classification results in 25 different classes with more than 0.5% of the classified pixels, and 15 classes containing rare profiles (between 0.5% and 0.2%). Figure 2 shows some of the significant classes. We note, firstly, that the IR profiles in class 0 are similar to the most abundant class obtained by Khomenko et al. (2003, Fig. 1). Secondly, class 2 corresponds to network structures in our FOV, and it has characteristic broad IR profiles (also observed in plage regions by Stenflo et al. 1987). Finally, the most conspicuous fact is that some classes (e.g., # 4 and # 8 in Fig. 2) show opposite polarity in the visible and in the IR lines (the lobes have opposite sign).

In order to estimate magnetic field strengths, we performed Milne-Eddington (ME) inversions of the profiles of all individual pixels, as well as the mean profiles among those in each class. We use MILK, a code which is part of the Community Inversion Codes444http://www.hao.ucar.edu/public/research/cic/ developed by HAO (see Socas-Navarro 2001). Each inversion provides 11 parameters, including the longitudinal component of the magnetic field strength, and the fill factor (fraction of resolution element covered by magnetic fields). (For a description of the procedure and approximations, see, e.g., Skumanich & Lites 1987.) The inversion is separately done for the two pairs of lines, i.e., we carry out two independent inversions per pixel, one for the visible pair and another for the IR pair. In each inversion we use the same atmospheric and magnetic field parameters for the two lines, whereas the ratio of absorption coefficients is fixed according to the expected values. Note that our ME inversion does not fit line asymmetries, although they are clearly present among the observed profiles.

## 4 Results

The PCA classification indicates that 25% of the pixels with coincident visible and IR signals give visible and IR profiles with opposite polarity. It follows that the visible and IR lines are tracing different magnetic structures in the same resolution element. Some classes have IR profiles with three lobes while the visible profiles have only two (not included in Fig. 2). Actually, the observed Stokes $`V`$ profiles present a wide variety of asymmetries, in agreement with previous observations (e.g., Sánchez Almeida & Lites 2000; Khomenko et al. 2003).

The ME inversions show that most of the classes have visible Stokes $`V`$ profiles characteristic of kG fields, while, simultaneously, the IR profiles indicate the presence of sub-kG field strengths. This also happens when we analyze all individual profiles. The result cannot be significantly biased by noise, since the mean profile of a class results from averaging hundreds of profiles and, therefore, bears negligible noise (see Fig. 2). On the other hand, the systematic difference between the IR and the visible is insensitive to the specific technique used to diagnose the field strength. The ratio between the Stokes $`V`$ extrema of the two visible lines can be used to diagnose the regime of field strength in the atmosphere. A ratio close to 1 indicates kG magnetic fields, whereas a value of some 0.6 is expected for sub-kG fields (see Sect. 5 in Socas-Navarro & Sánchez Almeida 2002). Most of the classes have ratios between 0.9 and 1, which evidences IN magnetic structures in the regime of kG fields. After the inversion of profiles in the various classes, and also from individual profiles, we found a mean magnetic field strength of 1100 G from the visible lines. They fill only 1% of the resolution elements. The mean field strength in the IR is $`300`$ G, and covers some 2% of the part of FOV studied here.

The ME inversions also show a few classes with weak magnetic fields both in visible and the IR. There are also classes with strong fields in both spectral ranges corresponding to network points in our FOV (class 2 in Fig. 2).

Longitudinal magnetograms measure flux density. One can infer it from ME inversions as the product of the longitudinal component of the field strength times the fill factor. When averaged over the 30% of the FOV that we study, the mean unsigned flux density turns out to be $`11`$ G in the visible, and $`6`$ G in the IR. Due to the complexity of the IN magnetic fields, the observed unsigned flux is expected to increase with increasing angular resolution (see, e.g., Sánchez Almeida et al. 2003; Domínguez Cerdeña et al. 2003b). The fact that the IR, with higher resolution, has lower unsigned flux, seems to indicate that this is a true solar effect. Specifically, most of the unsigned magnetic flux of the IN seems to be in the form of kG fields. On the other hand, the signed flux of the region in the visible is $`1.5`$ G, while it is less than $`0.2`$ G in the IR. A similar net flux has been observed before in visible lines (e.g. Lites 2002; Domínguez Cerdeña et al. 2003b). We find it to be co-spatial with very low signed flux in IR lines.

Using magnetic fields from the ME inversion, we have estimated the Probability Density Function (PDF) of finding a given longitudinal magnetic field strength in our IN quiet Sun region. We assume that the visible and the IR polarization signals trace different structures (a reasonable hypothesis, in view of the findings described above). The range of observed field strengths is divided in bins of 50 G. Then, treating visible and IR observations as independent measurements, we assigned each measured longitudinal magnetic field to one of the bins. The probability of finding magnetic field strengths in a given bin (i.e., the PDF) will be proportional to the sum of fill factors of all those measurements in the bin. The result is shown in Fig. 3, the solid line. This PDF has been normalized so that its integral is 0.3, i.e., the fraction of FOV under analysis. In other words, the normalization assumes that the remaining 70% of the surface is field free. Figure 3 includes the PDF for the IR data alone (the dashed line), and for the visible data alone (the dotted line). We also plot the PDF suggested by Socas-Navarro & Sánchez Almeida (2003, the dash-dotted line). We find that 75% of the total flux is concentrated in fields stronger than $`500`$ G, which, however, occupy only 0.4% of the full FOV.

A final comment is in order. The visible and IR maps that we employ have a time lag of one minute and were coaligned with a finite accuracy (however better than the angular resolution). The lack of strict co-alignment and simultaneity may produce some differences between the visible and the IR spectra, however, they would be of random nature and, consequently, they can hardly explain the observed systematic differences.

## 5 Conclusions

The (almost) simultaneous visible and IR Stokes $`V`$ spectra described above have several conspicuous properties. They reveal the coexistence of kG and sub-kG magnetic fields in most of the pixels with polarization signals above the noise (30% of the FOV). The IR lines indicate sub-kG fields whereas visible lines trace kG. This demonstrates the conjecture by Sánchez Almeida & Lites (2000) and Socas-Navarro & Sánchez Almeida (2002, 2003).

At least 25% of pixels with coincident visible and IR signals have one polarity in the visible and the opposite in the IR. This second result offers a direct proof of the existence of two polarities in 1$`\stackrel{}{\mathrm{.}}`$5 resolution elements. They have been inferred from line asymmetries (Sánchez Almeida & Lites 2000), as well as from the skewness of the distribution of line-of-sight magnetic field inclinations in measurements away of the solar disk center (Lites 2002).

The Probability Density Function for the longitudinal magnetic fields, that we infer from ME inversions, has a long tail of kG fields. Some 75% of the observed (unsigned) magnetic flux is in field strength larger than 500 G. This excess of kG flux is in agreement with the PDF suggested by Socas-Navarro & Sánchez Almeida (2003) to reconcile discrepancies between field strengths obtained with different methods (see § 1). We have also found the extended tails in the IR line Fe I $`\lambda `$15648 Å predicted by such PDF (see classes 0 and 5 in Fig. 2). They carry the information on the strong fields, but are easily misinterpreted due to the large lobes formed by weak fields.

Our observations prove that the use of only visible lines or only IR lines to study the IN fields biases the results. Diagnostics based on simultaneous visible and IR observations offer a viable and reliable alternative.

We thank M. Collados for training with TIP operation, and for assistance during the data reduction. Thanks are also due to H. Socas-Navarro, who provided the IC, to T. Eibe, who supported us with TIP at the VTT, and to the team of support astronomers and operators of THEMIS. Comments on the manuscript by E. Landi Degl’Innocenti and the referee were very useful. IDC acknowledges support by the Deutsche Forschungsgemeinschaft (DFG) through grant 418 SPA-112/14/01. The VTT is operated by the Kiepenheuer-Institut für Sonnenphysik, Freiburg, and the French-Italian telescope THEMIS is operated by CNRS-CNR, both at the Observatorio del Teide of the Instituto de Astrofísica de Canarias, Spain. The work has partly been funded by the Spanish Ministry of Science and Technology, project AYA2001-1649.


text:
# The Pierre Auger Observatory

## I Introduction

The observation of ultra high energy cosmic rays (UHECR) with energies above $`10^{20}`$ eV ev , is puzzling in at least two ways. First of all, it is quite difficult to conceive of acceleration mechanisms which are adequate to impart such enormous energies to cosmic ray particlesbh ; ol  – energies comparable to those carried by everyday objects like tennis balls or golf balls! Secondly, even if such a mechanism is found, it is difficult to see how such high energy particles would make it through the background radiation: the celebrated GZK cutoffgzk  predicts that protons over about $`5\times 10^{19}`$ eV should rapidly lose energy in inelastic collisions with the cosmic microwave background photons with similar energy degradation mechanisms being present for most other particles, including heavy nuclei.

A very important point to make early in this talk is that while there is a puzzle concerning so-called “super-GZK events”, the Pierre Auger Obseratory will provide interesting information regardless of how well the GZK cutoff holds out. If there is no GZK cutoff, this will be an unambiguous sign of new physics. If the GZK cutoff is found, then we can rest assured of that piece of the physics and confidently use the data to try to understand the nature of the sources. This point is often missed by people who see the whole physics motivation as “is there a GZK cutoff or not?”.

Above $`10^{15}`$ eV, cosmic ray primaries are not detected directly, but rather through the effects that such particles produce when they strike the upper atmosphere. There they initiate a cascade of reactions, some nuclear, but most forming an electromagnetic shower made of repeated bremsstrahlung and $`e^+e^{}`$ pair creation events. This shower can be detected experimentally through the fluorescence it produces in the atmosphere ( due to excited nitrogen) or via the particles that reach the ground.

Our understanding of the highest energy part of the spectrum above the so-called “ankle” ($`5\times 10^{18}`$ eV) is poor due to a combination of low statistics, uncertain energy resolution, uncertainties in energy conversion arising from models, and a lack of knowledge of the mass composition and the fluorescence yield efficiency. Clearly there is something interesting going on, but to fully understand the situation we need more statistics, better control over systematic uncertainties, and full sky coverage: enter the Pierre Auger Obervatory. A comprehensive review of the state of the art prior to the Pierre Auger Observatory can be found in bigreview .

## II The Pierre Auger Observatory

The Pierre Auger Observatory (PAO) is actually comprised of two sub-observatories, one currently under construction in Mendoza, Argentina since 2000. This site is especially interesting, since in addition to being in the wine-growing district, it offers a view of the centre of the galaxy. Another is planned for the Northern hemisphere, and while for the remainder of this talk I will concentrate on the Southern site it is important to understand that the full observatory is comprised of two sites. This will eventually allow full-sky coverage which is very important to allow good studies on anisotropiesanisotropies .

The PAO is designed to measure the energy, arrival direction and primary species with excellent precision and very high statistics. A unique feature of the design is the combination of both fluorescence detection and ground-based particle detectors which can be operated independently as well as together in “hybrid” mode.

The scale of the observatory was determined by the requirement that we can collect high statistics in and around the region of the expected GZK cutoff, with 1600 particle detectors separated from each other by 1.5 km and covering an area of 3000 km2, overlooked by four fluorescence detectors which can only operate when ambient conditions offer a clear, dark sky, which leads to a roughly 10% duty cycle. Figure 1 shows, together with one of the surface detectors, a photograph of one the fluorescence detectors where both the mirror and the box of 440 photomultipliers which register the fluorescence light can be seen. The fluorescence measurements are complemented by a a very comprehensive atmospheric monitoring system. The surface array stations are water Čerenkov detectors and can operate continuously.

The fact that about 10% of the showers detected by PAO will be observed by both surface and fluorescence detectors, offers the possibility of doing calibrations and understanding systematic errors in a manner that has never been possible before. Access to a large-dimensional parameter space of observables should allow not only determination of the direction and energy of incoming primaries, but also the disentanglement of information about composition from the notoriously difficult systematic errors associated with the choice of hadronic interaction models.

Each ground-based detector is a cylindrical, opaque tank of 10 m2 and a water depth of 1.2 m, where particles produce light by Čerenkov radiation. The filtered water is contained in a bag which diffusely reflects the light collected by three photomultipliers (PMT’s) installed on the top. The large diameter PMT’s ($``$ 20 cm ) are mounted facing down and look at the water through sealed polyethylene windows that are integral part of the internal liner. The signals are processed locally and a second level trigger is identified before transmitting the data to the central acquisition system Tiina . The fact that the tanks are quite deep enables showers to be detected efficiently over a wide angular range. Due to the size of the array, the stations have to be able to function independently and yet in communication with the central data acquisition system. The stations operate on battery-backed solar power and communicate with a central station by using wireless LAN radio links Paul . Absolute timing information is obtained from the Global Positioning System (GPS) clejer  and is used to reconstruct the direction from which the primary came. Figure 1 shows a water Čerenkov detector installed in the Southern Observatory as well as one of the fluorescence detectors. Mounted on top of the tank are the solar panel, electronic enclosure, mast, radio antenna and GPS antenna. The battery is housed in a box attached to the tank.

The expected angular resolution for the ground array of the Southern Auger Observatory is less than 1 for all energies, and better for large events above $`10^{20}`$eV. The expected energy resolution is estimated to be 12%, averaged over all energies (assuming a proton-iron primary mixture), falling to 10% at $`10^{20}`$eV. The limiting aperture for the full Southern Observatory array and for zenith angle less than 60 is 7350 km2sr. The detection efficiency at the trigger level should reach 100% for energies above $`10^{19}`$eV Jere . Additionally, if events above 60 degrees can be analyzed effectively, the aperture will increase by about 50% .

In hybrid mode, the Pierre Auger Observatory is expected to have 6% energy resolution and an angular precision of 0.5 at $`10^{20}`$eV where only statistical errors are taken into account in these estimates. The detector is optimized for energies above $`10^{19}`$eV, with good reconstruction expected at energies down to 1 EeV. The hybrid data set will provide the best evaluation of primary species, allowing a simultaneous fit to all parameters sensitive to mass composition.

The first cosmic ray event detected by one of the two prototype telescopes installed at Los Leones is displayed in Figure 2. A twenty pixel track, produced by light from a shower, with a length of 8 $`\mu `$s can be seen. The angular velocity of the shower image across the sky allows the distance of the shower core to be established as 5 km. The time duration of the signal in the field of view of the telescope corresponds to a track length of 2.4 km. The mirror inverts the image, and particles from the sky appear to be going upwards as seen by the camera.

Construction continues, and we look forward to the large amount of very high quality which will appear in the years to come, and the light it will shed on one of nature’s great mysteries.

## Acknowledgments

I would like to thank the organizers of this conference for a most interesting meeting, as well as all my collaborators in the Pierre Auger Observatory. Special thanks are due to Maria Teresa Dova for assistance in the preparation of this talk, and, as always, I would like to thank the National Science Foundation for its continued support.


text:
# The RHMC Algorithm for 2 Flavours of Dynamical Staggered Fermions

## 1 Introduction

Computations with two flavours of dynamical staggered quarks are quite popular at present. There are a number of possible problems with such calculations such as flavour symmetry breaking and non-locality of the square-root of the four-flavour action. In this investigation we shall ignore these and consider only the possible errors introduced through algorithmic approximations.

We propose the use of the Rational Hybrid Monte Carlo (RHMC) algorithm . This method is stochastically exact, in the sense that it is free from molecular dynamics (MD) stepsize errors. It is comparable to the usual $`R`$ algorithm  in performance, but without the need for extrapolation in the MD stepsize $`\delta t`$.

## 2 Two Flavour Algorithms

All Hybrid Molecular Dynamics (HMD) algorithms have the same underlying structure: a fictitious momentum field is introduced, and the gauge field is integrated along classical trajectories in fictitious time, interleaved with refreshment of the momenta from a Gaussian heatbath. When integrating Hamilton’s equations the evaluation of the fermionic contribution to the force acting on the gauge fields is the costliest part of generating full QCD gauge field configurations. Most algorithmic developments are techniques to calculate the fermionic force more efficiently.

The desired probability distribution for the gauge fields $`U`$ with $`n_\mathrm{f}`$ flavours of staggered fermions is

$$P(U)e^{S_G(U)}det\left[(U)\right]^{n_\mathrm{f}/4},$$
(1)
where $``$ is the staggered fermion kernel and $`S_G`$ is the gauge action. Thus for $`n_\mathrm{f}=2`$ flavours of fermion we require the square root of the fermion determinant. Choosing a suitable normalisation the spectrum of the staggered fermion kernel $``$ is contained in the interval $`[\epsilon ,1]`$, where $`\epsilon =\left(1+\frac{16}{m^2}\right)^1=\frac{m^2}{16}+O(m^4)`$.

### 2.1 The $`R`$ algorithm

The identity $`det^{n_\mathrm{f}/4}=\mathrm{exp}trln^{n_\mathrm{f}/4}`$ allows us to express the determinant as a term in the action, and the number of flavours just appears as a factor in front of this term, $`S_{\text{PF}}=\frac{n_\mathrm{f}}{4}trln`$. In the $`R`$ algorithm  the evaluation of the force corresponding to this trace is performed stochastically, as computing it exactly would be prohibitively expensive. To do this without introducing $`O(\delta t)`$ errors we introduce an auxilliary field that is evaluated at time $`(1\frac{n_\mathrm{f}}{4})\frac{\delta t}{2}`$ along each integration step. For two flavours this breaks time-reversal invariance, violates Liouville’s theorem, and leads to an irreversible and non area-preserving algorithm: as such, it cannot be made exact by the inclusion of a Metropolis step. A detailed analysis of the errors in the probability distribution produced by this algorithm was presented in .

## 3 Rational Hybrid Molecular Dynamics

RHMD  uses a uses a rational approximation to fractional powers $`\alpha `$ of a matrix. This is analogous to the use of polynomial approximations introduced in , but optimal (Chebyshev) rational approximations give a much closer approximation over a given interval than the corresponding polynomials of similar degree (Fig. 1). The best choice for this rational approximation would seem to be given by a relative minimax approximation: this has the property that it is guaranteed to deviate by at most an amount $`\mathrm{\Delta }`$ relative to the correct value over a given interval, $`max_{\epsilon x1}\left|1x^\alpha r(x)\right|=\mathrm{\Delta }`$. The errors are of the same kind as those introduced by the use of floating-point arithmetic. The maximum error $`\mathrm{\Delta }`$ falls exponentially with the degree of the rational approximation used.

We approximate the determinant in (1) as a Gaussian integral over a bosonic (pseudofermion) field $`\varphi `$,

$$det^{n_\mathrm{f}/4}detr()^2𝑑\varphi 𝑑\varphi ^{}e^{\varphi ^{}r()^2\varphi },$$

where $`1r()^{n_\mathrm{f}/8}\mathrm{\Delta }`$ in the spectral norm. The RHMD algorithm proceeds as follows

* A momentum refreshment heatbath using Gaussian noise: $`P(\pi )e^{\frac{1}{2}|\pi |^2}`$.
* A pseudofermion heatbath using Gaussian noise: $`\varphi =\left[r\left((U)\right)\right]^1\xi `$, where $`P(\xi )e^{\frac{1}{2}|\xi |^2}`$. Note that the inverse $`1/r(x)`$ is itself a rational function.
* An MD trajectory consisting of $`\tau /\delta t`$ steps with Hamiltonian $`\overline{H}=\frac{1}{2}|\pi |^2+S_G+\varphi ^{}\overline{r}()\varphi `$ with $`1\overline{r}()^{n_\mathrm{f}/4}\overline{\mathrm{\Delta }}`$ in the spectral norm.

This leads to an algorithm which has finite stepsize errors of $`O(\delta t^2)`$ and errors of $`O(\overline{\mathrm{\Delta }})`$ and $`O(\mathrm{\Delta })`$ incurred from the use of rational approximations.

Rational functions can be expressed as a product or as a partial fraction expansion,

$$\overline{r}(x)=\overline{c}_0\underset{i=1}{\overset{d}{}}\frac{(x\overline{\gamma }_i)}{(x\overline{\beta }_i)}=\overline{\alpha }_0+\underset{i=1}{\overset{d}{}}\frac{\overline{\alpha }_i}{x\overline{\beta }_i}.$$

In partial fraction form  the pseudofermionic force takes the form

$$\frac{\overline{H}}{U}=\varphi ^{}\frac{\overline{r}()}{U}\varphi =\underset{i=1}{\overset{d}{}}\overline{\alpha }_i\chi _i^{}\frac{}{U}\chi _i$$

where $`\chi _i=(\overline{\beta }_i)^1\varphi `$. A multishift solver  can be used to compute all the $`\chi _i`$ in a common Krylov space. The computational cost of generating the appropriate Krylov space depends upon the smallest shift, and the only extra cost is that of updating the extra $`d1`$ solution vectors.

## 4 Rational Hybrid Monte Carlo

The RHMC algorithm is similar to RHMD but with the addition of a Metropolis accept/reject step. The acceptance probability for this is given by

$$P_{\text{ACC}}=min(1,e^{\delta H}\frac{det_{}^{}{}_{}{}^{n_\mathrm{f}/4}r(^{})^2}{det^{n_\mathrm{f}/4}r()^2})$$
(2)
where $`H=\frac{1}{2}|\pi |^2+S_G+\varphi ^{}r()^2\varphi `$.

If we use a high enough degree rational approximation $`r`$ such that $`\mathrm{\Delta }1ulp`$ (unit of least precision for the floating point arithmetic used) in the computation of the pseudofermion heatbath and of $`\delta H`$ we can ignore the explicit determinants in (2) without introducing any systematic errors beyond the ever-present rounding errors. In practice we find that it is sufficiently cheap to use the same (machine accuracy) rational approximation for the MD integration as well ($`\overline{\mathrm{\Delta }}=\mathrm{\Delta }`$), although this is not logically required.

## 5 Noisy Rational Hybrid Monte Carlo

RHMCN  allows the use of a lower degree rational approximation $`r`$ while keeping the algorithm exact. The algorithm replaces the Metropolis step with a Kennedy-Kuti (KK)  noisy accept/reject step. A stochastic summation is used to estimate the determinant ratio in (2). The acceptance probability is defined as $`P_{\text{ACC}}=\lambda _\pm +\lambda _{}Q(U,U^{})`$ for $`U>U^{}`$ and $`U<U^{}`$ respectively, where $`Q(U,U^{})`$ is an unbiased (noisy) estimator of the determinant ratio occurring in (2) and $`\lambda _\pm `$ are parameters used to ensure the resultant probability distribution lies in the range $`[0,1]`$. When $`\lambda _+=\lambda _{}=\lambda `$ the average acceptance rate is $`P=2\lambda `$.

## 6 Comparison of RHMC and $`R`$

The use of multishift solvers in the implementation of the rational algorithms results in a lower computational cost than might otherwise be expected. On a $`16^3\times 32`$ lattice with $`\beta =5.26`$, $`m=0.01`$, and $`\delta t=0.01`$ the time to perform one trajectory of length $`\tau =0.5`$ on a single node Pentium 4 processor with no assembler optimisations is 274 minutes for the $`R`$ algorithm. RHMC takes 318 minutes when a degree 10 rational function is used. Although a single $`R`$ trajectory is faster, RHMC compares very favourably when it is remembered that an $`O(\delta t^2)`$ extrapolation ought to be carried out when using the inexact $`R`$ algorithm.

The computational cost of using RHMCN is unfavourable compared with RHMC. Although a lower degree rational approximation can be used this does not give a large benefit because of the efficacy of the multishift solver. To ensure negligible probability violations occur when $`\mathrm{\Delta }>1ulp`$ the KK acceptance test becomes increasingly expensive for large volumes and small quark masses. The cost of performing a single trajectory using the parameters given above with a degree 6 rational approximation and typical stochastic summation parameters is 385 minutes. We found that setting $`\lambda `$ to give $`P=70\%`$ gives as large an acceptance rate as feasible without incurring excessive violations, but for large $`V`$ and small $`m`$ it was found necessary to reduce this to about $`50\%`$.

## 7 Conclusion

We have found that it is easy and cheap to compute rational powers of the staggered fermion kernel to within machine accuracy using Chebyshev rational approximations expressed as partial fractions, and applied using a multishift solver. This form of the RHMC algorithm thus enables the exact HMC algorithm to be extended to the case of an arbitrary number of flavours. The ability to use Krylov space solvers makes RHMC faster than the PHMC algorithm . In terms of computational cost there is very little to chose between the $`R`$ and RHMC algorithms. Since this seems to be the only possible advantage of $`R`$ over RHMC, we conclude that there is no reason for the continued use of the $`R`$ algorithm.

## Acknowledgments

We would like to thank Robert Edwards, Ivan Horváth, Bálint Joó, and Stefan Sint for discussions.


text:
# A study of Complexity in Gamma Ray Burst using the Diffusion Entropy Approach

## 1 Introduction

Gamma-Ray Bursts (GRB), discovered in 1973 by the Vela satellites (Klebesadel et al. kle73 (1973)), are among the most powerful phenomena in the universe. Their typical emission is of the order of $`10^{52}`$ erg and the nature of these events seems to be destructive, maybe related to explosions. The energy emission during the flaring phase (prompt emission) is in the gamma-ray range, with a $`\nu F_\nu `$ peaked around few hundreds of keV.

The GRB durations are distributed following a bimodal distribution covering three decades (Kouvelietou et al. kou93 (1993)). Short bursts last less than milliseconds while the longest ones can last up to hundreds of seconds. Moreover, GRB temporal series (light curves) are typically very different from one burst to another. The light curves have many peaks, and in spite of extensive statistical studies, the temporal behavior of GRB remains a puzzle. Contrary to the very diverse behavior in the time domain, GRB have a simple behavior in the Fourier domain: the Power Density Spectrum, averaged over many bursts, shows a universal power law decay $`1/f^\beta `$ (with a break at $`1Hz`$) with $`\beta 5/3`$ (Beloborodov et al. bel00 (2000)). This characteristic decay has been observed in many other complex systems such as fully developed turbulence in a fluid flow (Frisch fri96 (1996)), self-organized criticality (Bak et al. bak87 (1987)), and, more generally, in systems with long-range correlations. One of the methods used to study the “complexity” of these systems is Diffusion Entropy (DE), introduced by Scafetta et al. (scaf01 (2001)). In our analysis we make use of the available light curves with 64 ms time resolution obtained by the Burst and Transient Source Experiment (BATSE) in the four Large Area Detector (LAD) energy channels.

The paper is organized as follows: in Sect. 2 we describe the method, in Sect. 3 we consider different analyses by applying the method to the BATSE light curves. We first apply the DE method to two GRB that show quite different light curves. We then apply the DE method to the long GRB population (with $`T902`$ seconds), building up the statistics of the scaling indices in the four energy channels. We find a relationship between the scaling index $`\delta `$ and the power measured in GRB emission (fluence). We show, with the support of simple simulations, how this relation represents a key study for constraining models. Finally, in order to investigate the temporal evolution of the complexity in the GRB light curve, we study the Diffusion Entropy as a function of time for the two selected bursts of the first analysis. Our conclusions are drawn in Sect. 4.

## 2 Diffusion Entropy Method

Diffusion Entropy was introduced in Scafetta et al. (scaf01 (2001)) and it has been applied to many different fields and, in astrophysics, to Solar Flares (Grigolini et al. grigo02 (2002)). It is based on the central limit theorem and allows the correlation of a given sequence to be measured (Feller fel71 (1971)). In the standard language of the diffusive processes, given a signal $`\left\{x_i\right\}_{i=1,N}`$ we can think of each $`x_i`$ as the length of a jump of a walker moving in a one dimensional space. We create a set of many different diffusion trajectories by means of a moving window of size $`\mathrm{\Delta }t`$ with $`1<\mathrm{\Delta }t<N`$ which is the same as a boxcar smoothing forward in time from point $`t`$ to any point $`t+\mathrm{\Delta }t`$. In practice, we generate $`N\mathrm{\Delta }t+1`$ trajectories $`y_j`$ considering the sums:

$$y_j(\mathrm{\Delta }t)=\underset{i=j}{\overset{j+\mathrm{\Delta }t}{}}x_i.$$
(1)
Each $`y_j(\mathrm{\Delta }t)`$ can be thought of as the final position of the walker after $`\mathrm{\Delta }t`$ time steps. The normalized histogram that contains the final positions at time $`\mathrm{\Delta }t`$ ($`y_i(\mathrm{\Delta }t)`$) is the probability distribution $`p(y,\mathrm{\Delta }t)`$ of the final positions $`y`$ of the walker at time $`\mathrm{\Delta }t`$ for which the Shannon entropy is:

$$S(\mathrm{\Delta }t)=_{\mathrm{}}^+\mathrm{}p(y,\mathrm{\Delta }t)\mathrm{ln}(p(y,\mathrm{\Delta }t))𝑑y.$$
(2)
From the central limit theorem, if the time series $`\left\{x_i\right\}_{i=1,N}`$ has finite variance and is uncorrelated, the probability $`p(y,\mathrm{\Delta }t)`$ satisfies the scaling condition:

$$p(y,\mathrm{\Delta }t)=\frac{1}{\mathrm{\Delta }t^\delta }F(\frac{y}{\mathrm{\Delta }t^\delta });$$
(3)
then $`\delta =\delta _\text{N}0.5`$ and $`F`$ is the Gaussian function for normal (Brownian) diffusion (Feller fel71 (1971)).

More generally, a larger class of diffusive processes can be characterized by the scaling condition (3) with $`\delta \delta _\text{N}`$ and non-Gaussian $`F`$. Such processes have been widely observed in nature (Shlesinger et al. shles93 (1993)) and are known collectively as “anomalous” processes. For these the central limit theorem hypothesis no longer holds: the signal $`x_i`$ maintains long time correlation (not Markovian) and/or it is not stationary. In other words, the DE algorithm smoothes the signal by summing a series of data contained in the window $`\mathrm{\Delta }t`$ and characterizes how the information loss increases with increasing window size. The rate of information loss for an “anomalous” process is faster than for a normal diffusion. The DE method provides a discrimination between normal and anomalous diffusion by determining the scaling index $`\delta `$. When the scaling condition (3) holds, the Shannon entropy $`S(\mathrm{\Delta }t)`$ can be written in the general form

$$S(\mathrm{\Delta }t)=A+\delta \mathrm{ln}(\mathrm{\Delta }t),$$
(4)
independent of the value of $`\delta `$ and of the shape of the function $`F`$ (Scafetta et al. scaf01 (2001)).

As a first step we apply DE method to the “synthetic” light curve, where the counts are distributed according to the Poisson distribution. We set the binning of the data to 64 ms and we generate 200 seconds of signal. As expected from the theory, the $`\delta `$ index is $`\delta _\text{N}=0.5`$. We will refer to this reference signal as “noise”. Fig. 1 shows the noise signal (top) and the corresponding entropy $`S(\mathrm{\Delta }t)`$ as a function of the window $`\mathrm{\Delta }t`$ (bottom).

In the figure each point corresponds to the value of the Shannon entropy $`S(\mathrm{\Delta }t)`$ as a function of the size $`\mathrm{\Delta }t`$ of the window, varying from 64 ms to 20 seconds with steps of 64 ms. The choice of this value comes from the bin size of the experimental data and from the duration of the light curve. We choose the maximum size of the window as $`1/10`$ of the total duration of the light curve. This condition is required to have enough trajectories $`y_i(\mathrm{\Delta }t)`$ to compute the probability $`p(y,\mathrm{\Delta }t)`$. However, the peculiar feature of GRB time series, i.e. the short time duration of the prompt emission with respect to the whole data set, makes the analysis more difficult. Indeed, the index $`\delta `$ measures the non-stationarity/memory “intensity” of the light curve, and, therefore, the diffusion process is governed by the balance between the GRB signal and the background noise.

Notice that, if we consider the whole BATSE data series, each light curve is composed of several tens of seconds of noise and the GRB signal over-imposed. If, for instance, the GRB is short and weak, the measured $`\delta `$ will be close to $`\delta 0.5`$. From this point of view, the DE could mistake a signal with a small contribution of non-stationary behavior for a noisy time series. For this reason, we restrict our analysis to the long burst population by selecting only those bursts with $`T902`$ seconds. This choice ensures that $`\delta `$ indices are uncorrelated with respect to the GRB duration and it also ensures that the measure of the DE is uncorrelated from the experimental acquisition time.

## 3 Diffusion Entropy and BATSE data

We apply the DE method to data available in the CGRO/BATSE archive111ftp://cossc.gsfc.nasa.gov/compton/data/batse/ascii\_data/64ms/. These data are 64 ms binned time series, containing the counts for the four LAD energy channels, 1-4: (1) 20-50 keV, (2) 50-100 keV, (3) 100-300 keV, (4) 300-1000 keV. We perform three different studies. We first compute the DE of two very different GRB: the scaling indices $`\delta `$ are indeed different from the reference value of noise signals $`\delta _\text{N}`$. We compute the DE indices for the entire catalog, and we study the dependence of the $`\delta `$ indices on the four energy channels. We then study the correlation between the $`\delta `$ indices and the fluence, finding the relationship between the complexity and the power in GRB. We also propose a schematic model to understand this relation. Finally we show that DE varies during the GRB time evolution.

### 3.1 The case of a single burst

In Fig. 2 we plot the count rate of two different GRB (910429, top and 910601, bottom) as a function of time for each of the four BATSE energy channels.

The two selected burst exhibit different characteristics in terms of temporal variability, and peak intensities. The first is a faint, smooth burst while the second is a strong, spiky burst. We adopt the model that the structure of a GRB can be viewed as the non-stationary increase of the probability of a walker to make longer jumps. Since the activity of the GRB lasts only a finite amount of time, this is what we have called non-stationarity. This behavior is evident in the first three channels of GRB 910429 (lowest energies) where the signal exceeds the noise background. In the fourth channel the signal seems to be dominated by the noise, and the not stationarity of the signal is less evident. What we mean by complexity is essentially the same as anomalous diffusion: it is an estimate of how far an anomalous diffusive process is from completely uncorrelated noise. In the case of GRBs the high complexity arises from the high time variability, the nonperiodicity, and the typical non-stationarity of the signal. Yet our claim is the complexity of a GRB signal derives from the highly structured environment that surround the burst. Indeed, there are reasons to think that the high power released in the GRB phenomena, combined with the high variability, requires a site where very efficient energy conversion processes take place in a structured, perhaps clumpy, environment with substantial variation from one source to another. The measurement of the Diffusion Entropy by means of the scaling index $`\delta `$, is indeed an estimation of the complexity of the GRB light curve: the difference between $`\delta `$ and $`\delta _\text{N}`$ quantifies the “distance” between the light curve and the noise.

The results of the application of the DE algorithm to the time series of Fig. 2 are shown in Fig. 3. The lines are the fits to the data with the function $`A+\delta \mathrm{ln}(\mathrm{\Delta }t)`$ for the different energy channels. The variable $`\mathrm{\Delta }t`$ is the size of the temporal window for which the Shannon entropy $`S(\mathrm{\Delta }t)`$ is evaluated. The temporal scale over which the statistics is “good” is wide since we analyze an order of $`N10^3`$, corresponding to $`\mathrm{\Delta }t`$ from 64 ms to 24 seconds in 375 steps. The resulting indices are shown in Fig. 3. GRB 910601 shows higher values in each channel than GRB 910429, reflecting the greater complexity of the former.

In order to estimate the error on the $`\delta `$ indices due to Poisson fluctuations of photon counting of the original data, we proceeded as follows. A set of fictitious light curves $`\left\{x_i^j\right\}_{i=1,N;j=1,L}`$ was generated from an original $`\left\{x_i^0\right\}_{i=1,N}`$ light curve, where $`x_i^j`$ is a random number extracted from a Poisson distribution of mean $`x_i^0`$. We find a Gaussian distribution of $`\delta `$ indices. The dispersion of which gives an estimate of the uncertainty in $`\delta `$. This analysis with $`L=500`$ gives an error $`\delta _{\text{err}}0.05`$. Therefore, studying the $`\delta `$ indices we can discriminate a signal from the noise also in “critical” situations, as in GRB 910429 (top) channel 4, where the signal-to-noise ratio is close to one. Indeed, in this particular case, $`\delta =0.61\pm 0.03`$. Even if the fourth channel series is dominated by stationary noise, it has a non-stationary component (i.e. a memory effect): the computation of the scaling index $`\delta `$ reveals this “anomaly” and gives a unique measure of its value. We notice that in many cases the background is itself non-stationary due to, for instance, fluctuations of instrumental noise or the presence of other transient sources in the field of view of the instrument (like pulsars). In these cases, to treat only the GRB contribution to the non-stationarity, we have to detrend the time series. To do this, we fit the whole time series with a second order polynomial, taking into account only the portion of the data set far from the GRB prompt emission (i.e. the ending tail of the data series). We then remove the obtained polynomial background from the original data set.

### 3.2 DE statistics of the GRB catalog

In the previous section we have shown how the complexity of a GRB can be characterized by the scaling index $`\delta `$. For GRBs we notice that the DE is correlated both with the intensity of the signal and with the duration of the prompt emission. The analysis of the entire catalog can give us an indication of how the two quantities (duration and intensity) affect the value of the $`\delta `$ indices. The DE index scales linearly with the $`T90`$ parameter. This trend is due to the balance between noise and GRB signal given by the ratio of the duration of the noise to the duration of the signal. More interesting is the relationship between the DE and the intensity of the signal. To remove the correlation with the duration of the signal, we select a portion of the light curve proportional to the the duration of the GRB. We cut the data series at twice the $`T90`$ value, to ensure that only the GRB signal is analyzed. Short bursts cannot be analyzed, since the portion of the light curve equal to 2 $`T90`$ is too short and the measurement of the scaling index $`\delta `$ will be drastically affected by the lack of data. On the other hand, considering a longer portion of the light curve would exclude the longest bursts, for which the duration is similar to the duration of the data sets. We then apply the DE analysis to the population of long GRB (1224 drawn from the BATSE catalog). The distributions of long bursts are shown in Fig. 4 for the four energy channels.

The mean values of $`\delta `$ for the whole catalog in the four energy bands are $`0.80`$, $`0.82`$, $`0.79`$, and $`0.62`$ respectively. The RMS are approximately the same in the four channels and are indicated in the upper right corner of each panel. The distributions are, in general, not symmetric, and the $`\delta `$ indices are distributed in the range between 0.4 and 1. The mean scaling parameter $`\overline{\delta }`$, and, more generally, the shape of the distributions are very similar for the first three channels, while $`\overline{\delta }`$ is much smaller in the fourth channel. From this analysis, the first three channels seem to have the same degree of complexity, while the contribution of the non-stationarity/memory is, on average, weaker in the fourth channel. The number of counts in the highest energy range is small compared to the number of counts in channel 1,2, and 3. The case of completely correlated noise, also known as ballistic motion for which the walker always jumps in the same direction with a constant jump length, is the upper bound for the $`\delta `$ value (=1). The distribution of the values in the first three channels indicates that the diffusive process is close to the ballistic case. This means that, in general, GRB signals are highly correlated signals.

### 3.3 The power-complexity relation

In order to understand how the degree of complexity depends on the intensity, we perform another analysis. We consider the time series resulting from summing the time series in the four channels, using the sum of the fluences in the four channels as indicator of the intensity of the burst. We call $`\delta _{\text{tot}}`$ the $`\delta `$ index resulting by analyzing this time series. In Fig. 5 it is shown the scaling parameter $`\delta _{\text{tot}}(F)`$ of the light curve, as a function of the total fluence (filled dot) and the scaling parameter $`\delta _i(F)`$ for the four different channels as a function of the total fluence. Notice that the fluence is a measure of the energy emitted by the burst; thus, even if the count rate in the fourth channel is lower than in the first three channels, the fluences are similar. The plot shows that, at low fluences, the dominant contribution to the diffusion entropy is given by the first three channels where the non-stationarity of the light curve is already pronounced and the complexity is already developed. In the fourth channel the values of the $`\delta `$ index are closer to $`0.5`$. As the fluence increases, the complexity in channel four increases and at $`F10^4\text{ erg/s}`$ is comparable with the values of the first three channels. The balance between the increment of the complexity with the fluence in the four channels results in the logarithmic relation found by fitting the $`\delta _{\text{tot}}`$ as a function of the total fluence $`F`$.

We find a logarithmic dependence of the complexity of the data with respect to the fluence of the GRB: the higher the measured power, the greater the complexity of the GRB system. We notice that this result qualitatively agrees with the variability- luminosity relationship reported by Fenimore & Ramirez-Ruiz (feni00 (2000)), and Reichart et al. (reic01 (2001)). The fluence is a measure of the energy conversion efficiency in the GRB source, so the relationship between $`\delta `$ and the fluence can be interpreted as the increase of the complexity of the energy conversion process. This new relationship can be used as a testing ground for GRB scenarios or to constrain models, as, for instance, the internal shocks picture (Rees & Meszaros rees94 (1994), Kobayashi, Piran & Sari kob97 (1997), Daigne & Mochkovitch daig98 (1998), Piran piran99 (1999)), or the external shock scenario (Rees & Meszaros rees92 (1992), Chevalier & Li chev99 (1999), Dermer et al. derm95 (1995)).

### 3.4 GRB simulations: a simple (Poisson) model

We developed a schematic model by defining a pulse shape, and assuming that a GRB is made by the superposition of several pulses. We chose the simplest model by limiting the number of parameters, so we adopted a schematic fixed pulse shape: the simulated counts as a function of time are expressed by an exponential function:

$$F(t)=\{\begin{array}{cc}A\mathrm{exp}((tt_\text{p})/t_\text{d})\hfill & tt_\text{p}\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$$
(5)
the pulse shape $`F(t)`$ depends on the value of $`A`$ which gives the amplitude of the spikes, on the peak time $`t_\text{p}`$ which gives the peak center, and on the time constant $`t_\text{d}`$. The duration of the spikes is given by the relation:

$$T90=t_\text{d}\mathrm{ln}0.1,$$
(6)
obtained by computing the temporal window which contains $`90\%`$ of the spike. For a single spike, the value of the index $`\delta `$ increases with $`A`$, and with the duration of the pulse. To reproduce the variability observed in GRB light curves we randomly extracted the amplitude of the spike (varied the value of $`A`$) and we summed several peaks by extracting the waiting time between pulses from an exponential law with time constant $`\tau `$. To simulate analogous count rates we simulated the signal in the four channels, by taking into account the spectral shape of the typical GRB signal. In particular we used the mean values of the ratios of the fluences of the BATSE catalog in the different channels to normalize the simulated channels. We then added Poissonian noise to the synthetic light curves in order to simulate the background. What we obtained are plausible light curves. Each individual structure has a “Fast Rise Exponential Decay” (FRED) shape, in agreement with the observations. The count rates are comparable with those measured by BATSE, and the power density spectrum is a power law, in agreement with Beloborodov et al. (bel00 (2000)). In Fig. 6 are shown two typical light curves.

In this model there are two time scales that determine the properties of the light curves. One is the duration of a single spike, which depends on the rise time and on the decay time. The second is the characteristic waiting time between pulses. These two characteristic times reflect the typical time scales at the sources for various physical models. For example, in the internal shocks scenario the rise and the decay times of the pulse are determined by the hydrodynamic time scale and by the angular spreading time scale (Sari & Piran sari97 (1997), Nakar & Piran naka02 (2002)). The separation between pulses reflects the activity of the “central engine” (Kobayashi et al. kob97 (1997)). In the external shock scenario, the variability of the light curve arises from the characteristic dimensions of the circumburst medium (Dermer derm99 (1999)). The duration of the peaks and the separation between consecutive pulses are related to the length scales of the irregularities of the surrounding medium.

In Fig. 7 we present the relations between the computed fluence and the scaling index $`\delta `$ for synthetic light curves obtained from different combinations of the parameters. Each curve is the result of 100 GRB, with fixed number of peaks $`N_\text{p}`$, fixed peak duration $`t_\text{d}`$ and fixed mean waiting time between peaks $`\tau `$, and varying the amplitude of the signal (by multiplying the light curve by a scale factor). The relation found analyzing the BATSE catalog, and plotted with the dashed dotted line, is the reference behavior for the “intensity-complexity” relation.

From this analysis we can see that if the characteristic time between pulses ($`\tau `$) is longer than the duration of the pulses ($`t_\text{d}`$), the growth of the scaling index $`\delta `$ is too fast with respect to the observed logarithmic relation (model 1 in Fig. 7). On the other hand, the solutions with an expectation time between the pulses shorter than the duration of the pulses follow better the observed relation. This suggests that the preferred structure of the temporal series of GRB is composed by overlapped spikes while well-separated spikes make the $`\delta `$ index growth too fast. Models 2, 3, and 4 follow, at high fluences, the observed relation. For these models the spikes are overlapping and the temporal duration of the structures is proportional to the number of peaks. Under the assumption that $`t_\text{d}\tau `$ the different choices of the parameters give similar types of behavior. Model 5 shows that, if the number of the elementary spikes increases, the $`\delta `$ index grows slowly, and reaches the upper limit at higher fluences. At low fluences the growth of the index $`\delta `$ for the simulated light curves is rapid compared to the observed relation. This is also because the DE method is very sensible to low signal-to-noise ratios, and, as soon as the fluence increases the light curve reveals its behavior in terms of complexity. This shows that a tuning between the number of the spikes, the duration of the spikes, and the time between spikes is necessary to obtain the observed dependence. The relation between the fluence and the scaling index $`\delta `$ is a key study to understand the mechanism of production of gamma- ray and the environmental set up for GRB.

### 3.5 Diffusion Entropy as a function of time

Our final application of the DE algorithm is the analysis of the time evolution of the $`\delta `$ indices during a single GRB. This can be done by dividing the total light curve in overlapping intervals and applying the DE method to each of them. We proceeded as follows: we divided the GRB light curve in $`M`$ overlapping windows; each of fixed size $`W`$. We applied the DE method to each new window separately. We chose overlapping windows to study the continuous time evolution of the $`\delta `$ index. We refer to this procedure as the computation of the diffusion entropy index as a function of time (the reference time is the starting point of each series).

Figure 8 shows the results of this analysis for the two selected bursts. The value of the scaling index $`\delta `$ starts from one, indicating the presence of flaring activity, and then decreases, relaxing to the value 0.5 that indicates the end of the burst. The top panel of Fig. 8 shows the behavior of the diffusion entropy index as a function of time for the GRB 910429. The activity of the burst is clearly visible in the first eighty seconds, in agreement with the recorded data of Fig. 2. After 75 seconds the DE index reaches the dashed line at $`\delta _\text{N}`$, representing the noise value. This estimate of the burst duration is in agreement with the $`T90`$ value for this burst.

The bottom panel of Fig. 8 displays the DE index time evolution for GRB 910601. The duration of the burst using the Diffusion entropy approach, is longer than the duration given by the $`T90`$ parameter. This is due to the fact that DE is sensitive to low signal-to-noise ratio, and the presence of a low signal after the main burst can keep the $`\delta `$ index high. Furthermore, the scaling index $`\delta `$ in the fourth energy channel approaches $`\delta _\text{N}`$ more rapidly than in other channels. Using the simple model developed in the previous section we tested this behavior and we found that similar results can be achieved by adding to the spiky simulated burst of Fig. 6 a synthetic peak with a long exponential decay time ($`t_\text{d}50`$) and small amplitude (comparable with the noise level). This artificial spike represents either a weak pulse, part of the burst main emission but longer, or a fluctuation of the instrumental background. For example, a saturation effect can cause a long-range memory effect. Since the fourth channel has a lower signal-to-noise ratio, the $`\delta `$ index in this channel systematically reaches the noise value ($`\delta _\text{N}`$) earlier than in the other channels, where the signal-to-noise ratio is higher. The DE method is sensitive to the balance between uncorrelated noise and non-stationary signal, and a shorter duration of the non-stationary signal in the fourth channel could be the cause of a lower $`\overline{\delta }`$. This could explain why the mean $`\delta `$ value for the GRB catalog in the fourth channel is slightly smaller than in the first three channels.

Moreover, the DE as a function of time shows that GRB 910601 seems to have a re-flaring from 90 to 140 seconds. In this time interval the $`\delta `$ index grows again even if the signal intensity remains almost constant. This behavior is easy to see in Fig. 8 (bottom) while it is almost invisible looking at the light curve (bottom panel of Fig. 2). This shows that the DE could provide a new method to detect re-flaring activity and long time correlation in GRB time history. We have to stress that the presence of re-flaring activity has been observed in several cases in both short and long bursts. This effect could, in principle, be considered when evaluating the duration of GRB, in particular in the case of short bursts when the prompt emission lasts less than two seconds but the anomalous activity seems to last longer. However, we notice that, in some cases, the increase of DE after the prompt phase is due to fluctuation of the background. Both the orbiting motion of the instrument and the presence of other flaring sources in the field of view of the instrument can give an increment of the counts, and an increment of the $`\delta `$ indices. These cases have been confirmed by the BATSE catalog 222http://gammaray.msfc.nasa.gov/batse/grb/catalog/4b/4br\_comments.html. The DE approach detects these.

## 4 Conclusions

In this paper we performed a new analysis on the GRB phenomena applying a technique devoted to study the anomalous behavior in diffusive processes. The Diffusion Entropy method provides a measure of the non-stationarity/memory in GRB time series. The scaling index $`\delta `$ gives a well-defined value that represents the degree of complexity in a GRB light curve, with a precision of the order of few percent. In other words, the scaling index $`\delta `$ estimates how far a GRB signal is from uncorrelated noise, for which $`\delta =\delta _\text{N}=0.5`$. To remove the dependence of the scaling index $`\delta `$ on the duration of the burst, and to obtain an estimation of the degree of complexity of the signal independently of its duration we selected a portion of the light curve proportional to the $`T90`$ value of the BATSE catalog. This also allows a general description of the GRB phenomena that does not depend on the particular experimental set up. The DE method can be successfully applied in the case of long bursts, while, in the case of short bursts, the analysis is more complicated due to the lack of statistics.

For long bursts it has been possible to build up statistics applying the DE method on the BATSE catalog, obtaining a new description of the GRB phenomena in terms of their scaling index $`\delta `$. The $`\delta `$ indices are distributed between the value $`\delta _\text{N}`$ and 1 with mean $`0.8`$. The mean values of the distributions of the $`\delta `$ indices are almost the same for the first three energy channels, and smaller for the fourth channel. Furthermore, using the time series obtained by summing the light curves in the four channels, we studied the dependence of the scaling index $`\delta `$ on the measured power expressed by the total fluence. We found a logarithmic growth of the complexity with respect to the fluence. This feature connects a new variable (complexity) with an observable quantity (fluence), and the relationship between those quantities can be used to test GRB models and scenarios.

Using a simple statistical model where the parameters are the number and the duration of elementary spikes and the mean separation between them, we tried to reproduce the growth of the complexity with the fluence for synthetic bursts. We found that, in order to satisfy the observed relation between the $`\delta `$ index and the fluence, the duration of the single spikes should not be shorter than the mean separation between spikes. The pulses of the synthetic light curves cannot be completely separated and the elementary structures have to be superimposed. The connection between the typical time scales of the signal and the typical dimensions at the source suggests a well defined scenario: the emitting regions are structured and the spikes of the light curve are the results of the superimposition of many elementary spikes. In the internal shock scenario, for example, the crossing time of the shells (duration of the spikes) has to be of the same order or even longer than the travel time between shocks (waiting time between pulses), which also reflects the activity of the central engine. In the external shock scenario, the variability is related to the irregularity of the surrounding environment and the sizes of the clouds have to be bigger than the separation between them.

To study the evolution of the complexity during a single Gamma-Ray Burst, we proposed a more detailed analysis by means of the application of the DE method to subset of data, obtained from the data by selecting overlapping windows. In this case the DE method operates over smaller regions, compared to the global time history of the GRB: this implies that DE is more sensible to short time correlations. We found that the time that the $`\delta `$ index takes to reach the $`\delta _\text{N}`$ value (relaxation time) is shorter in the fourth energy channel compared with the relaxation time at lower energies. This behavior is due to the different signal-to-noise ratios of the different channels. We showed that the DE index as a function of time permits the detection of weak flares (or re-flaring activities) after the prompt GRB emission. In several cases we found that the $`\delta `$ index increased after the prompt phase. These spikes are weaker, compared to the prompt phase, and sometimes comparable with the background. They are indeed invisible looking at the light curve itself. The DE method is sensible to changes of statistics in the time series even if the intensity of the signal is comparable to the background level.

In further investigations we will use the available (first 1-2 seconds) high resolution (TTE) data for studying the short time correlation in GRB signal, especially in the case of short bursts, and to study the evolution of the scaling $`\delta `$ index with higher temporal resolution. The Gamma-Ray Large Area Telescope (GLAST) will be able to collect more data and, thanks to its short dead time, it will record GRB signals from the deep universe with a temporal resolution never reached before. The light curves of the GRB will be resolved up to the millisecond time scale, and their structure will be better understood.

###### Acknowledgements.

We are very grateful to Paolo Grigolini and Steve Shore for useful discussions during the preparation of this work.