text:
# FUSE Detection of Galactic O VI Emission in the Halo above the Perseus Arm

## 1 INTRODUCTION

The interstellar medium (ISM) consists of several phases. The hottest phase reaches temperatures of a few million degrees and is observed in X-ray emission. Gas at slightly cooler temperatures ($`3\times 10^5`$ K) is best traced by O VI, which has a strong resonance doublet at 1032/1038 Å. The O VI ion is generally assumed to trace collisionally-ionized gas, because photoionization by hot stars cannot easily create such high ionization states in the diffuse ISM (e.g., Martin & Bowyer 1990). Possible scenarios for the production of O VI include the cooling of gas originally heated by supernova shocks and turbulent mixing of hot and cold gas at the surfaces of infalling high-velocity clouds (HVCs) as they pass through the Galactic corona.

To date, Galactic O VI $`\lambda `$1032 has been observed in emission along six sight lines with the Far Ultraviolet Spectroscopic Explorer (FUSE) (Dixon et al. 2001; Shelton et al. 2001; Shelton 2002a; Welsh et al. 2002). The location of the O VI emitting gas along these sight lines is generally unknown. Welsh et al. (2002) note the narrow range of the observed intensities and suggest that the O VI emission originates in the boundary of the Local Bubble. Absorption-line measurements along about 100 sight lines (Savage et al., 2002) reveal only a poor correlation of O VI with other tracers of the ISM, such as neutral and ionized hydrogen or soft X-ray (SXR) emission. Savage et al. (2002) conclude that the observed O VI absorption is best explained by a patchy, absorbing disk with a scale height of about 2.3 kpc. This thick disk corotates with the Galactic plane and shows signs of outflow along about 20 % of the observed sight lines.

We have observed O VI emission along two sight lines, only one of which allows reliable measurements using line profile fitting. We compare the physical parameters of the O VI emission lines along this sight line with gaseous features observed at other wavelengths representing different phases of the ISM to determine the location of the O VI emitting gas and to investigate possible processes leading to the observed O VI emission.

## 2 OBSERVATIONS AND DATA REDUCTION

Our data sets (program IDs S4054801 and S4056101, hereafter SL5 and SL6, respectively) were obtained with FUSE as background observations recorded in time-tag mode on 2002 January 19–22. Our data are the first observations taken after the recovery of the satellite from a reaction wheel failure in December 2001 (Ake et al., 2002). We used the data obtained through the low-resolution (LWRS) aperture ($`30\mathrm{}\times 30\mathrm{}`$) on detecter segment LiF1A (see Moos et al. (2000) and Sahnow et al. (2000) for a complete description of FUSE and its on-orbit performance). The total exposure time of the 55 exposures of SL5 was 217 ksec, with 95 ksec recorded during orbital night. The eight exposures of SL6 yielded a total exposure time of 32 ksec, with 13 ksec during orbital night. During these observations, the satellite was pointed near the orbital pole (see Table 1 for observational data).

The data were processed using the CALFUSE pipeline version 2.0.5. Only the first steps of the pipeline were executed, i.e. the data were screened for pulse height (using limits 4–15) and corrected for Doppler effects of the satellite’s motion. Bursts were also removed, but no walk or drift correction was applied. Then the data were combined using the program ttag\_combine. Nearby airglow lines were used to determine the constant height of the extraction window. The wavelength calibration was derived from the position of the airglow lines. The resulting spectra for sight line SL5 are shown in Figure 1. Lines of O VI $`\lambda \lambda `$1032,1038 and C II$`{}_{}{}^{}\lambda `$1037 are clearly seen in the SL5 data. Due to the lower signal-to-noise ratio in the SL6 spectra, we will focus on the SL5 measurements in our analysis and refer only briefly to the SL6 data in section 3.

We followed the arguments of Shelton et al. (2001) to exclude terrestrial airglow, solar O VI emission, and scattered light inside the satellite as other possible explanations for the observed emission lines. Possible contamination (1$`\sigma `$ upper limit) of O VI $`\lambda `$1032 by fluorescing H2 is less than 200 LU (LU = photons s-1 cm-2 sr-1) for SL5. An unknown feature appears in the day-plus-night spectra on the blue side of O VI $`\lambda `$1032 but is absent in the night-only data. This feature is assumed to be airglow and has been seen in other data sets (e.g., Shelton et al. 2001). Due to its proximity to the O VI $`\lambda `$1032 line, reliable measurements for O VI could be obtained only from orbital night-time data.

## 3 RESULTS

We measured the line fluxes for SL5 in two ways. In the photon-counting method, we summed the counts in the wavelength region occupied by the emission line and subtracted a constant continuum derived from the surrounding spectrum. One-sigma uncertainties were derived assuming Gaussian statistics. In the second method, we used a line fitting algorithm, which determined the minimum $`\chi ^2`$ for a range of synthetic line profiles using three free parameters for each line: central wavelength, line width (FWHM), and peak intensity. The background was fitted by a straight line. We assumed that the emission line profiles were a convolution of a Gaussian with the instrumental line spread function, which for filled LWRS slit emission is a 106 km s-1 wide flat-top profile. One-sigma uncertainties for each parameter were derived by varying the parameter and optimizing $`\chi ^2`$ with the other parameters until the newly calculated $`\chi ^2`$ had increased by 1 relative to the optimal solution (Avni, 1976).

The results of both methods are given in Table 2 for the SL5 data. $`I_{\mathrm{tot}}`$ is the integrated intensity of the emission line, FWHM is the intrinsic line width, and $`\lambda _{\mathrm{cen}}`$ is the central wavelength of the line. The derived velocity was converted to the local standard of rest ($`v_{\mathrm{LSR}}`$). The large uncertainties of $`I_{\mathrm{tot}}`$ for O VI $`\lambda `$1038 and C II$`{}_{}{}^{}\lambda `$1037 in the $`\chi ^2`$ method reflect the close proximity of these features: reducing the strength of one line was compensated by increasing the strength of the other, allowing a larger parameter space with only slow changes in $`\chi ^2`$. The extinction toward SL5 is approximately $`E_{\mathrm{B}\mathrm{V}}=0.024`$ (Schlegel, Finkbeiner, & Davis 1998), corresponding to an attenuation in the O VI intensity of $`30`$ %, if all of the reddening occurs in front of the O VI emitting gas and the extinction parameterization of Cardelli, Clayton, & Mathis (1989) is assumed. The line ratio O VI $`\lambda `$1032/O VI $`\lambda `$1038 would approach 2 in optically thin gas, whereas in an optically thick medium the line ratio would approach 1 due to self-absorption. We derived a line ratio of $`2.0\pm 0.8`$; the large uncertainty makes it impossible to constrain the optical depth of the O VI emitting gas.

Our second sight line, SL6, lies near SL5 and shows O VI $`\lambda `$1032 emission in the day-plus-night spectrum. However, the emission line is only marginally detected in the night-only data and does not allow reliable measurements. The day-plus-night spectrum suggests an O VI $`\lambda `$1032 intensity and velocity comparable to those measured in the SL5 data.

## 4 DISCUSSION

Figure 2 shows a schematic picture of the Milky Way. The positions of the four major spiral arms are adopted from Taylor & Cordes (1993). The model of HVC Complex C is based on H I contours (Wakker, 2001) assuming a roughly constant thickness and $`z`$ height for Complex C. No assumptions for the depth of Complex C are made. All O VI emission sight lines observed with FUSE are shown. While the other sight lines have Galactic latitudes of about 60$`\mathrm{°}`$ or more (except for sight line 2, hereafter SL2), SL5 has a latitude of only about 40$`\mathrm{°}`$, thus probing a different part of the Milky Way.

SL5 intersects HVC Complex C at its southern boundary. Due to the large extinction of Complex C, it is unlikely that the observed O VI emission originates in gas beyond these HVCs, i.e. the location of Complex C can be considered an upper limit on the distance to the O VI emitting gas. A (weak) lower limit on the distance to Complex C at $`l90\mathrm{°}`$ is 6 kpc (Wakker, 2001). If the source of the O VI emission were the surface of infalling clouds, we could expect O VI emission at the same high velocities as measured in the Leiden-Dwingeloo H I Survey in this direction. The H I HVC component is observed at $`v_{\mathrm{LSR}}=115`$ km s-1 (B. Wakker 2002, private communication). The SL5 data show no O VI counterpart at this velocity. The 1$`\sigma `$ upper limit for HVC O VI is 300 LU. We conclude that the observed O VI emission does not come from infalling HVCs.

The SL5 sight line passes $``$3 kpc above the Perseus arm (Fig. 2). The velocity of the Perseus arm in the direction of SL5 ($`l=95\stackrel{}{\mathrm{.}}4`$) is $`70\pm 10`$ km s-1 (Kepner, 1970). Assuming corotation of associated gas at higher altitudes, the Perseus arm velocity projected onto sight line SL5 at $`b=36\stackrel{}{\mathrm{.}}1`$ is $`58\pm 8`$ km s-1. This is consistent with the velocities observed in the SL5 spectrum and implies that the O VI emitting gas has no velocity component perpendicular to the disk. Figure 3 shows a map of H$`\alpha `$ emission from the Wisconsin H-Alpha Mapper (WHAM) Survey (courtesy of M. Haffner & G. Madsen). A group of strong H$`\alpha `$ filaments extends between $`l=85\mathrm{°}`$ and $`110\mathrm{°}`$ from the disk into the halo indicating an outflow of gas. Sight line SL5 intersects an H$`\alpha `$ filament above this region of strongest H$`\alpha `$ emission. The WHAM spectra closest to our SL5 pointing show H$`\alpha `$ emission at velocities ranging from $`50`$ to $`+20`$ km s-1 (M. Haffner 2002, private communication). The H$`\alpha `$ filaments are most pronounced in the velocity range from $`45`$ to $`25`$ km s-1 at higher latitudes and $`25`$ to 0 km s-1 closer to the disk. This suggests that the outflow velocity decreases with $`z`$ in H$`\alpha `$. A simple Galactic rotation model (a corotating halo with $`v=220`$ km s-1) places gas with the observed O VI velocity at $`z5.5`$ kpc and at a distance of $`9.3`$ kpc from the Sun along sight line SL5 in the neighborhood of the Perseus arm H$`\alpha `$ filaments. The similar velocities of these filaments, the Perseus arm, and the O VI emitting gas suggest that these three components are associated with each other: The shock-heated gas ejected by supernova explosions in the Perseus arm and observed in SXR emission (Snowden et al., 1997) is cooling down or mixing with cooler gas and emitting the O VI photons that we see. The O VI emitting gas is on the verge of falling back to the disk, while the underlying outflow supplies new, hot gas.

The measured intensities in the SL5 data ($`1600\pm 300`$ LU at 1032 Å and $`800\pm 300`$ LU at 1038 Å) are lower than all previously measured O VI emission line intensities (2000–3300 and 1100–2000 LU, respectively). All the previous sight lines point toward high Galactic latitudes inside what Welsh et al. (2002) call the Local Interstellar Chimney. Only SL2 and SL5 (both at $`|b|40\mathrm{°}`$) appear to intersect the walls of the Local Bubble and the Chimney (see Sfeir et al. 1999 for contours). A new study by Shelton (2002b) yields 2$`\sigma `$ upper limits of 420 LU (1032 Å) and 540 LU (1038 Å) for O VI emission originating in the Local Bubble, intensities much lower than any of the O VI detections. This suggests that a significant fraction of the observed O VI emission along SL2 and SL5 originates beyond the Local Bubble. The lower intensities in the SL5 data can be explained by more extinction due to the larger distance to the O VI emitting region, if it in fact lies above the Perseus arm as suggested above. Although the O VI intensities differ by a factor of about 2 between SL2 and SL5, the SXR fluxes are similar ($`930\pm 60`$ and $`944\pm 15\times 10^6`$ counts s-1 arcmin-2, respectively; Snowden et al. 1997111Values obtained using the NASA HEASARC X-Ray Background Tool at http://heasarc.gsfc.nasa.gov/cgi-bin/Tools/xraybg/xraybg.pl). It is not clear whether most of the observed SXR emission comes from the same region as the O VI emission.

The correlations among the velocities of spiral arms, H$`\alpha `$ filaments, and O VI emitting gas that we find for SL5 are not present for any of the previously-observed O VI emission sight lines. For example, the spiral arm and O VI velocities along SL2 are consistent only if the O VI emitting gas has an outflow velocity of $`100170`$ km s-1 at $`z=19`$ kpc. No H$`\alpha `$ filaments or other signs of outflows were found near the previously-observed sight lines. (Note: SL2 points toward a region not observed by the WHAM survey.) The corotating halo model would place the O VI emitting gas $`1020`$ kpc above the plane (sight lines 2, 4a, 4b) or even higher ($`>600`$ kpc and $`>70`$ kpc for sight lines 1a and 1b); sight line 3 yields positive O VI velocities, whereas the model requires negative velocities. We therefore conclude that SL5 is the only sight line for which a reasonable and consistent explanation of the O VI emission and its parameters can be found, and that SL5 and possibly SL6 probe halo gas that originated in the underlying spiral arm.

## 5 CONCLUSIONS

We have observed Galactic O VI emission along two low-latitude sight lines, only one of which yields reliable measurements. Based on the observed velocities of the O VI emitting gas, the H$`\alpha `$ filaments, and the underlying spiral arm, we conclude that the O VI emission seen along sight line SL5 (and probably SL6) originates in the halo about $`3`$ kpc above the Perseus arm. Our observations are consistent with a Galactic fountain scenario in which gas is heated and ejected by supernova explosions in the spiral arms, becomes boyant in the colder disk gas, and rises into the halo where it cools down and sooner or later falls back to the disk. We find no O VI counterpart to the HVC component observed in H I.

This research is supported by NASA contract NAS5-32985 to the Johns Hopkins University. The authors thank L. M. Haffner and G. J. Madsen for creating the WHAM map. The Wisconsin H-Alpha Mapper is funded by the National Science Foundation. The authors are also thankful to B. P. Wakker for his comments on HVC Complex C and the H I line profiles. This work made use of the NASA Extragalactic Database (NED) and the NASA Astrophysics Data System (ADS).


format: text

text:
# Using the Scaling Analysis to Characterize Financial Markets

## 1 Introduction

The scaling concept is increasingly applied outside the traditional physical sciences domain . In the recent years, its application to financial markets, initiated by Mandelbrot in the 1960 , has largely increased also in consequence of the abundance of available data . Two types of scaling behaviors are studied in the finance literature:

1. The behavior of some forms of volatility measure (variance of returns, absolute value of returns) as a function of the *time interval* on which the returns are measured. (This study will lead to the estimation of a *scaling exponent* related to the Hurst exponent.)
2. The behavior of the tails of the distribution of returns as a function of *the size of the movement* but keeping the time interval of the returns constant. (This will lead to the estimation of the *tail index* of the distribution .)

Although related, these two analysis lead to different quantities and should not be confused as it is often the case in the literature as can be seen in the papers and debate published in the November 2001 issue of Quantitative Finance. For more explanations about this and the relation between the two quantities, the reader is referred to the excellent paper by Groenendijk et al. . In this study, we are interested in the first type of analysis. Until now, most of the work has concentrated in studies of particular markets: Foreign Exchange , Stock  or Fixed Income . These studies showed that empirical scaling laws hold in all these markets and for a large range of frequencies: from few minutes to few months. In a recent book , the hypothesis of heterogeneous market agents was developed and backed by empirical evidences. In this view, the agents are essentially distinguished by the frequency at which they operate in the market. The scaling analysis, which looks at the volatility of returns measured at different time intervals, is a parsimonious way of assessing the relative impact of these heterogeneous agents on price movements. Viewing the market efficiency as the result of the interaction of these agents , brings naturally to think that it is the presence of many different agents that would characterize a mature market, while the absence of some type of agents should be a feature of less developed markets. Such a fact should then reflect in the measured scaling exponents. The study of the scaling behaviors must therefore be an ideal candidate to characterize markets. To further explore this issue, we perform an empirical analysis of daily data across different financial markets to examine the similarities or differences in the scaling properties.

Until the 1960s, the only stochastic and scaling model in finance was the Brownian motion, originally proposed by Bachelier in 1900 , and developed several decades later . This theory predicts that the returns of market prices should follow a normal distribution with stable mean and finite variance. However, there are ample empirical evidences that the returns are not normally distributed but have higher peak around the mean and fatter tails . Moreover, it is also observed that volatility clustering is a general characteristic of financial markets . Generalizations of the classical Brownian motion were made by Mandelbrot and followers involving either fractional Brownian motions , or Lévy motion . Closely related additive scaling models have also been developed : Brownian, fractional Brownian, Lévy processes. The above approaches generally involve additive monofractal processes and analyses; but, in contrast, several scaling systems appear to be more complex. Recently, a controversy has erupted  to know if the processes that describe financial data are truly scaling or simply an artifact of the data. Moreover, several publications propose new scaling models or empirical analyses that better describe empirical evidences . It should be however noted that - as underlined by Stanley et al.  \- in statistical physics, when a large number of microscopic elements interact without characteristic scale, universal macroscopic scaling laws may be obtained independently of the microscopic details.

In this paper we address the question of the scaling properties of financial time series by empirically analyzing daily data for Foreign Exchange rates, Stock Market indices and Bond futures (described in Section 2). We study very developed as well as emerging markets in order to see if the scaling properties differ between the two and if they can serve to characterize and measure the development of the market. Here the scaling law is not used to conclude anything on the theoretical process but on the contrary we use it as a “stylized fact” that any theoretical model should also reproduce. The purpose of this paper is to point out how a relatively simple statistics gives us indications on the market characteristics, very much along the lines of the review paper by Brock . In Section 3, we recall the theoretical framework and we introduce the generalized Hurst exponents analysis. In Section 4 we describe the methodology utilized to empirically analyze the data. In Section 5 we compute and compare the scaling spectral exponents and the Hurst exponents. Finally some conclusions are given in Section 6.

## 2 The Studied Markets

We study several financial markets which are at different development stage: mature and liquid markets, emerging and less liquid markets. These markets deal with different instruments: equities, foreign exchange rates, fixed income futures. In particular, the data that we analyze are: Foreign Exchange rates (FX) (see Table 1), Stock Market indices (SM) (see Table 2), Treasury rates corresponding to twelve different maturity dates (TR) (see Table 3) and Eurodollar rates having maturity dates ranging from $`3`$ months to $`4`$ years (ER) (see Table 4). Hereafter we give a brief description of the time-series studied in this paper.

* FX: The Foreign Exchange rates (Table 1) are daily rates of exchange of major currencies against the U.S. dollar. The time series that we study go from $`1990`$ to $`2001`$ and $`1993`$ to $`2001`$. These rates have been certified by the Federal Reserve Bank of New York for customs purposes. The data are noon buying rates in New York for cable transfers payable in the listed currencies. These rates are also those required by the Securities and Exchange Commission (SEC) for the integrated disclosure system for foreign private issuers. The information is based on data collected by the Federal Reserve Bank of New York from a sample of market participants.
* SM: The Stock Market indices (reported in Table 2) are 32 of the major indices of both very developed markets and emerging markets. These daily time series range from $`1990`$ or $`1993`$ to $`2001`$.
* TR: The Treasury rates (Table 3) are daily time series going from $`1990`$ to $`2001`$. The yields on Treasury securities at ‘constant maturity’ are interpolated by the U.S. Treasury from the daily yield curve. This curve, which relates the yield on a security to its time to maturity, is based on the closing market bid yields on actively traded Treasury securities in the over-the-counter market. These market yields are calculated from composites of quotations obtained by the FD Bank of New York. The constant maturity yield values are read from the yield curve at fixed maturities, currently 3 and 6 months and 1, 2, 3, 5, 7, 10, and 30 years. The Treasury bill rates are based on quotes at the official close of the U. S. Government securities market for each business day. They have maturities of 3 and 6 months and 1 year.
* ER: The Eurodollar interbank interest rates (Table 4) are bid rates with different maturity dates and they are daily data in the time period $`1990`$-$`1996`$ .

As an example, the behaviors of FX rates and the SM rates for Japan (JPY/USD and Nikkei 225) and Thailand (THB/USD and Bangkok SET) as a function of time $`t`$ are shown in Fig. 1 in the time period $`1997`$-$`2001`$. Another example is given in Fig. 2 that shows the TR time series as a function of $`t`$ at different maturities dates in the time period $`1997`$-$`2001`$ and the ER time series as a function of $`t`$ at different maturities dates in the time period $`1990`$-$`1996`$.

## 3 The theoretical framework

The scaling properties in time series have been studied in the literature by means of several techniques such as the rescaled range statistical analysis $`R/S`$  and modified $`R/S`$ analysis , multiaffine analysis, detrended fluctuation analysis (DFA) , periodogram regression (GPH method) , the ARFIMA estimation by exact maximum likelihood , the $`(m,k)`$-Zipf method , the moving-average analysis technique , the Average Wavelet Coefficient Method , and -in this respect- connection to multi-fractal/multi-affine analysis (the $`q`$ order height-height correlation) have been made . The rescaled range statistical analysis ($`R/S`$ analysis) was first introduced by Harold Edwin Hurst  to describe the long-term dependence of water levels in rivers and reservoirs. It provides a sensitive method for revealing long-run correlations in random processes. This analysis can distinguish random time series from correlated time series and gives a measure of a signal “roughness”. What mainly makes the Hurst analysis appealing is that all these information about a complex signal are contained in one parameter only: the *Hurst exponent*. The original approach of Hurst to the scaling properties of the time series is recalled in details in Appendix. One of the weaknesses of the original method is that it relies on maximum and minimum data, which makes it very sensitive to outliers. In order to study the multi-fractal features of the data we here use an alternative method to the original approach of Hurst.

### 3.1 Generalized Hurst exponent

The Hurst analysis brings to light that some statistical properties of time series $`X(t)`$ (with $`t`$=$`\nu `$, $`2\nu `$, …, $`k\nu `$, …, $`T`$) scale with the observation-period ($`T`$) and the time-resolution ($`\nu `$). Such a scaling is characterized by an exponent $`H`$ which is commonly associated with the long-term statistical dependence of the signal. A generalization of the approach proposed by Hurst should therefore be associated with the scaling behavior of statistically significant properties of the signal. To this purpose we analyze the $`q`$-order moments of the distribution of the increments  which is a good characterization of the statistical evolution of a stochastic variable $`X(t)`$. It is given by the:

$$K_q(\tau )=\frac{|X(t+\tau )X(t)|^q}{|X(t)|^q},$$
(1)
where the time-interval $`\tau `$ can vary between $`\nu `$ and $`\tau _{max}`$. (Note that, for $`q=2`$, the $`K_q(\tau )`$ is proportional to the autocorrelation function $`a(\tau )=X(t+\tau )X(t)`$.)

The generalized Hurst exponent $`H(q)`$111We use $`H`$ without parenthesis as the original Hurst exponent and $`H(q)`$ as the generalized Hurst exponent. can be defined from the scaling behavior of $`K_q(\tau )`$ which can be assumed to be given by the relation

$$K_q(\tau )\left(\frac{\tau }{\nu }\right)^{qH(q)}.$$
(2)
Within this framework, we can distinguish between two kind of processes: (i) a process where $`H(q)=H`$, constant independent of $`q`$; (ii) a process with $`H(q)`$ not constant. The first case is characteristic of uni-scaling or uni-fractal processes and its scaling behavior is determined from a unique constant $`H`$ that coincides with the Hurst exponent $`H`$. This is for instance the case for self-affine processes where $`qH(q)`$ is linear ($`H(q)=H`$) and fully determined by its index $`H`$. (Recall that, a transformation is called affine when it scales time and distance by different factors, while a behavior that reproduces itself under affine transformation is called self-affine . A time-dependent self-affine function $`X(t)`$ has fluctuations on different time scales that can be rescaled so that the original signal $`X(t)`$ is statistically equivalent to its rescaled version $`\lambda ^HX(\lambda t)`$ for any positive $`\lambda `$ , i.e. $`X(t)\lambda ^HX(\lambda t)`$. Brownian motion is self-affine by nature.) In the second case, when $`H(q)`$ depends on $`q`$, the process is commonly called multi-scaling (or multi-fractal)  and different exponents characterize the scaling of different $`q`$-moments of the distribution. In this ‘curve’ of exponents $`H(q)`$, some values of $`q`$ are associated with special features. For instance, when $`q=1`$, $`H(1)`$ describes the scaling behavior of the absolute values of the increments. The value of this exponent is expected to be closely related to the original Hurst exponent, $`H`$, that is indeed associated with the scaling of the absolute spread in the increments. The exponent at $`q=2`$, is associated with the scaling of the autocorrelation function and is related to the power spectrum . A special case is associated with the value of $`q=q^{}`$ at which $`q^{}H(q^{})=1`$. At this value of $`q`$, the moment $`K_q^{}(\tau )`$ scales linearly in $`\tau `$ . Since $`qH(q)`$ is *in general* a monotonic growing function of $`q`$, we have that all the moments $`H_q(\tau )`$ with $`q<q^{}`$ will scale slower than $`\tau `$, whereas all the moments with $`q>q^{}`$ will scale faster than $`\tau `$. The point $`q^{}`$ is therefore a threshold value. In this paper we focalize the attention on the case $`q=1`$ and $`2`$. Clearly in the uni-fractal case $`H(1)=H(2)=H(q^{})`$. Their values will be equal to 1/2 for the Brownian motion and they would be equal to $`H0.5`$ for the fractional Brownian motion. However, for more complex signals, these coefficients do not in general coincide. We thus see that the non-linearity of the empirical function $`qH(q)`$ is a solid argument against Brownian, fractional Brownian, Lévy, and fractional Lévy models, which are all additive models therefore giving for $`qH(q)`$ straight lines or portions of straight lines. The curves for $`qH(q)`$ vs. $`q`$ are reported in Fig. 3 for some of the data. One can observe that, for all these time series, $`qH(q)`$ is not linear in $`q`$ but slightly bending below the linear trend. The same behavior holds for the other data. This is a clear sign of deviation from Brownian, fractional Brownian, Lévy, and fractional Lévy models, as already seen in FX rates . (Other cases showing marked deviations from Brownian motion have been discussed elsewhere .)

### 3.2 Scaling spectral density and Hurst exponent

For financial time series, as well as for many other stochastic processes, the spectral density $`S(f)`$ is empirically found to scale with the frequency $`f`$ as a power law: $`S(f)f^\beta `$. Here we use a simple argument to show how this scaling in the frequency domain is related to the scaling in the time-domain. Indeed, it is known that the spectrum $`S(f)`$ of the signal $`X(t)`$ can be conveniently calculated from the Fourier transform of the autocorrelation function (Wiener-Khinchin theorem). On the other hand, the autocorrelation function of $`X(t)`$ is proportional to the second moment of the distribution of the increments which, from Eq. 2, is supposed to scale as $`K_2\tau ^{2H(2)}`$. But, the components of the Fourier transform of a function which behaves in the time-domain as $`\tau ^\alpha `$ are proportional to $`f^{\alpha 1}`$ in the frequency-domain. Therefore, we have that the power spectrum of a signal that scales as Eq. 2 must behave as:

$$S(f)f^{2H(2)1}.$$
(3)
Consequently, the slope $`\beta `$ of the power spectrum is related with the generalized Hurst exponent for $`q=2`$ through: $`\beta =1+2H(2)`$ . Note that Eq. 3 is obtained only assuming that the signal $`X(t)`$ has a scaling behavior in accordance to Eq. 2. Here we are not making any hypothesis on the kind of underlying mechanism that might lead to such a scaling behavior.

## 4 Methodology

Let us here recall that the theoretical framework we presented in the previous section is based on the assumption that the signal has the scaling property described in Eq. 2. Moreover, we have implicitly assumed that the scaling properties associated with a given time series stay unchanged across the observation time window $`T`$. On the other hand, it is well known that financial time series show evidences of variation of their statistical properties with time, and depend on the observation time window $`T`$. The simplest case which shows such a dependence is the presence of a linear drift ($`\eta t`$) added to a stochastic signal ($`X(t)=\stackrel{~}{X}(t)+\eta t`$) with $`\stackrel{~}{X}(t)`$ satisfying Eq. 2 and the above mentioned properties of stability within the time window. Clearly, the scaling analysis described in the previous section must be applied to the stochastic component $`\stackrel{~}{X}(t)`$ of the signal. This means that we must subtract the drift $`\eta t`$ from the signal $`X(t)`$. To this end one can evaluate $`\eta `$ from the following relation:

$$X(t+\tau )X(t)=\eta \tau .$$
(4)
Other more complex deviations from the stationary behavior might be present in the financial data that we analyze. In this context, the subtraction of the linear drift can be viewed as a first approximation. The accuracy of this approximation has been verified by varying the observation time window. We verify that the results obtained for different time window sizes are all comparable, with fluctuations within a range of $`10\%`$.

Our empirical analysis is performed on the daily time series TR, ER, FX and SM (described in Section 2) which span typically over periods between $`1000`$ and $`3000`$ days. In particular, we analyze the time series themselves for the TR and ER, whereas we compute the returns from the logarithmic price $`X(t)=\mathrm{ln}(P(t))`$ for FX and SM. Moreover, all of these variables are ‘detrended’ by eliminating the linear drift as described above (Eq. 4).

We compute the $`q`$-order moments $`K_q(\tau )`$ (defined in Eq. 1) of the ‘detrended’ signals and their logarithms with $`\tau `$ in the range between $`\nu =1`$ day and $`\tau _{max}`$ days. In order to test the robustness of our empirical approach, for each series we analyze the scaling properties varying $`\tau _{max}`$ between $`5`$ and $`19`$ days. The resulting exponents computed using these different $`\tau _{max}`$ are stable in their values within a range of $`10\%`$.

We verify that the scaling behavior given in Eq. 2 is well followed and we compute the associated generalized Hurst exponent $`H(q)`$ whose values are given in Section 5.

In order to test that our method is not biased we estimate the generalized Hurst exponents for simulated random walks produced by using three different random numbers generators. We perform $`100`$ simulations of random walks with the same number of data points as in our samples ($`991`$ and $`3118`$) and estimate the generalized Hurst exponents $`H(1)`$ and $`H(2)`$ and the power spectra exponents $`\beta `$. The results are reported in Table 5. In all the cases, $`H(1)`$ and $`H(2)`$ have values of $`0.5`$ within the errors. Only when we consider uniformly distributed random numbers in the interval ($`0`$,$`1`$) (Rand uses a lagged Fibonacci generator combined with a shift register random integer generator, based on the work of Marsaglia.) we obtain for $`H(1)`$ of $`0.47\pm 0.01`$, but also in this case $`H(2)`$ is 0.5 within the errors. This shows that our method is powerful and robust and does not suffer of bias as other methods do. On the other hand, the estimations of $`\beta `$ from the power spectrum have values around $`1.8`$ (instead of $`2`$), showing therefore that this other method is affected by a certain bias.

## 5 Results and Discussion

### 5.1 Computation of the generalized Hurst exponent

In this section we report and discuss the results for the scaling exponents $`H(q)`$ computed for $`q=1`$ and $`q=2`$. These exponents $`H(1)`$ and $`H(2)`$ for all the assets and different markets (presented in Section 2) are reported in Figs. 4 and 5 respectively. Figures 4 (a) and 5 (a) refer to the Treasury and Eurodollar rates in the time period from $`1990`$ to $`1996`$. Whereas Figures 4 (b) and 5 (b) are relative to the Stock Market indices and Foreign Exchange rates in the time period reported in Tables 1 and 2. The data points are the average values of $`H(1)`$ and $`H(2)`$ computed from a set of values corresponding to different $`\tau _{max}`$ and the error bars are their standard deviations. Let us first observe that, for fixed income instruments (Figs. 4 (a) and 5 (a)), $`H(2)`$ is close to $`0.5`$ while $`H(1)`$ is rather systematically above $`0.5`$ (with the $`3`$ months Eurodollar rate that shows a more pronounced deviation because it is directly influenced by the actions of central banks). On the other hand, when Stock markets are concerned, we find that the generalized Hurst exponents $`H(1)`$, $`H(2)`$ show remarkable differences between developed and emerging markets. In particular, the values of $`H(1)`$, plotted in Fig. 4 (b), present a differentiation across $`0.5`$ with high values of $`H(1)`$ associated with the emerging markets and low values of $`H(1)`$ associated with developed ones. Moreover, we can see from Fig. 5 (b) that the different assets can be classified into three different categories: First the ones that have an exponent $`H(2)>0.5`$ which includes all indices of the emerging markets and the BCI $`30`$ (Italy), IBEX $`35`$ (Spain) and the Hang Seng (Hong Kong). A second category concerns the data exhibiting $`H(2)0.5`$ (within the error bars). This category includes: FTSE $`100`$ (UK), $`AEX`$ (Netherlands), $`DAX`$ (Germany), Swiss Market (Switzerland), Top $`30`$ Capital (New Zealand), Telaviv $`25`$ (Israel), Seoul Composite (South Korea) and Toronto SE $`100`$ (Canada). A third category is associated with $`H(2)<0.5`$ and includes the following data: Nasdaq $`100`$ (US), $`S\&P500`$ (US), Nikkei $`225`$ (Japan), Dow Jones Industrial Average (US), CAC $`40`$ (France) and All Ordinaries (Australia). We find therefore that all the emerging markets have $`H(2)0.5`$ whereas all the well developed have $`H(2)0.5`$.

For what concerns the Foreign exchange rates, we find that they show $`H(1)>0.5`$ quite systematically which is consistent with previous results computed with high frequency data , although the values here are slightly lower. An exception with pronounced $`H(1)<0.5`$ is the HKD/USD (Hong Kong) (Fig. 4 (b)). This FX rate is or has been at one point pegged to the USD, that is why its exponent differs from the others. Whereas in the class $`H(1)0.5`$ we have: $`ITL/USD`$ (Italy), $`PHP/USD`$ (Philippines), $`AUD/USD`$ (Australia), NZD/USD (New Zealand), ILS/USD (Israel), $`CAD/USD`$ (Canada), $`SGD/USD`$ (Singapore), $`NLG/USD`$ (Netherlands) and $`JPY/USD`$ (Japan). On the other hand the values of $`H(2)`$ (Fig. 5 (b)) show a much larger tendency to be $`<0.5`$ with some stronger deviations such as: $`HKD/USD`$ (Hong Kong), $`PHP/USD`$ (Philippines), $`KRW/USD`$ (South Korea), $`PEN/USD`$ (Peru) and $`TRL/USD`$ (Turkey). Whereas values of $`H(2)>0.5`$ are found in: $`GBP/USD`$ (United Kingdom), $`PESO/USD`$ (Mexico), $`INR/USD`$ (India), $`IDR/USD`$ (Indonesia), $`TWD/USD`$ (Taiwan) and $`BRA/USD`$ (Brazil).

Let us remind that $`H(2)>0.5`$ is commonly associated with a persistent behavior in the fluctuations of the returns, whereas an exponent $`H(2)<0.5`$ indicates anti-persistence .

These analysis has been done also on different time periods and the values are reported in Table 6 for the exponents $`H(1)`$ and $`H(2)`$ for the time period from $`1997`$ to $`2001`$ for Foreign Exchange rates, Stock Market indices and Treasury rates. In order to verify the stability of these results over different time periods we calculate and compare, for the Stock market data, the generalized Hurst exponents $`H(1)`$ and $`H(2)`$ in the whole time period (shown in Table 2) and in time periods of $`250`$ days. Moreover, we tested the numerical robustness of our results by using the Jackknife method  which consists in taking out randomly $`1/10`$ of the sample and iterates the procedure $`10`$ times (every time taking out data which were not taken out in the previous runs). We observe (see Fig. 6) that the generalized Hurst exponents computed on these Jackknife-reduced time series are very close to those computed on the entire series with deviations inside the errors estimated by varying $`\tau _{max}`$ (as described in Section 4). On the other hand, the analysis on sub-periods of $`250`$ days shows fluctuations that are larger than the previous estimated errors (and larger than the variations with the Jackknife method) indicating therefore that there are physically-significant changes in the market behaviors over different time periods (Fig. 6 (a)). This phenomenon was also detected in  when studying Exchange rates that were part of the European Monetary System. It seems that $`H(1)`$ is particularly sensitive to institutional changes on the market. This study confirms it for Stock indices. The scaling exponents cannot be assumed to be constant over time if a market is experiencing major institutional changes. Nevertheless, well developed markets have values of $`H(2)`$ that are on average smaller than the emerging ones. Moreover, the weakest markets have oscillation bands that stay above $`0.5`$ whereas the strongest have oscillation bands that contain $`0.5`$.

### 5.2 Spectral analysis

In order to empirically investigate the statistical properties of the time series in the frequency domain we perform a spectral analysis computing the power spectral density (PSD)  by using the periodogram approach, that is currently one of the most popular and computationally efficient PSD estimator. This is a sensitive way to estimate the limits of the scaling regime of the data increments. The results for some SM, FX, TR and ER data in the time periods $`1997`$-$`2001`$ and $`1990`$-$`1996`$, are shown in Figs. 789. For SM and FX we compute the power spectra of the logarithm of these time series. As one can see the power spectra show clear power law behaviors: $`S(f)f^\beta `$. This behavior holds for all the other data.

The non-stationary features have been investigated by varying the window-size on which the spectrum is calculated from $`100`$ days to up to the entire size of the time series. The power spectra coefficients are calculated through a mean square regression in log-log scale. The values reported in Fig. 10 are the average of the evaluated $`\beta `$ over different windows and the error bars are their standard deviations. Fig. 10 (a) refers to a time period between $`1990`$ to $`1996`$ whereas the Stock Market indices and Foreign Exchange rates (Fig. 10 (b)) are analyzed over the time periods reported in Tables 1 and 2. Moreover, the averaged $`\beta `$ values in a different time period, namely from $`1997`$ to $`2001`$ are reported in Table 7 for Foreign Exchange rates, Stock Market indices and Treasury rates. These values differ from the spectral density exponent expected for a pure Brownian motion ($`\beta =2`$) . However, we have shown in Section 4 that this method is biased and we have indeed found power spectra exponents around $`1.8`$ for random walks using three different random numbers generators.

It must be noted that, the power spectrum is only a second order statistic and its slope is not enough to validate a particular scaling model: it gives only partial information about the statistics of the process.

### 5.3 Comparison between the generalized Hurst exponent and the power spectra

We here compare the behavior of the power spectra $`S(f)`$ with the function $`f^{2H(2)1}`$ which - according to Eq. 3 \- is the scaling behavior expected in the frequencies domain for a time series which scales in time with a generalized Hurst exponent $`H(2)`$. We performed such a comparison for all the financial data and we report in Figs. 78 those for Foreign Exchange rates and Stock Market indices for Thailand and JAPAN (in the time period $`1997`$-$`2001`$). In Fig. 9 are reported the comparison for the Treasury and the Eurodollar rates having maturity dates $`\theta =10`$ years and $`\theta =1`$ year respectively. As one can see the agreement between the power spectra behavior and the prediction from the generalized Hurst analysis is very satisfactory. This result holds also for all the other data. Note that the values of $`2H(2)+1`$ do not in general coincide with the values for the power spectral exponents evaluated by means of the mean square regression. The method through the generalized Hurst exponent appears to be more powerful in catching the scaling behaviour even in the frequency domain.

## 6 Conclusion

Scaling behaviors are rather universal across financial markets. By analyzing the scaling properties of the $`q`$-order moments (Eq. 1) we show that the generalized Hurst exponent $`H(q)`$ (Eq. 2) is a powerful instrument to characterize and differentiate the structure of such a scaling properties. We show that $`qH(q)`$ has a non-linear dependence on $`q`$ which is a clear signature of deviations from pure Brownian motion and other additive models. The empirical analysis across a wide variety of stock markets shows that the exponent $`H(2)`$ is sensitive to the degree of development of the market. On one end we find: Nasdaq $`100`$ (US), $`S\&P500`$ (US), Nikkei $`225`$ (Japan), Dow Jones Industrial Average (US), CAC $`40`$ (France) and All Ordinaries (Australia); all with $`H(2)<0.5`$. Whereas, on the opposite side, we find the Russian AK&M, the Indonesian JSXC, the Peruvian LSEG, etc. (Fig. 5 (b)); all with $`H(2)>0.5`$. We observe emerging structures also in the scaling behaviors of interest rates and exchange rates. The robustness of the present empirical approach is tested in several ways: by varying the maximum time-step ($`\tau _{max}`$); by using the Jackknife method; by varying the time-window sizes; by comparing with three distinct simulated Brownian motions. We verify that the observed differentiation among different degrees of market development is clearly emerging well above the numerical fluctuations. Finally, from the comparison between the empirical power spectra and the prediction from the scaling analysis (Eq. 3, Figs. 78 and 9) we show that the method through the generalized Hurst exponent describes well the scaling behavior even in the frequency domain.

## Acknowledgments

T. Di Matteo wishes to thank Sandro Pace for fruitful discussions and support. M. Dacorogna benefited from discussions with the participants to the CeNDEF workshop in Leiden, June 2002.

## Appendix:The Hurst exponent

Let us consider a time series $`X(t)`$ defined at discrete time intervals $`t=\nu ,2\nu ,3\nu ,\mathrm{}k\nu `$. Let us define the average over a period $`T`$ (which must be an entire multiple of $`\nu `$) as

$$X_T=\frac{\nu }{T}\underset{k=1}{\overset{T/\nu }{}}X(k\nu ).$$
(5)
The difference between the maximum and the minimum values of $`X(t)`$ in the interval $`[\nu ,T]`$ is called the range $`R`$, which is defined as:

$$R(T)=max[X(t)]_{\nu tT}min[X(t)]_{\nu tT}.$$
(6)
The Hurst exponent $`H`$ is defined from the scaling property of the ratio:

$$\frac{R(T)}{S(T)}\left(\frac{T}{\nu }\right)^H,$$
(7)
where $`S(T)`$ is the standard deviation:

$$S(T)=\sqrt{\frac{\nu }{T}\underset{k=1}{\overset{T/\nu }{}}[X(k\nu )X_T]^2}.$$
(8)
The Hurst exponent is sensitive to the long-range statistical dependence in the signal. It was proved by Hurst  and Feller  that the asymptotic behavior for any *independent random process* (Poisson process) with finite variance is given by:

$$\frac{R(T)}{S(T)}=\left(\frac{\pi }{2\nu }T\right)^{1/2}$$
(9)
which implies $`H=1/2`$. However, many processes in nature are not independent random processes, but on the contrary show significant long-term correlations. In this case the asymptotic scaling law is modified and $`R/S`$ is asymptotically given by the power law behavior in Eq. 7 with $`H0.5`$. It must be noted that the lack of robustness of the original Hurst $`R/S`$ approach in the presence of short memory, heterocedasticity, multiple scale behaviors has been largely discussed in the literature (see for instance ) and therefore several alternative approaches have been proposed. Also the fact that the range relies on maxima and minima makes the method error prone because any outlier present in the data would have a strong influence on the range.


text:
# 1 Universality prediction for the screening masses

## 1 Universality prediction for the screening masses

QCD predicts the existence of a deconfinement temperature above which the fundamental degrees of freedom cease to be confined in color singlets and can propagate in isolation. This prediction is being experimentally tested in experiments where heavy nuclei are collided at high energies to produce such deconfined state of matter, the quark-gluon plasma. It is therefore imperative from the theoretical point of view to develop a better understanding of the behavior of gauge theories at temperatures above the deconfinement transition (for a recent review see e.g. Ref. ).

Since a non-zero temperature for a quantum field theory is equivalent to a finite size and periodic boundary conditions in Euclidean time, in many important respects gauge theories at high temperature can be treated as effectively living in one dimension less, a phenomenon generally referred to as dimensional reduction. More precisely, dimensional reduction can be defined through a specific pattern of pairwise degeneracies in the spectrum, characteristic of theories living in one less dimension than the original one, so that the occurrence of dimensional reduction can be quantitatively verified e.g. by Monte Carlo simulations. For $`SU(2)`$ and $`SU(3)`$ gauge theories in $`3+1`$ dimensions, this comparison was successfully performed in , where convincing evidence was provided that dimensional reduction occurs for temperatures all the way down to the deconfinement one.

Much less is known about the actual form of the dimensionally reduced theory, especially for temperatures not too far from $`T_c`$, where perturbative methods are not reliable. The simplest hypothesis, that the effective theory is the pure gauge theory with the same gauge group in one less dimension and in the confined phase, is ruled out by numerical data on the spectrum of screening masses for $`SU(2)`$ , while for $`SU(3)`$ the agreement is better. A more sophisticated proposal was made in , where the effective 3D theory is taken to be a three-dimensional gauge theory coupled to an adjoint scalar. This effective theory has been shown to correctly reproduce the low-lying states of the spectrum for both $`SU(2)`$  and $`SU(3)`$  for temperatures down to $`2T_c`$.

On the other hand, in the temperature range just above $`T_c`$, we have the possibility of making a solid conjecture about the form of the dimensionally reduced theory, at least for the case of $`SU(2)`$ or, in general, for gauge groups such that the transition is of second order: this is the scaling region above the transition, where universality arguments hold and the effective theory can be conjectured to belong to the universality class of theories whose global symmetry group coincides with the center of the gauge group. For $`SU(2)`$, the relevant universality class is the one of the 3D Ising model. This application of universality arguments to the deconfinement transition was introduced by Svetitsky and Yaffe  and many of its predictions have been accurately verified: the values of the critical indices (the most recent results can be found in Refs. ), universal amplitude ratios , and correlation functions at criticality . In this approach, the fundamental degree of freedom of the effective theory is the Polyakov loop, which is the order parameter of the deconfinement transition.

It is far from obvious a priori that the higher masses in the spectrum should be universal, since only the lowest one contributes to the free energy. However, it has recently been shown that for the three dimensional Ising universality class, in the broken symmetry phase, there is strong evidence for such universality: the Ising model and lattice regularized $`\varphi ^4`$ theory both exhibit a rich spectrum of massive excitations, and the mass ratios coincide in the scaling region . Particularly interesting is the existence of an excited mass in the scalar sector at less than twice the mass of the fundamental excitation: being the interaction attractive in the broken symmetry phase, it is natural to interpret such state as a bound state of two fundamental particles, an interpretation that was confirmed by studying the Bethe-Salpeter equation for 3D $`\varphi ^4`$ theory in the continuum .

It is therefore natural to ask whether the same spectrum characterizes $`SU(2)`$ gauge theory in the scaling region above the deconfinement temperature. This paper is devoted to a numerical verification of this ansatz. In particular we will focus on the scalar sector: we will show that a bound state of the fundamental screening mass can indeed be extracted from the correlation functions of suitably defined operators, that are constructed from Polyakov loops in a manner similar to that employed in  to study the spectrum of the Ising model. We also evaluated the lowest screening mass with angular momentum 2 in the scaling region and found a value compatible with the one of the Ising model. Preliminary results were reported in Ref. .

In the following section we describe in detail the method we used and the numerical results, while Section 3 is devoted to the discussion of the results and some speculation about the case of $`SU(3)`$.

## 2 Numerical determination of the spectrum of scalar screening masses

As a first step we considered the connected wall-wall correlation function of the Polyakov loop on a lattice with $`N_x\times N_y\times N_z\times N_t`$ sites and lattice spacing $`a`$

$$G(|z_1z_2|)\overline{P}(z_1)\overline{P}(z_2)\overline{P}(z_1)\overline{P}(z_2),$$

where

$$\overline{P}(z)\frac{1}{N_xN_y}\underset{n_1=1}{\overset{N_x}{}}\underset{n_2=1}{\overset{N_y}{}}P(n_1a,n_2a,z)$$

is the average of the Polyakov loop $`P(x,y,z)\text{Tr}_{n_4=1}^{N_t}U_4(x,y,z,n_4a)`$ over the $`xy`$-plane at a given $`z`$. The wall average implies the projection at zero momentum in the $`xy`$-plane. The correlation function $`G(|z_1z_2|)`$ takes contribution from all the screening masses in the $`0^+`$ channel, i.e. from the lowest (fundamental) mass $`m_1`$, from the first excited mass $`m_2`$, and so on. For a periodic lattice, we have

$`G(|z_1z_2|)`$ $`=`$ $`C_0e^{m_1L_z}`$
$`+`$ $`C_1\left[e^{m_1|z_1z_2|}+e^{m_1(L_z|z_1z_2|)}\right]`$
$`+`$ $`C_2\left[e^{m_2|z_1z_2|}+e^{m_2(L_z|z_1z_2|)}\right]+\mathrm{},`$
with $`L_z=aN_z`$. Here the dots represent the contribution from higher mass excitations and from multi-particle cuts. It should be noticed that, since $`\overline{P}(z)0`$ in the high temperature phase, a $`z`$-independent exponential term appears in the r.h.s. of Eq. (2. It is easy to verify that such term is sub-dominant with respect to the contribution of the fundamental mass and of the excited masses $`m_i`$, $`i>1`$ for which $`m_i<2m_1`$.

We discretized the theory on the lattice with the Wilson action and generated Monte Carlo configurations by the overheat-bath updating algorithm  with Kennedy-Pendleton improvement . Measurements were taken every 10 upgradings and error analysis was performed by the jackknife method applied at different blocking levels.

We made simulations at two values of $`\beta `$ slightly above the critical value $`\beta _c=2.29895(10)`$ determined in Ref. : at $`\beta =2.33`$ (corresponding to $`T/T_c1.10`$) we collected 1M measurements on a $`18^2\times 36\times 4`$ lattice and at $`\beta =2.36`$ (corresponding to $`T/T_c1.21`$) we collected 1.5M measurements on a $`18^3\times 4`$ lattice.

To make sure that our Monte Carlo ensembles were not contaminated by configurations with mixed phase, we plotted histograms of the Polyakov loop averaged over the lattice configuration. In all cases, histograms showed two non-overlapping peaks located on opposite sides with respect to the origin, thus indicating that no configurations with mixed phase were present in the thermal equilibrium ensemble.

We determined the effective mass through

$$G(z)=C_0e^{m_{\text{eff}}L_z}+C\left[e^{m_{\text{eff}}z}+e^{m_{\text{eff}}(L_zz)}\right]$$
(2)
and found for both $`\beta `$ values that data for $`m_{\text{eff}}`$ as a function of $`z`$ reach a plateau value, corresponding to the fundamental mass $`m_1`$, as $`z`$ is increased towards $`N_z/2`$. The deviations from the plateau value at small $`z`$ can be attributed to lattice artifacts and to the possible effect of other physical states. In order to single out the latter contribution, we rescaled the values of $`m_{\text{eff}}`$ and $`z`$ by $`m_11/\xi `$ and $`\xi `$, respectively, and put together on the same plot data from $`\beta =2.33`$ and $`\beta =2.36`$. As shown in Fig. 1, data from the two different $`\beta `$’s fall on the same curve before reaching the plateau, indicating that the pre-asymptotic behavior of the effective mass is not a lattice artifact, but it is instead a physical effect, due to higher mass excitations.

Universality arguments hold in the scaling region of the deconfinement transition: in order to check whether our choice of $`\beta `$ is in fact within this region, one can compare the values of the lowest screening mass $`m_1`$ at $`\beta _1=2.33`$ and $`\beta _2=2.36`$, whose ratio is predicted by universality to scale as:

$$\frac{am_1(\beta _1)}{am_1(\beta _2)}\left(\frac{\beta _1\beta _c}{\beta _2\beta _c}\right)^\nu $$
(3)
where $`\nu 0.63`$ is the correlation length critical index for the three dimensional Ising model. The values we measured111The value quoted for $`\beta _1`$, which is the point where the analysis of the whole spectrum was performed, was obtained with the variational method described below. The value at $`\beta _2`$ was obtained simply from the long-distance behavior of $`m_{\text{eff}}`$ as defined in Eq. (2). The value obtained with the latter method at $`\beta _1`$ is compatible with the one quoted. give

$$\frac{am_1(\beta _1)}{am_1(\beta _2)}=\frac{0.3654(31)}{0.5431(61)}=0.673(13)$$
(4)
to be compared to

$$\left(\frac{\beta _1\beta _c}{\beta _2\beta _c}\right)^\nu 0.65$$
(5)
We conclude that the scaling region includes, within the accuracy of our data, these two values of $`\beta `$, and hence that such region extends at least up to temperatures $`1.2T_c`$. On the other hand, the value of the mass at $`T=2T_c`$ reported in Ref.  is rather far from the scaling prediction, therefore the scaling region does not extend up to $`2T_c`$.

In order to detect possible excited states, we adopted the so-called “variational” method . The method consists of the following steps: (i) defining a basis of (wall-averaged) operators $`\{O_\alpha \}`$, (ii) computing the connected cross-correlation matrix among these operators,

$$C_{\alpha \beta }(z)=O_\alpha (z)O_\beta (0)O_\alpha (z)O_\beta (0),$$
(6)
(iii) diagonalizing $`C_{\alpha \beta }(z)`$ to obtain the eigenvalues $`\lambda _i(z)`$ and finally (iv) extracting the masses $`m_i`$ through the relation

$$\lambda _i(z)=c_0+c\left[e^{m_iz}+e^{m_i(L_zz)}\right].$$
(7)
The fundamental mass $`m_1`$ corresponds to the leading eigenvalue $`\lambda _1`$, the first excited mass $`m_2`$ corresponds to the next-to-leading eigenvalue $`\lambda _2`$, and so on. In practical simulations, for each value of $`z`$ we have a numerical determination of an “effective” mass $`m_i(z)`$ through Eq. (7). As $`z`$ increases, $`m_i(z)`$ reaches a plateau value which we identify with the screening mass $`m_i`$. In our simulations, the signal-to-noise ratio allowed to extract a screening mass only from the two leading eigenvalues.

The effectiveness of this method relies on the choice of a “good” set of operators; it is convenient to define operators living on different length scales (for instance, by use of recursive definitions). Moreover, by defining operators with non-trivial transformation under spatial rotation, it is possible to look for states with non-zero angular momentum.

As a first set of operators with $`0^+`$ quantum numbers, we considered those which can be built by adapting to the present case the recursive “smoothing” procedure used in Ref.  for the 3D Ising model and the 3D lattice regularized $`\varphi ^4`$ theory:

$$P^{(0)}(x,y,z)=P(x,y,z),P^{(n+1)}(x,y,z)=\text{sign}(u)\left[(1w)|u|+wy\right],$$
(8)
where $`u`$ is the average of four nearest values of the Polyakov loop on the $`xy`$-plane,

$$u=\frac{1}{4}\left(P^{(n)}(xa,y,z)+P^{(n)}(x,ya,z)+P^{(n)}(x+a,y,z)+P^{(n)}(x,y+a,z)\right),$$

and we chose $`w=0.1`$ and $`y=P(x,y,z)`$. We considered a subset of six such (wall-averaged) operators corresponding to $`n=0`$, 3, 6, 9, 12 and 15 smoothing steps. At $`\beta =2.33`$ ($`T/T_c1.10`$) on a $`18^2\times 36\times 4`$ lattice, we collected $``$1.6M measurements.

The relevant observable to check our ansatz is the ratio between the first excited and the fundamental screening masses in the scalar channel, which in the 3D Ising and $`\varphi ^4`$ theories assumes the value

$$\frac{m_2}{m_1}=1.83(3).$$
(9)
The effective values of this ratio for $`SU(2)`$ gauge theory are shown in Fig. 2: while the determination of this ratio is much more difficult than in the Ising model, the compatibility of the result with the expectations from universality is satisfactory. We have also observed that the excited scalar state cannot be detected when only the first four operators of our set (0, 3, 6 and 9 smoothing levels) are included in the analysis. This fact suggests that the first excited state in the $`0^+`$ channel couples with operators living at rather large length scales.

Finally, in order to look for states with non-zero angular momentum, we considered a second set of $`2^+`$ operators, inspired by Ref. ,

$`\overline{P}^{(n)}(z)`$ $`=`$ $`{\displaystyle \frac{1}{N_xN_y}}{\displaystyle \underset{n_1=1}{\overset{N_x}{}}}{\displaystyle \underset{n_2=1}{\overset{N_y}{}}}P(n_1a,n_2a,z)`$ (10)
$`\times `$ $`\left[P(n_1a+na,n_2a,z)P(n_1a,n_2a+na,z)\right],`$
with $`n=1,\mathrm{},5`$. By the same procedure adopted in the case of the screening mass in the $`0^+`$ channel, we got for the fundamental mass in the $`2^+`$ channel the value $`m_1^{(2^+)}=0.86(21)a^1`$. Again, we studied the observable $`m_1^{(2^+)}/m_1^{(0^+)}`$, i.e. the ratio between the fundamental mass in the $`2^+`$ channel, determined from this second basis of operators, and the fundamental mass in the $`0^+`$ channel, determined from the first basis of operators. The result of our determination, as a function of $`z`$, is shown in Fig. (3). Despite the noisiness, one can again see that the result is compatible with the value found  in the 3D Ising universality class222This mass ratio was actually evaluated in Ref.  for the 3D Ising gauge model: however duality implies  that in the broken symmetry phase the spectrum coincides with the one of the 3D Ising model.

$$\frac{m_1^{(2^+)}}{m_1^{(0^+)}}=2.56(4).$$
(11)
## 3 Discussion

The numerical results we have obtained are in agreement with our starting hypothesis: the spectrum of screening masses for $`SU(2)`$ gauge theory near $`T_c`$ is indeed the one predicted by universality arguments, i.e. it coincides with the spectrum of the transfer matrix of the 3D Ising model in the scaling region of the broken-symmetry phase. In particular a bound state of two elementary screening masses can be detected. An unambiguous verification of the coincidence of the two spectra will require both the analysis of more states and a substantial increase in the accuracy in the masses we have determined.

At first sight, application of these methods to the case of $`SU(3)`$ does not seem possible, given that the deconfinement transition is first order. However, the transition is rather weak, both in the gauge theory and in the 3-state Potts model that would be in the same universality class if universality arguments applied. Therefore the three-dimensional correlation length at criticality, while not diverging, assumes rather high values: hence, there exists a range of temperatures in which the correlation length is large in terms of lattice spacing, and it is not unreasonable to expect the spectra of the two theories to be similar in such region. We plan to investigate this issue studying the spectrum of the three-dimensional 3-state Potts model in the broken symmetry phase, and comparing it to the spectrum of $`SU(3)`$ pure gauge theory near $`T_c`$.

Acknowledgement. We are grateful to M. Caselle for many useful discussions.


text:
# Baryon Loading of Gamma Ray Bursts by Pick-up Neutrons

## 1 Introduction

An outstanding question in GRB fireballs is the fraction of baryons within them. It is suspected that this fraction is low, because they seem to expand with ultra-relativistic Lorentz factors, $`\mathrm{\Gamma }10^2`$. While the striking paucity of baryons could be accounted for by invoking energy release on field lines where baryons are confined, say by an event horizon or strong binding to strange quark matter, the question would then arise as to whether the mechanism that enforces baryon purity would be so effective that there would be none whatsoever in the fireball.

This is possible; the fireball could consist of just pairs, gamma rays, and low frequency Poynting flux, but then another question would arise: How do the pairs survive recombination while expanding from an extremely compact region? Were the fireball adiabatically expanding, baryon-free and thermal (Paczynski 1986, Goodman 1986), the $`e^+e^{}`$ pairs would mostly annihilate at an internal temperature of about 15 KeV, corresponding to a radius of not much more than $`10^9`$ cm. By contrast, proton-electron pairs would face no such problem, but would raise the first question: If there are so few, why are there any at all? \[Note that the afterglow from the giant Aug. 27, 1998 flare (Frail, Kulkarni, and Bloom 1999) from SGR 1900+14 was sufficiently weak that the outburst seems to have not put much of its energy into escaping pairs (Eichler 2002), in marked contrast to long, cosmologically distant GRB’s. This serves as a reminder that a material fireball that mostly survives the compact regions of its origin is not to be taken for granted.\]

A theory of baryon content in GRB fireballs could answer these questions. In this letter we attempt to estimate the baryon content that would arise if the fireball originated entirely on magnetic field lines that connect to an event horizon (or anything equally effective in enforcing total baryon purity). We assume that the origin of the baryons is leakage in the form of neutrons crossing from baryon rich field lines to the baryon free outflow (Eichler and Levinson, 1999, hereafter EL). The freshly picked-up neutrons are extremely energetic in the frame of the outflow, and this creates the possibility of creating very energetic ($`10^9E10^{15}`$eV) pions and neutrinos. The analysis of neutron pick-up and neutrino production rates and energies is somewhat similar to that of EL, but here we note a particular instability that exists in this situation: a collisional avalanche. This leads to particularly efficient neutron pick-up and neutrino production.

Neutrons and neutrinos have also been discussed by other authors. Derishev et al. (1999a,b) argued that roughly equal numbers of neutrons and protons should be present in the fireball, and studied the consequences for the hydrodynamics and emission during the prompt and afterglow phases (see also Beloborodov 2002). They proposed that decoupling of the neutrons during the acceleration of the fireball may lead to a burst of $`10`$ GeV neutrinos through inelastic scattering (see also Bahcall and Meszaros 2000). In all the work cited above, the mechanism of baryon contamination is not addressed. The authors merely assume that the baryon content of the fireball is fixed by some mechanism on the smallest scales, prior to its acceleration. Lemoine (2002) reexamined the conditions under which the neutrons in the fireball may remain free and concluded that fusion of neutrons and protons to 4He should precede decoupling if the dynamical timescale exceeds $`10^4`$ s (cf. Beloborodov 2002). He then argued that this would lead to a dramatic suppression of the density of free neutrons and, consequently, neutrino emission. However, he did not discuss the possibility of viscous heating at the interface between the fireball and its surroundings. In what follows we consider this possibility and argue that baryonic contamination can occur on much larger scales ($`10^{10}`$ cm) via diffusion of neutrons from a baryon rich wind into the baryon free fireball.

Additional baryon loading may arise from mixing, owing to hydrodynamic instabilities at the interface between the baryon poor jet (BPJ) and the baryon rich outflow (BRW) or the enveloping material of the host star. However, the degree of mixing may depend strongly on changes in the magnetic topology, which cannot be easily assessed. Moreover, it is unlikely that it could load all but the outmost layers of the BPJ when the hydrodynamic crossing time scale for the dense, outer material exceeds the proper expansion timescale of the BPJ.

## 2 Wind ejection

In GRB scenarios that invoke a stellar type progenitor, the ejection of a baryon poor fireball follows a catastrophic event that results in the formation of a system consisting of a compact object surrounded by a hot disk or torus. Here we assume that the fireball emanates from the immediate vicinity of the compact object (a black hole, say) essentially devoid of baryons. This degree of baryon purity is natural if the fireball is generated on field lines that thread the black hole, either by neutrino annihilation or by extraction of the black hole rotational energy (Levinson and Eichler 1993).

During the ejection of the fireball, the matter in the central parts is maintained at temperatures well above 1 MeV, at which emission of MeV neutrinos is the dominant cooling mechanism. Neutrinos escaping from the inner parts of the torus heat the surface layers and drive a powerful, baryon rich wind that propagates at a sub-relativistic speed, and confines the baryon poor fireball. Critical point analysis (van Putten & Levinson 2003, see also Levinson and Eichler 1993) yields an estimate for the baryon rich wind luminosity of

$$L_w10^{49.5}\beta _{sc}T_{t10}^4A_{13}\mathrm{erg}\mathrm{s}^1,$$
(1)
and for the corresponding mass loss rate of

$$\dot{M}=2L_w/3a_{sc}^210^{30}L_{w51}\beta _{sc}^2\mathrm{gr}\mathrm{s}^1,$$
(2)
where $`A=10^{13}A_{13}`$ cm2 is the surface area of the torus, $`T_t=10^{10}T_{t10}`$ K is the torus temperature and $`a_{sc}=\beta _{sc}c10^{10}(M/M_{\mathrm{}})^{1/2}(r_c/10^6\mathrm{cm})^{1/2}`$ cm $`s^1`$ is the sound speed at the wind critical point, with $`M`$ being the mass of the compact object, and $`r_c`$ the radius of the critical point. For a temperature of 4 MeV, a wind luminosity of $`10^{51}`$ erg s-1 is expected.

The large optical depth of the expelled wind ($`\tau _T10^8L_{w51}/r_{10}`$) and the extremely short cooling time render the wind pressure radiation dominated. The determination of the wind’s temperature and density profiles beyond the critical point requires quantitative treatment of wind acceleration in the supercritical region. Assuming for simplicity that the wind velocity $`v_w=\beta _wc`$ is constant, we obtain a wind temperature

$$T_w10^{8.5}L_{w51}^{1/4}\beta _w^{1/4}(\psi ^2\theta ^2)^{1/4}r_{10}^{1/2}\mathrm{K},$$
(3)
and baryon density

$$n_p=10^{23}L_{w51}/r_{10}^2\beta _w^3(\psi ^2\theta ^2)\mathrm{cm}^3,$$
(4)
where $`\psi `$ is the opening angle of the baryon rich wind (which, as argued below, may be collimated).

If the explosion occurs inside a star, as in the collapsar scenario for GRB’s, both the fireball and the surrounding baryon rich wind must advance first through the stellar core before they break out. During this stage they both may be collimated by the kinetic pressure of the core. As they advance through the star, the BPJ and the baryon rich wind each drives a forward shock that accelerates the stellar material. It is reasonable to assume that the BPJ and the baryon rich wind have similar isotropic equivalent luminosities in which case we expect their head velocities to be roughly the same. The temperature of the shocked wind plasma (behind the reverse shock) is given approximately by eq. (3), upon replacing $`\beta _w`$ with $`\beta _h`$. Eventually, the BPJ breaks through the host star’s envelope if there is one. Otherwise we would not observe the GRB it is supposed to produce. ¿From the point of view of neutrino production, however, the pre-breakout stage may also be important.

## 3 The neutron content of the torus wind

If the torus is a remnant of a degenerate star, then it contains nearly half its mass in neutrons. At temperatures above the dissociation temperature of the nuclei (0.5 - 0.7 MeV) the neutrons and protons become free. The neutron to proton ratio $`\xi _n`$ would then either be the initial value, or, if the weak interaction equilibration time

$$t_{weak}5\times 10^2(T/3\mathrm{M}\mathrm{e}\mathrm{V})^5s$$
(5)
is shorter than the outflow time, $`t_{exp}`$, the equilibrium value, $`\xi _n\mathrm{exp}[(m_pm_n)/T]`$, where $`m_i`$ is the mass of species $`i`$ ($`i=n,p`$). As noted by Derishev et al. (1999b), in the latter case neutrons will be effectively produced through the reaction $`p+e^{}n+\nu `$, even if the initial torus composition is predominantly hydrogen. If freeze-out of the weak interaction occurs at $`T>(m_nm_p)1`$ MeV, then $`\xi _n1`$ is generally expected.

Above the critical point, the temperature in the wind declines with radius. Once it drops below about $`T_{rec}77`$ keV (depending weakly on density) the free protons and neutrons will recombine to form deuterium. The reaction rate for deuterium formation is $`<\sigma _dv>5\times 10^{20}`$ cm3 s-1. The corresponding recombination time is then given by

$$t_{rec}10^{3.5}(\psi ^2\theta ^2)\beta _w^3r_{10}^2L_{w51}^1s,$$
(6)
where eq. (4) has been used. Recombination will be effective at radii

$$r_{10}<10^3L_{w51}\beta _w^4(\psi ^2\theta ^2)^1$$
(7)
at which $`t_{rec}<<t_{exp}`$, provided the temperature there is below $`T_{rec}`$. Given $`T_{rec}77`$ keV and equation (3)) we conclude that the radius below which free neutrons can exist in the wind is

$$r_n2\times 10^9L_{w51}^{1/2}(\psi ^2\theta ^2)^{1/2}\beta _w^{1/2}\mathrm{cm}.$$
(8)
For reasonable wind parameters, $`r_n`$ should lie in the range between a few times $`10^9`$ to a few times $`10^{11}`$ cm.

While we have assumed that the BPJ is ensheathed by a neutron-rich outflow which we expect exists, we could, for the purposes of the following discussion, assume that the BPJ is in direct contact with the envelope of the host star that collimates it. Neutrons would be freed up at the interface because the inner wall of the envelope would be heated by the BPJ.

## 4 Neutron pick-up

The free protons and neutrons in the wind are coupled by nuclear elastic scattering. At the temperatures of interest the corresponding rate is $`<\sigma _{el}v><\sigma _0c>=10^{15}`$ cm3 s-1, independent of center of mass energy.

The flux of neutrons diffusing into the BPJ from the interface separating the BPJ and the baryon rich wind is given by $`J_D(r)=\lambda _{np}v_{is}n_n/x=\lambda _{np}v_{is}(n_n/l)`$, where $`v_{is}=\beta _{is}c=(kT/m_p)^{1/2}`$ the ion thermal speed (not to be confused with the sound speed of the radiation dominated plasma), $`\lambda _{np}=\beta _{is}/(n_p\sigma _0)`$ is the mean free path for np collisions, $`x`$ denotes the cylindrical radius, and $`l`$ the gradient length scale (which may vary with $`r`$).

Assuming that the boundary between the BPJ and baryon rich wind is very sharp at the BPJ injection radius ($`r_010^7`$ cm), the gradient length scale at some larger radius $`r`$ is $`l(\lambda _{np}v_{is}t_{exp})^{1/2}`$ with $`t_{exp}=r/v_w`$ being the wind expansion time. The total number of neutrons diffusing into the BPJ below some radius $`r`$ is given by $`N_{diff}(r)2\pi \theta r^2t_{exp}J_D(r)=2\pi \theta r^2ln_n`$. Combining the above results one finds,

$$N_{diff}(r)=2\pi \theta \frac{\beta _{is}}{\beta _w^{1/2}}n_n\frac{r^{5/2}}{\sqrt{n_p\sigma _0}}10^{50}\xi _n[(\psi /\theta )^21]^{1/2}(\beta _{is}/\beta _w^2)r_{10}^{3/2}L_{w51}^{1/2},$$
(9)
where $`\xi _n=n_n/n_p`$.

If the density gradient length scale, $`l`$, at the BPJ baryon rich wind interface is larger than assumed above, e.g., as a result of some mixing, then the gradient would be smeared and the diffusive flux would be reduced, but presumably advection would replace it and keep the neutron-proton mixture hot. It is hard to see how mixing could penetrate to the BPJ center without disrupting the BPJ entirely, so any mixing would be restricted by this consideration to load only the periphery of the BPJ.

Now the fraction of neutrons drifting through the BPJ that decay, $`t_{\mathrm{cross}}/\tau _n10^{3.5}\theta \beta _{is}^1r_{10}`$, where $`t_{cross}=\theta r/v_{is}`$ is the BPJ crossing time and $`\tau _n=900`$ s is the neutron lifetime, is small. However, each decay liberates a proton that generates more protons via collisions with the undecayed neutrons. The proton fraction thus grows exponentially in what we term a collision avalanche, until becoming comparable to the neutron fraction. To estimate the growth length of the shower, we note that the density of target neutrons inside the BPJ (static in the Lab frame) is $`(t_{cross}/t_{exp})(N_{diff}/\pi \theta ^2r^3)`$, and the optical depth for a collision of a picked-up baryon with the target neutrons is

$$\tau _{np}=\sigma _{np}r(t_{cross}/t_{exp})(N_{diff}/\pi \theta ^2r^3)=10^{4.5}\beta _w^{1/2}\psi ^1r_{10}^{1/2}L_{w51}^{1/2},$$
(10)
where a cross section for inelastic pn collisions of $`\sigma _{pn}=40`$ mb has been adopted (Hagiwara et al. 2002). Evidently, the growth length of the shower is much shorter than the injection radius and it will saturate already at the base of the BPJ. At this point every neutron diffusing into the BPJ is picked up via a collision with a fast baryon coming from below.

Combining equations (8) and (9), and taking $`\beta _{is}=10^2`$ (the ion thermal speed at the recombination temperature) yields the total number of neutrons picked up by the BPJ:

$$N_{cap}=N_{diff}(r=r_n)10^{47}\xi _n\theta (\psi ^2\theta ^2)^{5/4}\beta _w^{11/4}L_{w51}^{5/4}.$$
(11)
Adopting for illustration $`\theta =\psi /2=0.1`$, $`\xi _n=1`$, $`\beta _w=0.3`$, we obtain $`N_{cap}3\times 10^{49}L_{51}^{5/4}`$.

Now the fact that the avalanche growth is so rapid shows that the inwardly drifting neutrons may be converted back to having a 50 percent proton component shortly after crossing into the BPJ, and will merely line the BPJ outer boundary with a hot viscous sublayer (just as peripheral turbulent mixing would do). The short mean free path also means that relative Lorentz factor differences across it are likely to be much less than the total difference. The important assumption is that the inner side is exposed to contact with the more relaticvistic BPJ and is kept hot enough to have a free neutron component.

The picture we are led to by our basic model assumptions is a hot viscous sublayer with density decreasing inwards, and neutrons and protons coexisting in roughly equal numbers in the most tenuous regions. Below some density, the neutrons stream freely into the interior of the BPJ. Define the free streaming density $`n_{fs}`$ to be that at which the proper hydrodynamic time $`r/c\mathrm{\Gamma }`$ equals the proper collision time $`\mathrm{\Gamma }/n<\sigma v>`$, i.e.

$$n_{fs}\mathrm{\Gamma }^2c/r<\sigma v>.$$
(12)
If the bulk Lorentz factor $`\mathrm{\Gamma }`$ is determined by the density of picked up neutrons then

$$\mathrm{\Gamma }_{fs}=L_j/(n_{fs}mc^2\pi \theta ^2r^2ch)=26r_{12}^{1/3}L_{j50}^{1/3}\theta ^{2/3}h^{1/3}$$
(13)
where $`L_j`$ is the BPJ luminosity and h is the specific enthalpy of the fluid in units of $`m_pc^2`$, and where we have used equation (12) in the last part of eq. (13). The free streaming density is

$$n_{fs}=(L_j/\pi \theta ^2m_pc^3h)^{2/3}r^{5/3}\sigma ^{1/3}.$$
(14)
The number of neutrons per unit time crossing the free streaming boundary inward within radius r is given roughly by

$$dN_{cr}/dt=\pi \theta r^2n_{fs}c/\mathrm{\Gamma }=8\times 10^{49}r_{12}^{2/3}L_{j50}^{1/3}\theta ^{1/3}h^{1/3}s^1$$
(15)
where we assume that the random component of the neutron velocity at the free streaming boundary is close to c. Thus, at $`r_{12}1`$, most of the neutrons that diffused into the BPJ at $`r_{10}1`$ are already free streaming.

When $`nn_{fs}`$, h is assumed to be close to unity. In the free streaming zone, where $`\mathrm{\Gamma }\mathrm{\Gamma }_{fs}`$, the neutrons are subjected to large shear in the BPJ, and as they move axisward, find themselves moving relative to the local frame at lab angle $`\chi `$ and nearly backward with a local Lorentz factor $`\gamma ^{}\mathrm{\Gamma }\mathrm{\Gamma }_{fs}(1\beta \beta _{fs}cos\chi )`$. Thus h may be estimated as $`\gamma ^{}\mathrm{\Gamma }\mathrm{\Gamma }_{fs}(1\beta \beta _{fs}cos\chi )`$ in the free streaming zone. At the free streaming surface, h is very geometry dependent, and can be between 1 and $`\mathrm{\Gamma }_{fs}^2`$

The local spread of neutron velocities at the free streaming boundary is of order $`1/\mathrm{\Gamma }_{fs}`$, so if $`\theta 1/\mathrm{\Gamma }`$, the transverse velocity of an ”inwardly” free streaming neutron may point away from the axis, as long as it does so less than the average local velocity, and $`\chi 1/\mathrm{\Gamma }_{fs}`$. Some fraction($`1/e`$) of the inwardly streaming neutrons will encounter further collisions, thus loading a ”collisional annulus” of thickness $`1/\mathrm{\Gamma }_{fs}`$ and radius $`\theta `$, while the rest may continue further inward until they decay. If $`\mathrm{\Gamma }_{fs}^2\mathrm{\Delta }t`$, where $`\mathrm{\Delta }t`$ is the GRB duration, exceeds the neutron lifetime $`900s\mathrm{\Gamma }_{fs}`$, then most of the neutrons decay within the fireball. Otherwise they decay behind it and leave a baryon pure core ahead.

The fluid parameters in the collisional annulus are similar to those described by Derishev and co-workers, where the radial acceleration just happens to proceed on a scale comparable to the mean collision time, inducing a modestly relativistic relative (radial) velocity between the neutrons and protons. In our picture, this apparent numerical coincidence is in fact natural for the annulus defined by marginally freely streaming neutrons. The annulus occupies $`2/\theta \mathrm{\Gamma }_{fs}0.03r_{12}^{1/3}L_{j50}^{1/3}\theta ^{1/3}h^{1/3}`$ of the beam solid angle. For $`\theta =0.1`$ and $`h^{1/3}6`$, this is comparable to the solid angle of the core, and suggests the possibility that many of the GRB we see are just these annuli. Nevertheless, the observational effects of a much higher $`\mathrm{\Gamma }`$ core are worth considering, especially because its elements can overtake the denser, peripheral neutron outflow.

It may be that the BPJ is collimated, and that its outer parts - or at least the inwardly free streaming neutrons - may eventually converge at the BPJ axis. The density can increase downstream and neutrons that have already passed through the free streaming boundary (according to the formal local definition) may then with high probability collide further downstream with impact angle $`\chi \mathrm{\Gamma }_{fs}^1`$in a region of $`\mathrm{\Gamma }\mathrm{\Gamma }_{fs}`$. (Even in the case of an asymptotically conical BPJ (Levinson and Eichler 2000), the fact that $`\mathrm{\Gamma }_{fs}`$ decreases with r means that neutrons free streaming inward from a point $`r_110^{12}`$cm, will graze those free streaming from further upstream at point $`r_210^{12}`$cm where $`r_2r_1`$. The impact angle as viewed in the lab frame is of order $`\mathrm{\Gamma }_{fs,2}^1`$.)

When $`\chi 1/\mathrm{\Gamma }_{fs}`$, freely axisward streaming neutrons are ”broadsided” by faster interior plasma at impact angle $`\chi 1/\mathrm{\Gamma }_{fs}`$, and their optical depth to baryons in the interior plasma is enhanced by a factor $`\theta ^2\mathrm{\Gamma }_{fs}^2`$. Thus, the optical depth presented by axisward streaming neutrons to the interior plasma is at least of order unity and neutrino production by picked up particles is efficient. The energy of the neutrinos released in this avalanche is of order 0.05 of the typical proton energy.

The ”top down” nature of pickup suggests that much energy can be dissipated in extremely energetic collisions. A high $`\mathrm{\Gamma }`$ flow that is slowed down by neutron pickup can be viewed as slowing down in stages, such that at each stage $`ϵ(r)dN/dtL_j`$, where $`ϵ(r)`$ is the average energy per baryon at radius r and dN/dt, also a function of r, is the rate of pickup within radius r. For each new collisional pickup, about 50 percent of the original energy is dissipated into pion decay products, mostly neutrinos. (In each collision between a moving baryon and a target neutron, many pions are produced. The leading pion has about 0.2 of the original baryon energy, and the muon neutrino from its decay will have about 0.05 of the original, provided that the impact angle exceeds $`1/\mathrm{\Gamma }_{fs}`$. Thus neutrinos of up to $`0.05ϵ`$ will emerge, and they will contain about 5 percent of the jet energy. Another 10 percent or so is in the soft pions, which have a lab frame Lorentz factor somewhat larger but of order the center of mass Lorentz factor.) In the crude approximation that the neutrino losses are a small fraction of the collision energy, the spectrum of neutrinos is then $`dN/dEE^2`$ over a dynamical range that depends on how many optical depths the avalanche proceeds through; including neutrino losses leads to an even harder spectrum. Even if neutron pickup is extremely inefficient, this would still allow efficient neutrino emission in fewer, higher energy neutrinos. In the extreme scenario where there is only neutron decay and a single optical depth for collisions, there are about $`6\times 10^{49}c\mathrm{\Delta }t/\tau _n\mathrm{\Gamma }_{fs}10^{48}/\mathrm{\Gamma }_{fs}`$ neutrons that decay within the fireball. This suggests an asymptotic $`\mathrm{\Gamma }`$ of $`10^4`$ and a maximum proton energy of $`10^8m_pc^2`$. Most of this energy will be liberated as neutrinos over an optical depth of several. Note that for BPJ proton energies of order $`10^8m_pc^2`$, even the soft pions can decay into $`E1`$ TeV neutrinos. Such extremely hard spectra would be a highly distinctive signature of the model, and the optical depth crossed by the avalanche would be manifested in the neutrino spectrum.

It is also instructive to view the problem from the point of view of a free streaming neutron. It has a better than even chance of not making another collision once entering the free streaming zone. However, the rare collision is with extremely relativistic plasma, and the energy liberated per collision goes as $`\mathrm{\Gamma }^2`$. Writing the expected energy release per path length as $`<dE/ds>=n<\sigma v>\mathrm{\Gamma }^2`$, and estimating $`\mathrm{\Gamma }^2`$ as $`L/\pi r^2\theta ^2n`$, we notice that $`<dE/ds>`$ is then independent of density, and much of the energy release can be in the form of rare but very energetic collisions.

A neutrino flux of $`10^{52}F_{52}`$ ergs could be detected at the one count level with a $`1km^2`$ neutrino detector at a distance of $`(F_{52}f)^{1/2}`$ Gpc, where f is the beaming factor. As noted earlier (Eichler and Levinson 1999, Meszaros and Waxman 2001), smothered GRB’s could also be detected in (and only in) neutrinos.

## 5 Light Element Production

Light element production by spallation was discussed by Eichler and Letaw (1987) in regard to scenarios of UHE particle acceleration within supernovae ejecta. It was estimated that each cosmic ray would produce about 10 $`{}_{}{}^{10}B`$ atoms from interaction with somewhat heavier nuclei, and that light element abundances constrain the average supernova to produce no more than about $`10^{50}`$ prompt cosmic rays in the midst of heavy nuclei. For a shock accelerated spectrum, the number density is dominated by GeV protons, and this would constrain the energy to less than $`10^{48}`$ ergs in cosmic rays. However, this would not apply to a spectrum of pickup ex-neutrals that is dominated by protons above 1 TeV. Also, an intense burst of particles on the same surface would produce ”overkill” - i.e. light elements produced by spallation would then be further degraded by repeated spallation events. This is significant at a high energy particle fluence exceeding $`10^{26}cm^2`$, though Rayleigh-Taylor instabilities at the head of the jet could lessen the overkill by constantly providing fresh surface.

Most important, however, is that a GRB-associated supernova is not typical. It may therefore be worth searching for excess light element lines in the young supernova remnants associated with GRB’s. Enhancements of two orders of magnitude above the cosmic abundance could obtain. The original GRB-associated supernova, SN1998bw, had a remarkably large outflow velocity, c/6, suggesting that the ejecta may not have been particularly massive, and may have been preferentially close to the collimated fireball within.

While energetic protons which could be channeled away from enveloping material by electromagnetic forces might never interact with the surrounding ejecta, neutrons would easily penetrate the envelope if the ejecta collimate the GRB fireball, for then the walls of the BPJ curve inward.

## 6 Conclusions

We have worked out the details of an earlier suggestion that baryon loading of GRB fireballs is accomplished by pick-up ex-neutrons that crept across magnetic field lines into the path of the collimated fireball from the collimating material. For a GRB lasting, say, 30 seconds ($`r_{12}1`$), and having conical collimation $`\theta `$, we find that there is a free streaming annulus at which the flow has a Lorentz factor $`\mathrm{\Gamma }_{fs}`$ of about $`35L_{50}^{1/4}(\theta /0.1)^{1/2}`$. (In this estimate we assumed that the specific enthalpy h, which is rapidly varying at this point, is roughly $`\mathrm{\Gamma }_{fs}`$, and note that the result is only weakly dependent on it.) The thickness of the annulus is about $`1/\mathrm{\Gamma }_{fs}`$, and $`0.6(\theta /0.1)^{1/2}`$ of the solid angle within the cone of opening angle $`\theta `$ is subtended by the annulus. Well inside the annulus, i.e. when $`\theta 1/\mathrm{\Gamma }`$, the bulk Lorentz factor may be considerably higher, and the spectrum could be considerably harder. This suggests the possibility of extremely hard GRB’s that might yield UHE photons and/or neutrinos, yet be relatively inconspicuous in soft $`\gamma `$-rays. Outside the annulus, the baryon loading is greater, $`\mathrm{\Gamma }`$ is considerably lower and this part of the outflow could be responsible for X-ray flashes (Berezinsky and Prilutsky, 1985), when the viewing angle happens to coincide with it. Prediction of a general angular distribution for observed GRB’s is frustrated by our ignorance of the intrinsic distribution for all GRB’s, as well as of the dependence of soft gamma ray efficiency on the amount of mass loading.

If the GRB is collimated by surrounding material, such as the envelope of a host star, enough that $`\theta \mathrm{\Gamma }_{fs}`$, then the transverse structure is much closer to being uniform.

The baryons in the fireball can, of course, then go on to generate further neutrinos downstream of their point of origin as in several previous discussions (Eichler 1994, Paczynski and Xu 1994, Waxman and Bahcall 1997). In this paper, on the other hand, we have presented a scenario for baryon loading whereby UHE neutrinos are a logical consequence. The neutrinos that result have individual energies of order $`\mathrm{\Gamma }^2m_\pi 10^{1215}`$ eV, which are easier to detect with large underwater and under-ice neutrino detectors than those at several GeV, which could result from differential acceleration of protons and neutrons by fireball pressure. They would have a very hard spectrum and be a highly distinctive feature of the pick-up model. Remarkably, the total energy output in neutrinos can in principle as high as that of the observed fireball, or even higher.

The ability to catch GRB-associated supernovae at an early stage should be greatly enhanced by SWIFT, and it may be possible to search for spallation-induced light element enhancement in the young supernova ejecta.

The emission of $`\gamma `$-rays from an annulus just inside an optically thick wall should give rise to a strongly polarized reflected component, as noted by EL. A quantitative discussion of this will be given in a subsequent paper.

We thank E. Derishev and the referee for useful discussions. This research was supported by the Arnow Chair of Astrophysics at Ben Gurion University and by an Adler Fellowship awarded by the Israel Science Foundation.