text:
# Radio observations of HD 93129A: The earliest O star with the highest mass loss ?

## 1 Introduction

Massive, early-type O stars at the upper end of the stellar main sequence are known to have strong stellar winds, with velocities up to $``$3000 km s-1 (Prinja et al. 1990). Their evolution into Wolf-Rayet (WR) stars is marked by high mass loss rates of $`10^7`$ to $`10^5`$ M yr-1 (e.g., Lamers & Leitherer 1993). The transition can be described by three Of phases (O – O((f)) – O(f) – Of – WR), where “f” stands for the presence of certain N iii and He ii spectral lines (see Maíz Apellániz & Walborn 2002 for details; also Table 1). Once the H-rich outer layers of the star are lost, the high temperature core is exposed, and the evolution from Of to WR star is complete. The action of the stellar wind on the ambient interstellar medium creates a cavity of hot, rarefied gas, surrounded by a slowly expanding envelope (Benaglia & Cappa 1999).

Wolf-Rayet stars have anomalously strong and broad emission lines, which are thought to come from the rapidly expanding stellar wind, but few absorption lines. With an average mass of 10 M, their mass loss rate can reach values up to 10-4 M yr-1. Van der Hucht et al. (2001, 2002) cataloged well over 200 WR stars in our Galaxy, more than half of which belong to the WN class, i.e. they show strong helium and nitrogen emission lines. A large fraction of the cataloged WR stars ($``$40%) are known to be binaries.

According to Crowther et al. (1995), the most massive ($`>`$ 60 M) O stars evolve from the Of stage directly into highly luminous late WN(+abs) and WN7 stars without passing through an intermediate Luminous Blue Variable or red supergiant phase. The slightly less massive O stars evolve through either of those intermediate phases into slightly less luminous WN8 stars (see also Langer et al. 1994). The quantitative spectroscopic properties show a smooth progression from Of stars to lower excitation or late WN stars (WNL) (Bohannan 1990). Walborn et al. (1992), Lamers & Leitherer (1993) among others, suggested that transition types Ofpe/WNL, and WNL stars, form a sequence of increasing wind density.

Examples of transition objects are Sk–670 22 (type O3If\*/WN6-A, Walborn 1982), Sk–710 34 (type O4f/WN3, Conti & Garmany 1983), and, lately, the ones found in 30 Doradus (Massey & Hunter 1998). The mass loss rates of the transition objects are expected to lie between those for Of and WR stars.

Almost 400 Galactic O stars have recently been compiled in the new ‘Galactic O Star Catalogue’ (GOS) by Maíz-Apellániz & Walborn (2002). There are only about a dozen stars with spectral types O3.5 or earlier; five of these belong to the clusters Trumpler 14 and Trumpler 16, in the Carina region. Among them one star, HD 93129A, stands out as the earliest, most luminous, and most massive star known in our Galaxy. Walborn et al. (2002) re-analyzed spectra of HD 93129A, previously classified as O3, and find that it is the prototype of a new subcategory and the sole representative of type O2 If\* (see Table 1). HD 93129A appears to have the most extreme stellar winds regarding ionization and terminal velocities (Walborn et al. 2002): $`v_{\mathrm{}}=3200\pm 200`$ km s-1, an effective temperature $`T_{\mathrm{eff}}=52,000\pm 1000`$ K, and a stellar mass of $`130\pm 15`$ M (Taresch et al. 1997). Only $`2.6\mathrm{"}`$ from HD 93129A, Walborn (1973) discovered a second star, HD 93129B, of type O3.5 V((f+)). Another star of the same type, HD 93128, is found at a projected distance of $`20\mathrm{"}`$. The stellar parameters for all three stars are listed in Table 1.

In the radio regime, free-free thermal emission from the ionized wind of early-type OB and Wolf-Rayet stars can be modeled (Wright & Barlow 1975, Panagia & Felli 1975) yielding a spectral index ($`\alpha `$) of typically $`0.60.7`$. Generally, the radio flux densities are of the order of mJy for stars at distances up to 2 kpc. About 80 WR stars (Abbott et al. 1986, Chapman et al. 1999) and $``$120 O stars (Bieging et al. 1989, Scuderi et al. 1998, Benaglia et al. 2001b, etc.) have been observed. Roughly half of them were detected at least at one frequency. When spectral indices could be measured, $``$40% showed non-thermal or composite indices.

In recent years the hypothesis was proposed (Chapman et al. 1999, Dougherty & Williams 2000) that non-thermal emission could be indicative of binarity. In that case, part of the emission comes from the colliding wind region (CWR) of two early massive stars. Chapman et al. (1999) observed 36 WR stars with the ATCA at 3, 6, 13 and 20 cm, obtaining an angular resolution of $`1\mathrm{"}`$ at 3 cm. Dougherty & Williams (2000) discussed the spectral index behaviour of 25 WR stars, by taking into account observations carried out mainly with the VLA, ATCA, MERLIN and the WSRT. Benaglia & Romero (2003) presented a table with the characteristics of the non-thermal or composite WR emitters.

Unfortunately, in some particular cases where the non-thermal emission is strong, not enough spectroscopic data are currently available to investigate the binary status. For example, CD–47 4551 (CPD–47 2963), an O4 III(f) star (Walborn 1982), was observed with the ATCA (Benaglia et al. 2001b) and detected with flux densities of $`1.77\pm 0.05`$ mJy at 3 cm and $`2.98\pm 0.05`$ mJy at 6 cm ($`\alpha =0.89\pm 0.06`$). If the star is a binary with another WR or OB star, the non-thermal emission could be produced at the CWR of both stars.

Non-thermal radio emission can also be produced by the stellar wind of a single star, as described by White (1985). In his model, electrons are accelerated to relativistic energies by shocks in the wind near the star, and emit radio radiation. For strong shocks, the expected spectral index is +0.5. Though this is similar to the thermal index value, the non-thermal flux density coming from a single stellar wind will suffer a break at a frequency of a few GHz, and behave as $`S(\nu )\mathrm{log}(\nu )`$ for larger frequencies, while the thermal flux will obey a $`\alpha =0.60.7`$ law up to the infrared range.

Radio observations are very useful to derive the mass loss rates of early-type stars (Leitherer et al. 1995). The source HD 93129A represents one of the most interesting candidates because it is the only known example of the hottest supergiant class. With the help of high resolution radio data, properties of the wind can be unveiled. In detecting the stellar wind, we intend to identify its radio spectrum, look for a possible colliding wind region, and compare the mass loss rate with those derived from other methods. The observations will also allow us the study of other early-type stars in the region.

The paper is organized as follows: Section 2 summarizes the characteristics of the target sources relevant to this study. Section 3 describes the observations and data reduction. In Section 4 we present the results, followed by the discussion in Section 5. Section 6 contains our conclusions.

## 2 The target stars in Trumpler 14

The field of view of the present data (see Section 3) corresponds to the cluster Trumpler 14. Lists of the catalogued members can be found in Massey & Johnson (1993), and recently DeGioia-Eastwood et al. (2001). A wide range of spectroscopic or photometric distances can be found for Tr 14 in the literature. Humphreys & McElroy (1984) derived a distance of 3.47 kpc. Morrell et al. (1988) discussed two possible distances, 2.8 (variable $`R`$, from Tapia et al. 1988) and 3.45 kpc ($`R`$ = 3.2), depending on the value taken for the ratio of total to selective absorption $`R`$. Massey & Johnson (1993) and Vázquez et al. (1996) found similar values of 3.2 and 3.1 kpc, respectively. Most recently, Tapia et al. (2003) computed spectroscopic parallaxes to derive a distance of 2.8 kpc. We adopt the latter value throughout this paper. At this distance, the projected cluster size of $`140\mathrm{"}`$ (Tapia et al. 2003) corresponds to 1.9 pc. Chandra observations towards the $`\eta `$ Carina region (Evans et al. 2003) revealed strong X-ray emission from the stars HD 93129A/B. Their combined X-ray luminosity is $`L_\mathrm{X}(0.52.4\mathrm{keV})=2.9\times 10^{33}`$ erg s-1, adjusted to a distance of 2.8 kpc.

Penny et al. (1993) presented radial velocity measurements for seven stars in Tr 14, with mean values of $`4.4\pm 7.4`$ km s-1 for HD 93129A, $`7.2\pm 2.4`$ km s-1 for HD 93129B and $`7.3\pm 3.3`$ km s-1 for HD 93128. They concluded that none of them is a strong binary candidate (see also Conti et al. 1977, 1979). Based on the standard deviation in the mean stellar radial velocities, they estimated that the typical random motions in this cluster are 8.5 km s-1. For HD 93129A Penny et al. placed an upper limit of $``$10 km s-1 on the semiamplitude of any short-period binary. To search for possible long-period binaries they obtained speckle observations, with negative results. Surprisingly, recent HST observations indicate that HD 93129A is a binary (Walborn 2002) with a separation of 55 mas (Walborn 2003, priv. comm.; Nelan et al., in prep.), which could have been detected by Penny et al. (1993).

The optical positions of the three earliest stars in the observed field were taken from the GOS catalog, and come from three different sources, varying in precision: the Tycho-2 catalog (Hog et al. 2000) lists a very accurate position ($`\sigma 5`$ mas) for HD 93129A, whereas the position of HD 93129B originates from the Washington Double Star catalog (Worley & Douglass 1997) and that of HD 93128 was measured from NTT data, according to Maíz-Apellániz & Walborn (2002). For a summary of the stellar parameters see Table 1.

#### HD 93129A.

The spectrum of this early O2 If\* star is very similar to Wolf-Rayet stars of spectral type WN6 or WN7 (Walborn et al. 2002). Taresch et al. (1997) studied the spectrum from optical to FUV range, concluding that its strong wind is radiatively driven. By fitting an H$`\alpha `$ profile they obtained a mass loss rate of $`1.8\times 10^5`$ M yr-1 compared to $`2.2\times 10^5`$ M yr-1 derived from H$`\alpha `$ measurements by Puls et al. (1996; assuming $`T_{\mathrm{eff}}`$ = 50,500 K).

#### HD 93129B and HD 93128.

The effective temperature of HD 93128 is $`52000\pm 1500`$ K (Vacca et al. 1996). For HD 93129B, which has the same spectral type as HD 93128, we adopt the same value. The stellar luminosities were derived taking the absolute visual magnitudes from Walborn et al. (2002).

#### HD 93250.

This star is a member of the Tr 16 cluster, located $`10^{}`$ East of HD 93129A/B (well outside our primary beam). Leitherer et al. (1995) reported a 3 cm flux density of $`1.36\pm 0.17`$ mJy and a 3$`\sigma `$ upper limit of 3.57 mJy at 6 cm. Assuming a purely thermal spectrum, a distance of 2.2 kpc, a wind velocity of $`v_{\mathrm{}}`$ = 3200 km s-1 and an effective temperature of $`T_{\mathrm{eff}}`$ = 50,000 K ($`T_\mathrm{e}`$ = 0.4 $`T_{\mathrm{eff}}`$ = 20,000 K according to Drew 1989) Leitherer et al. derived a mass loss rate of $`4.1\times 10^5`$ M yr-1. Considering a re-classified spectral type of O3.5 V((f+)) (Walborn et al. 2002) and a distance of 2.5 kpc to Tr 16 (Tapia et al. 2003) we revise the mass loss rate to $`5\times 10^5`$ M yr-1 (see Section 4.1).

## 3 Observations and Data Reduction

The observations presented here were obtained with the Australia Telescope Compact Array (ATCA) in January 2003 using the 6B configuration. The target field, which contains three massive Of stars, was observed simultaneously at two frequencies: 8.64 GHz (3 cm) and 4.80 GHz (6 cm), with a total bandwidth each of 128 MHz over 32 channels. The maximum baseline was 6 km, resulting in angular resolutions around $`1\mathrm{"}2\mathrm{"}`$. The size of the primary beam is $`10^{}`$ at 6 cm and $`5^{}`$ at 3 cm. The pointing position was $`\alpha ,\delta `$(J2000) = $`10^\mathrm{h}\mathrm{\hspace{0.17em}43}^\mathrm{m}\mathrm{\hspace{0.17em}56}^\mathrm{s}`$, $`59^{}\mathrm{\hspace{0.17em}32}^{}\mathrm{\hspace{0.17em}51}\mathrm{"}`$, slightly offset from the three stars HD 93129A/B and HD 93128. The field was observed during intervals of 8 to 15 min, depending on the atmospheric phase stability, interleaved with short observations of the nearby phase calibrator 1045–62, and spanning an LST range of 12 h. The total time on the target field was nearly 300 minutes.

The data were reduced in miriad using standard procedures. The flux density scale was calibrated with observations of the primary calibrator PKS B1934–638, assuming flux densities of 2.84 Jy at 3 cm and 5.83 Jy at 6 cm (Reynolds 1994). After data calibration the visibilities were Fourier-transformed using ‘natural’- and ‘uniform’-weighting. To overcome confusion due to diffuse emission from extended sources in or towards the observed fields, the shortest baselines were removed. For the same reason ‘uniform’-weighting resulted in the lowest noise levels: 0.11 and 0.21 mJy beam-1 at 3 and 6 cm, respectively. The synthesized beam is $`1.1\mathrm{"}\times 1.0\mathrm{"}`$ at 3 cm and $`1.9\mathrm{"}\times 1.7\mathrm{"}`$ at 6 cm.

## 4 Results

Figs. 1 and 2 show the central part of the resulting 3 and 6 cm radio continuum images. The optical positions of the three target stars (see Table 1) are marked.

Only one radio source, HD 93129A, is detected. Using the miriad task imfit we fit the position and flux density (of a point source) at both frequencies resulting in $`\alpha ,\delta `$(J2000) = $`10^\mathrm{h}\mathrm{\hspace{0.17em}43}^\mathrm{m}\mathrm{\hspace{0.17em}57}^\mathrm{s}.45`$, $`59^{}\mathrm{\hspace{0.17em}32}^{}\mathrm{\hspace{0.17em}51.36}\mathrm{"}`$, $`2.0\pm 0.2`$ mJy at 3 cm, and $`\alpha ,\delta `$(J2000) = $`10^\mathrm{h}\mathrm{\hspace{0.17em}43}^\mathrm{m}\mathrm{\hspace{0.17em}57}^\mathrm{s}.47`$, $`59^{}\mathrm{\hspace{0.17em}32}^{}\mathrm{\hspace{0.17em}51.23}\mathrm{"}`$, $`4.1\pm 0.4`$ mJy at 6 cm. From the 3 and 6 cm flux densities we derive a spectral index of $`\alpha =1.2\pm 0.3`$ ($`S_\nu \nu ^\alpha `$) indicating predominantly non-thermal emission. The measured flux densities for HD 93129A and the 3$`\sigma `$ upper limits for HD 93129B and HD 93128 are listed in Table 2.

The separation between the fitted 3 and 6 cm radio continuum positions for HD 93129A is $`0.2\mathrm{"}`$, slightly larger than the positional uncertainty (which is approximately the size of the synthesized beam divided by signal-to-noise ratio of the source). The radio and optical positions agree within the given uncertainties.

The Simbad database lists 250 Galactic objects within a radius of $`2.5^{}`$ around our pointing position (see Section 3), at least 11 of which are O-type stars, including the three target stars. Table 3 lists the eight additional stars, none of which was detected. The same $`3\sigma `$ upper limits apply.

### 4.1 Mass loss rates

According to Wright & Barlow (1975), the stellar mass loss rate, $`\dot{M}`$, of an optically thick ionized wind can be expressed in terms of the measured thermal flux density, $`S_\nu `$, at a frequency $`\nu `$, from radio to IR ranges, as:

$$\dot{M}=5.32\times 10^4\frac{S_\nu ^{3/4}d^{3/2}v_{\mathrm{}}\mu }{Z\sqrt{\gamma g_\nu \nu }}\mathrm{M}_{}\mathrm{yr}^1.$$
(1)
Here $`d`$ is the stellar distance in kpc, $`v_{\mathrm{}}`$ the wind terminal velocity in km s-1, $`S_\nu `$ the thermal flux density in mJy, and $`\nu `$ the frequency in Hz; $`\mu `$ stands for the mean molecular weight of the ions, $`\gamma `$ is the mean number of electrons per ion and $`Z`$ the r.m.s. ionic charge. For a fully ionized H-He gas, $`\mu =1.38`$. In a first approximation, $`Z\gamma 1`$. We will take those values for all calculations throughout the text. The Gaunt factor $`g_\nu `$ can be approximated by:

$`g_\nu =9.77\left(1+0.13\mathrm{log}{\displaystyle \frac{(0.4\times T_{\mathrm{eff}})^{3/2}}{Z\nu }}\right).`$ (2)
From the derived spectral index of $`\alpha =1.2\pm 0.3`$ we know that the observed radio continuum emission consists of thermal and non-thermal components. Their relative contributions (and locations) are unknown since the detected source is essentially unresolved. If $`f_\mathrm{T}`$ is the fraction of the thermal emission at 3 cm with $`S_{3\mathrm{c}\mathrm{m}}^\mathrm{T}=f_\mathrm{T}\times S_{3\mathrm{c}\mathrm{m}}`$, $`f_\mathrm{T}=01`$, then the mass loss rate of HD 93129A is $`\dot{M}=(f_\mathrm{T})^{3/4}\times 8.6\times 10^5`$ M yr-1. For $`f_\mathrm{T}1`$ this would be one of the highest mass loss rates ever derived for a massive, early-type star. For HD 93129B and HD 93128 we derive 3$`\sigma `$ upper limits of $`\dot{M}<(f_\mathrm{T})^{3/4}\times 2.13\times 10^5`$ M yr-1.

Given a stellar mass, luminosity, effective temperature, and wind terminal velocity, the models of Vink et al. (2000) can be used to compute a theoretical — or expected — mass loss rate, $`\dot{M}_{\mathrm{exp}}`$. Using Vink’s program (see astro.imperial.ac.uk/$``$jvink/) and the values given in Table 1, we calculate $`\dot{M}_{\mathrm{exp}}=2.12\times 10^5`$ M yr-1 for HD 93129A. While this value is in good agreement with the H$`\alpha `$-derived mass loss rates (Puls et al. 1996, Taresch et al. 1997), it is a factor four smaller than that derived from the measured 3 cm radio continuum flux density, suggesting $`f_\mathrm{T}=0.15`$. For HD 93129B and HD 93128 we calculate $`\dot{M}_{\mathrm{exp}}`$ of $`0.3\times 10^5`$ M yr-1.

Table 3 lists the expected stellar mass loss rates for the other eight O stars in the field, estimated as for the earliest stars and assuming a distance of 2.8 kpc. The luminosities, masses, and effective temperatures were taken from Vacca et al. (1996). The terminal velocities were interpolated from Table 3 of Prinja et al. (1990) according to spectral type. The expected thermal flux density at 3 cm was computed inverting Eq. (1), and is listed in the last column.

## 5 Discussion

### 5.1 A CWR between HD 93129A and B ?

The distance of $`D=2.6\mathrm{"}`$ = 7300 AU between coordinates of HD 93129A and B is a lower limit to their absolute separation. If a colliding wind region exists between the two stars, it would be located at projected distances $`r_\mathrm{A}`$, $`r_\mathrm{B}`$ (see Eichler & Usov, 1993):

$$r_\mathrm{A}=\frac{1}{1+\eta ^{1/2}}D,r_\mathrm{B}=\frac{\eta ^{1/2}}{1+\eta ^{1/2}}D.$$
(3)
The parameter $`\eta `$ is defined in terms of the wind terminal velocities and the stellar mass loss rates $`\dot{M}`$:

$$\eta =\frac{\dot{M_\mathrm{B}}\times v_{\mathrm{},\mathrm{B}}}{\dot{M_\mathrm{A}}\times v_{\mathrm{},\mathrm{A}}}.$$
(4)
Using the derived theoretical values for $`\dot{M_\mathrm{A}}`$ and $`\dot{M_\mathrm{B}}`$, we find $`\eta `$ = 0.14. Thus $`r_\mathrm{A}=1.9\mathrm{"}5350`$ AU, and $`r_\mathrm{B}=0.7\mathrm{"}1950`$ AU. The putative CWR would be centered around $`\alpha ,\delta `$(J2000)CWR = $`10^\mathrm{h}\mathrm{\hspace{0.17em}43}^\mathrm{m}\mathrm{\hspace{0.17em}57.59}^\mathrm{s}`$, $`59^{}\mathrm{\hspace{0.17em}32}^{}\mathrm{\hspace{0.17em}52.9}\mathrm{"}`$, but no emission was detected with the ATCA. The same 3$`\sigma `$ upper limits apply (see Table 2).

Assuming a spectral index of $`\alpha =0.5`$ (Eichler & Usov 1993), we derive an upper limit to the non-thermal radio luminosity from a CWR between stars A and B of $`L_{\mathrm{obs}}<1.4\times 10^{28}\mathrm{erg}\mathrm{s}^1`$. The expected non-thermal luminosity from a possible A/B binary is difficult to estimate due to the poor knowledge of system parameters such as size of the CWR and the magnetic field. Moreover, according to Wright & Barlow (1975), the radii of the radio photospheres at 3 and 6 cm would be $``$30 and 60 AU, ($`<<r_\mathrm{A}`$, $`r_\mathrm{B}`$), respectively, if $`\dot{M}=5\times 10^5`$ M yr-1. The existence of a CWR between HD 93129A and B is very unlikely (see Section 5.4).

### 5.2 A CWR between the binary components of HD 93129A ?

Walborn (2002) recently reported that HD 93129A is a binary system. Using the HST Fine Guidance Sensor Interferometer, Nelan et al. (2003, in prep.) resolved it into a 55 mas binary with a magnitude difference of 0.9 mag (Walborn 2003, priv. comm.). In the binary system, we will call “a” the O2 primary, and “b”, the secondary. The total $`V`$ magnitude measured for HD 93129A is 7.2 mag (Hog et al. 1998). Given that $`V_\mathrm{a}=7.6`$ mag, and $`V_\mathrm{b}=8.5`$ mag (Walborn 2003, priv. comm.), the secondary has a magnitude similar to that of HD 93129B ($`V=8.9`$ mag) and HD 93128 ($`V=8.8`$ mag) (Walborn et al. 2002). The four stars are probable coeval, and their spectral types should be similar (E. Nelan et al., in preparation).

The binary separation (55 mas = 154 AU) implies that the stellar winds are likely to interact. In fact, there is a high probability that the system is an O+O colliding wind binary. Additional support to this idea can be found by studying the radio spectral index (see Section 5.3). Assuming that the companion is an O3.5 V star with the same terminal velocity, effective temperature, luminosity, and expected mass loss rate than HD 93129B (see Table 1), we find $`\eta =0.14`$, $`r_{\mathrm{Aa}}=112`$ AU, and $`r_{\mathrm{Ab}}=42`$ AU. The binary separation of HD 93129A is similar to that of WR 146 (210 AU, Setia Gunawan et al. 2000), thus it is reasonable to assume a similar extent for the CWR, i.e. $``$50 AU. We estimate the equipartition magnetic field for the CWR according to Miley (1980) and find $`B_{\mathrm{CWR}}4`$ mG. Using the magnetic field structure given by Eichler & Usov (1993), the surface magnetic field of HD 93129Aa can be extrapolated to $`B_{}100`$ G, if the stellar rotational velocity is $`0.1v_{\mathrm{}}`$ (Conti & Ebbets 1997). This estimate is consistent with the few observational values available (Donati et al. 2001, 2002).

Substituting the parameters of HD 93129Aa and Ab into Eqs. (15) and (16) from Eichler & Usov (1993) we obtain:

$`L_{\mathrm{rad}}`$ $`7.4\times 10^{31}\beta \left({\displaystyle \frac{ϵ}{10^2}}\right)\left({\displaystyle \frac{\xi }{10^3}}\right)\mathrm{erg}\mathrm{s}^1`$ (5)
for the expected non-thermal radio luminosity of the CWR, and $`\xi 5\times 10^5`$. Taking typical values of $`\beta 0.1`$, $`ϵ10^2`$, we calculate $`L_{\mathrm{rad}}4\times 10^{29}\mathrm{erg}\mathrm{s}^1`$.

From the observed flux densities we obtain a radio luminosity of $`L_{\mathrm{obs}}=1.12\times 10^{29}`$ erg s-1 for HD 93129A between 4.8 and 8.64 GHz, which is in good agreement with the total expected one.

In a colliding wind binary system like the one described here, a total X-ray luminosity due to Bremsstrahlung can be computed. Following Usov (1992) we obtain $`2.2\times 10^{33}`$ erg s-1, very close to the observed X-ray luminosity by Evans et al. (2003).

### 5.3 Spectral index considerations

The spectral index of $`1.2\pm 0.3`$ derived from the observations presented here is much steeper than that estimated by White (1985), $`\alpha =+0.5`$, for synchrotron emission from a single star. Steeper indices can result in the CW scenario as the effect of the softening of the particle spectrum due to inverse Compton (IC) losses (Benaglia & Romero 2003). At a CWR the accelerated electrons follow a power-law spectrum $`N(E)E^p`$, with $`p=2`$. Both synchrotron and IC processes help electrons to lose energy. At high energies, the radiative cooling would be dominated by IC losses, due to the strong photon fields at the CWR. At the energy $`E_\mathrm{b}`$ at which the cooling and escape times are equal (e.g. Longair 1997, p. 281), IC losses produce a break in the spectrum, from an index $`p`$ to $`p^{}=p1`$. The same population of relativistic electrons is capable of losing energy via synchrotron emission, producing a flux $`S(\nu )\nu ^\alpha `$, where $`\alpha =(p+1)/2`$. Thus, for $`E<E_\mathrm{b}`$, $`\alpha 0.5`$, and for $`E>E_\mathrm{b}`$, $`\alpha 1`$. This last value seems to be consistent with the spectral index derived from the present observations. The break energy can be computed if the following set of parameters is known (see Benaglia & Romero 2003): the mass loss rates, wind terminal velocities, and luminosities of both components, and the magnetic field at the CWR, together with the size $`s`$ of the non-thermal source. Replacing the values adopted or obtained above, the break frequency in the radio range will be a few GHz. ATCA observations at 1.4 and 2.4 GHz are necessary to test this scenario.

The relative importance of the processes already mentioned can be measured by estimating the time needed for an electron to lose all its energy, either in the form of synchrotron radiation ($`t_{\mathrm{syn}}`$) or IC scattering ($`t_{\mathrm{IC}}`$). The importance of synchrotron and inverse Compton radiation for electrons of the same energy can be measured computing the ratio between the corresponding lifetimes (Longair 1997, p.274):

$$\frac{t_{\mathrm{syn}}}{t_{\mathrm{IC}}}=\frac{U_{\mathrm{rad}}}{U_{\mathrm{mag}}},$$
(6)
where $`U_{\mathrm{rad}}`$ and $`U_{\mathrm{mag}}`$ are the energy densities of the radiation field and of the magnetic field respectively. Using the observed radio luminosity $`L_{\mathrm{obs}}`$, we get $`U_{\mathrm{rad}}=2.5\times 10^9`$ eV m-3. The variable $`U_{\mathrm{mag}}`$ is $`6.4\times 10^{15}`$ eV m-3. The corresponding time ratio is about $`4\times 10^5`$, thus $`t_{\mathrm{syn}}>>t_{\mathrm{IC}}`$. Following this last estimate, and under the conditions mentioned above, IC losses should dominate over synchrotron ones.

### 5.4 Colliding wind emission from other binaries

There are only a few examples (WR 146, WR 147, Cyg OB2 No. 5) of early-type binary systems in the literature where radio emission from the CWR has been observed. WR 146 is a WC6+O8 binary (van der Hucht 2001). Dougherty et al. (2000) observed the stars with MERLIN at 5 GHz twice, and with the VLA at 1.4, 5, 8.5 and 22 GHz. The combination of their MERLIN data with the 22 GHz VLA data allowed them to detect two thermal sources, coincident with the WR and OB stars, plus a non-thermal one in-between, nearer to the OB star, and possibly coming from the CWR. At a distance of 1.25 kpc (van der Hucht et al. 2001), the separation of the stellar components is $``$210 AU, and the size of the non-thermal source is $``$50 AU (Setia Gunawan et al. 2000). The flux densities quoted in the literature are: $`31.4\pm 0.4`$ mJy and $`28.5\pm 0.3`$ mJy at 5 GHz for the source identified as the CWR, $`7.0\pm 1.3`$ mJy at 22 GHz for the emission possibly related to the WR star, and $`10.4\pm 1.0`$ mJy from the OB star. The ratios of the non-thermal to thermal flux densities are $``$3 for the O star, and $``$4 for the WR star.

Contreras & Rodríguez (1999) presented 8.5 GHz-VLA observations of the star WR 147, a binary system formed by WN8 (h)+B0.5 V components (van der Hucht et al. 2001), at a distance of 650 pc (van der Hucht 2001). They detected a southern, thermal radio source coincident with the WR component, and a northern, non-thermal source. The latter is interpreted as originating at the interaction zone between the winds of the WR and OB stars. Their averaged measured flux densities are $`28.4\pm 0.5`$ mJy from the thermal source, and $`10.4\pm 0.5`$ mJy from the non-thermal source. The stellar separation is $``$420 AU (Setia Gunawan et al. 2001). The non-thermal to thermal flux density ratio is $``$0.4.

Cyg OB2 No. 5 seems to be formed by three stars: an O7 Ia+Ofpe/WN9 contact binary, and a B0 V star in a larger orbit (Contreras et al. 1997, Rauw et al. 1999), located at 1.8 kpc (e.g. Waldron et al. 1998). The separation between the contact pair and the third component is about 1700 AU. Radio observations by Persi et al. (1990) and Miralles et al. (1994) at 1.4, 5, 8.5 and 15 GHz, showed that the source coincident with the contact pair presented flux variations, defining a “high” emission state ($`S_{8.5\mathrm{GHz}}`$ = 5 – 7 mJy, and $`\alpha 0`$), and a “low” emission state ($`S_{8.5\mathrm{GHz}}`$ = 1 –2 mJy, and $`\alpha 0.6`$). Contreras et al. (1996, 1997) performed VLA observations at 5 and 8.5 GHz, detecting a non-thermal source ($`\alpha =2.4\pm 0.6`$) near the B0 V star, and proposed it was a CWR in the system. They measured a flux density at 8.5 GHz of $`0.31\pm 0.02`$ mJy. The ratio of the emission from the weaker non-thermal source to the emission of the contact pair at the lower state, is about 0.2.

The system HD 93129A presents a binary separation similar to that of WR 146. The ratio of non-thermal to thermal emission is unknown. If the magnetic field at the CWR is equal to the equipartition field derived above, $``$4 mG, then it is similar to those derived for WR 147 (5 mG) and WR 146 (25 mG); see Benaglia & Romero (2003). The scarcity of the examples and the huge differences in the ratios of CWR-non-thermal to stellar-thermal flux densities measured for the three cases from the literature, make the results of the comparison inconclusive.

### 5.5 The radio spectrum of HD 93129A

To separate the thermal and non-thermal contributions to the spectrum of HD 93129A, further radio observations are needed at lower and, particularly, at higher radio frequencies. These can be obtained with the ATCA at 1.4, 2.3, 16–25, and 85–91 GHz. The system can be resolved with the Long Baseline Array, but the flux densities of the thermal winds are likely too low for a detection, and an absolute position for the CWR would be difficult to obtain.

If the actual mass loss rate of HD 93129A is $`\dot{M}=2.1\times 10^5`$ M yr-1, i.e. $`f_\mathrm{T}`$ = 0.15 (see Section 4.1), the thermal and non-thermal contributions to the 3 cm radio flux density are $``$0.3 and $``$1.7 mJy, respectively. Assuming a spectral index of $`\alpha _\mathrm{T}`$ = 0.6 for the thermal component, we find the spectral index for the non-thermal component is $`\alpha _{\mathrm{NT}}1.4\pm 0.3`$. The expected total flux densities (see Fig. 3a) are $``$1.2, $`1`$, and $``$1.3 mJy at 16, 25 and 90 GHz, respectively.

In contrast, for $`f_\mathrm{T}`$ = 0.5 we find $`\alpha _{\mathrm{NT}}2.1\pm 0.3`$, similar to that derived for other Of systems (e.g., Cyg OB2 No. 5, Contreras et al. 1996, 1997). At 16, 25 and 90 GHz we expect $``$1.7, $``$2, and $``$4.1 mJy (see Fig. 3b). In this case, the radio derived mass loss rate of the system is $`5.1\times 10^5`$ M yr-1. The two scenarios can be distinguished with high-frequency ATCA synthesis observations.

## 6 Conclusions

The field of view towards the Galactic cluster Tr 14 has been observed with the ATCA in the radio continuum at 4.80 GHz (3 cm) and 8.64 GHz (6 cm). Only one source, the O2 If\* star HD 93129A, was detected with flux densities of $`2.0\pm 0.2`$ mJy at 3 cm, and $`4.1\pm 0.4`$ mJy at 6 cm, resulting in a very steep spectral index of $`\alpha =1.2\pm 0.3`$. This star, recently reported as binary, is the earliest known O-star and potentially the most massive system in the Galaxy, followed by $`\eta `$ Carina. The observed radio emission is most likely coming from (a) the thermal winds of both binary components and (b) the non-thermal emission from the colliding wind region. The negative spectral index of HD 93129A suggests that the latter is the dominant emission mechanism.

We derive a mass loss rate of $`(f_\mathrm{T})^{3/4}\times 8.6\times 10^5`$ M yr-1 for HD 93129A and 3$`\sigma `$ upper limits of $`(f_\mathrm{T})^{3/4}\times 2.13\times 10^5`$ M yr-1 for HD 93129B and HD 93128, where $`f_\mathrm{T}`$ is the thermal fraction of the 3 cm radio continuum emission.

Although the secondary star in the binary HD 93129A remains yet unclassified, the presence of intense non-thermal emission and high radio luminosity together suggest it is another massive, early-type star.

The lack of radio emission near HD 93129B, and the extremely large separation of HD 93129A/B (7300 AU at 2.8 kpc) conspires against the existence of a colliding wind region in-between the stars.

It is crucial to perform observations (i) at lower frequencies, where the thermal contributions should be negligible, and an accurate non-thermal index can be derived; this would help to investigate the processes responsible for the radiation detected; and (ii) at high frequencies, at which non-thermal emission decreases, and the estimation of the thermal flux will enable us to derive the value of the true mass loss rate of HD 93129A.

## Acknowledgments

We would like to thank the anonymous referee for comments and suggestions that help to improve the presentation of our investigation. We are indebted to E. Nelan and N. A. Walborn for allowing us to share their HST results before publishing, and to N. A. Walborn and G. E. Romero for a critical reading of the manuscript. P.B. is grateful to the ATNF staff at Sydney and Narrabri. This research has been supported mainly by Fundación Antorchas (P.B.), and has made use of the SIMBAD database operated at CDS, Strasbourg, France.

## References

Abbott, D.C., Bieging, J.H., & Churchwell, E. 1986, ApJ 303, 239
Benaglia, P., & Cappa, C.E. 1999, A&A 346, 979
Benaglia, P., Romero, G.E., Stevens, I.R., & Torres, D.F. 2001a, A&A 366, 605
Benaglia, P., Cappa, C.E., & Koribalski, B.S. 2001b, A&A 372, 952
Benaglia, P., & Romero, G.E. 2003, A&A 399, 1121
Bohannan, B. 1990, in “Properties of Hot Luminous Stars”, ed. C.D. Garmany, Brigham Young Press, 39
Chapman, J.M., Leitherer, C., Koribalski, B.S., et al. 1999, ApJ 518, 890
Chlebowski, T., & Garmany, C.D. 1991, ApJ 368, 241
Conti, P.S., Leep, E.M., & Lorre, J.J. 1977, ApJ 214, 759
Conti, P.S., Niemela, V.S., & Walborn, N.R. 1979, ApJ 228, 206
Conti, P.S., & Garmany, C.D. 1983, PASP 95, 411
Conti, P.S., & Ebbets, D. 1997, ApJ 213, 438
Contreras, M.E., Rodríguez, L.F., Gómez, Y., et al. 1996, ApJ 469, 329
Contreras, M.E., Rodríguez, L.F., Tapia, M., et al. 1997, ApJ 488, L153
Contreras, M.E., & Rodríguez, L.F. 1999, ApJ 515, 762
Crowther, P.A., Smith, L.J., Hillier, D.J., & Schmutz, W. 1995, A&A 293, 427
DeGioia-Eastwood, K., Throop, H., Walker, G., et al. 2001, ApJ 549, 578
Donati, J.-F., Wade, G.A., Babel, J., et al. 2001, MNRAS 326, 1265
Donati, J.-F., Babel, J., Harries, T.J., et al. 2002, MNRAS 333, 55
Dougherty, S.M., Williams, P.M., & Pollaco, D.L. 2000, MNRAS 316, 143
Dougherty, S.M., & Williams, P.M. 2000, MNRAS 319, 1005
Drew, J.E. 1989, ApJS 71, 267
Eichler, D., & Usov, V. 1993, ApJ 402, 271
Evans, N.R., Seward, F.D., Kraus, M.I., et al. 2003, ApJ 589, 509
García, B. Malaroda, S., Levato, H., et al. 1998, PASP, 110, 53
Hog, E., Kuzmin, A., Bastian, U., et al. 1998, A&A 335, L65 (Tycho-1)
Hog, E., Fabricius, C., Makarov, V.V., et al. 2000, A&A 355, L27 (Tycho-2)
Humphreys, R., & McElroy, D.B. 1984, ApJ 284, 565
van der Hucht, K.A. 2001, New Astronomy 45, 135
van der Hucht, K.A., Setia Gunawan, D.Y.A., Williams, P.M., et al. 2001, in “Interacting Winds from Massive Stars”, eds. A.F.J. Moffat & N. St-Louis, ASP Conf. Ser. 165, p. 267
van der Hucht, K.A. 2002, Astrophysics & Space Science 281, 199
Lamers, H.J.G.L.M., & Leitherer, C. 1993, ApJ 412, 771
Langer, N., Hamann, W.-R., Lennon, M. et al. 1994, A&A 290, 819
Leitherer, C., Chapman, J.M., & Koribalski, B.S. 1995, ApJ 450, 289
Levato, H., Malaroda, S., Morrell, et al. 2000, PASP 112, 359
Longair, M.S. 1997, High Energy Astrophysics, Vol. 2, Cambridge University Press, Cambridge
Maíz-Apellániz, J., & Walborn, N. R. 2002, in “A Massive Star Odyssey: from Main Sequence to Supernova”, eds. K. A. van der Hucht, A. Herrero, & C. Esteban, ASP Conf. Ser. 212, p. 560
Mason, B.D., Gies, D.R., Hartkopf, W.I., et al. 1998, ApJ 115, 821
Mason, B.D., Wycoff, G.L., Hartkopf, W.I., et al. 2001, AJ 122, 3466
Massey, P., & Johnson, J. 1993, AJ 105, 980
Massey, P., & Hunter, D.A. 1998, ApJ 493, 180
Miley, G.K. 1980, Ann. Rev. Astron. Astrophys. 18, 165
Miralles, M.P., Rodríguez, L.F., Tapia, M. et al. 1994, A&A 282, 547
Morrell, N., García, B., & Levato, H. 1988, PASP 100, 1431
Panagia, N., & Felli, M. 1975, A&A 39, 1
Penny, L.R., Gies, D.R., Hartkopf, W.I., et al. 1993, PASP 105, 588
Persi, P., Tapia, M., Rodríguez, L.F., et al. 1990, A&A 240, 93
Prinja, R.K., Barlow, M.J., & Howarth, I.D. 1990, ApJ 361, 607
Puls, J., Kudritzki, R.P., Herrero, A., et al. 1996, A&A 305, 171
Rauw, G., Vreux, J.-M., & Bohannan, B. 1999, ApJ 517, 416
Reed, B.C., & Beatty, A.E. 1995, ApJS 97, 189
Reynolds, J. 1994, “A Revised Flux Scale for the AT Compact Array”, ATNF Internal Report, AT/39.3/040
Scuderi, S., Panagia, N., Stanghellini, et al. 1998, A&A 332, 251
Setia Gunawan, D.Y.A., de Bruyn, A.G., van der Hucht, K.A., & Williams, P.M. 2000, A&A 356, 676
Setia Gunawan, D.Y.A., van der Hucht, K.A., Williams, P.M., et al. 2001, A&A 376, 460
Simon, K.P., Jonas, G., Kudritzki, R.P. et al. 1983, A&A 125, 34
Tapia, M., Roth, M., Vázquez, R.A. et al. 2003, MNRAS 339, 44
Taresch, G., Kudritzki, R.P., Hurwitz, M., et al. 1997, A&A 321, 531
Usov, V. 1992, ApJ 389, 635
Vacca, W.D., Garmany, C.D., & Schull, J.M. 1996, ApJ 460, 914
Vázquez, R.A., Baume, G., Feinstein, A., et al. 1996, A&AS 116, 75
Vink, J.S., de Koter, A., & Lamers, H.J.G.L.M. 2000, A&A 362, 295
Walborn, N.R. 1973, ApJ 179, 517
Walborn, N.R. 1982, AJ 87, 1300
Walborn, N.R. 1992, ApJ 393, L13
Walborn, N.R. 2002, in “A Massive Star Odyssey: from Main Sequence to Supernova”, eds. K. A. van der Hucht, A. Herrero, & C. Esteban, ASP Conf. Ser. 212, p. 13
Walborn, N.R., Howarth, I.D., Lennon, D.J., et al. 2002, AJ 123, 2754
Waldron, W.L., Corcoran, M.F., Drake, S.A. & Smale, A.P. 1998, ApJS 118, 217
Worley, C.E., & Douglass, G.G. 1997, A&AS 125, 523
White, R.L. 1985, ApJ 289, 698
Wright, A.E., & Barlow, M.J. 1975, MNRAS 170, 41


format: text

text:
# Probing the EOS of dense neutron-rich matter with high-energy radioactive beams

## 1 INTRODUCTION

With isospin symmetry conserving interactions, the $`\mathrm{EOS}`$ of neutron-rich matter of isospin asymmetry $`\delta (\rho _n\rho _p)/(\rho _n+\rho _p)`$ can be written as (see e.g., ), $`e(\rho ,\delta )=e(\rho ,0)+E_{sym}(\rho )\delta ^2+𝒪(\delta ^4)`$, where $`e(\rho ,\delta )`$ and $`e(\rho ,0)`$ is the energy per nucleon in isospin asymmetric and symmetric nuclear matter, respectively. The density dependence of the symmetry energy $`E_{sym}(\rho )`$ is the most important issue about the EOS of neutron-rich matter. In fact, to determine the $`E_{sym}(\rho )`$ has been a longstanding goal of extensive research with various microscopic and/or phenomenological models over the last few decades. However, the predicted $`E_{sym}(\rho )`$ are still rather divergent especially at high densities. Theoretical results can be roughly classified into two groups, i.e., a group where the $`E_{sym}(\rho )`$ rises monotonically with density and one in which it begins to fall above about $`2\rho _0`$. The high density behaviour of $`E_{sym}(\rho )`$ is actually considered by some astrophysicists as the most uncertain property of dense neutron-rich matter.

The importance of probing the EOS of neutron-rich matter lies in both nuclear physics and astrophysics. First of all, new structures of dripline nuclei, such as the appearance of new magic numbers, and the reaction mechanisms of radioactive nuclei are all strongly influenced by the isospin-dependence of the nuclear EOS. It also has profound consequences on many important issues in astrophysics. These include nucleosynthesis in pre-supernova evolution of massive stars, mechanisms of supernova explosions, neutrino flux, kaon condensations and the hadron-QGP phase transition in neutron stars, etc.

## 2 ISOSPIN ASYMMETRY IN NEUTRON STARS AND HEAVY-ION COLLISIONS

In a neutron star, several of its main properties are affected significantly by the correlation between the density and the proton fraction $`x_\beta `$ at $`\beta `$ equilibrium. The latter is determined by $`\mathrm{}c(3\pi ^2\rho x_\beta )^{1/3}=4E_{\mathrm{sym}}(\rho )(12x_\beta ).`$ The equilibrium proton fraction is therefore entirely determined by the $`E_{sym}(\rho )`$. We use the following two forms of the symmetry energy for an illustration $`E_{sym}^a(\rho )E_{sym}(\rho _0)u`$ and $`E_{sym}^b(\rho )0.5E_{sym}(\rho _0)u(3u),`$ where $`u\rho /\rho _0`$. The values of $`\delta _\beta =12x_\beta `$ corresponding to these two forms of symmetry energy are shown in the insert of Fig. 1. With the $`E_{sym}^b(\rho )`$, the $`\delta _\beta `$ is $`1`$ for $`\rho /\rho _03`$, indicating that the neutron star has become a pure neutron matter at these high densities. On the contrary, with the $`E_{sym}^a(\rho )`$, the neutron star becomes more proton-rich as the density increases. To see the connection between neutron stars and heavy-ion collisions, using a nuclear transport model we study in Fig. 1 the isospin asymmetry $`\delta _{like}`$ in central $`{}_{}{}^{132}Sn+^{124}Sn`$ reactions. Effects due to the different symmetry energies are clearly revealed especially at high densities. An astonishing similarity is seen in the resultant $`\delta \rho `$ correlations for neutron stars and heavy-ion collisions since the same underlying nuclear $`\mathrm{EOS}`$ is at work in both cases.

## 3 $`\pi ^{}/\pi ^+`$ PROBE OF THE EOS of DENSE NEUTRON-RICH MATTER

It is qualitatively easy to understand why the $`\pi ^{}/\pi ^+`$ ratio in heavy-ion collisions can be used to extract crucial information about the EOS of neutron-rich matter. On one hand, within the $`\mathrm{\Delta }`$ resonance model for pion production from first-chance independent nucleon-nucleon collisions, the primordial $`\pi ^{}/\pi ^+`$ ratio is $`(5N^2+NZ)/(5Z^2+NZ)(N/Z)^2`$. It is thus sensitive to the isospin asymmetry $`(N/Z)_{dense}`$ of dense matter in the participant region of heavy-ion collisions. On the other hand, within the statistical model for pion production, the $`\pi ^{}/\pi ^+`$ ratio is proportional to $`\mathrm{exp}\left[(\mu _n\mu _p)/T\right]`$, where T is the temperature. The difference in nucleon chemical potentials $`\mu _n\mu _p`$ is proportional to the strength of the symmetry potential. Therefore from both pictures the $`\pi ^{}/\pi ^+`$ ratio is useful for extracting the EOS of dense neutron-rich matter.

Shown in Fig. 2 are the isospin asymmetry of the high density region $`(n/p)_{\rho \rho _0}`$ (left) and the $`(\pi ^{}/\pi ^+)_{like}`$ ratio (right) $`(\pi ^{}/\pi ^+)_{like}(\pi ^{}+\mathrm{\Delta }^{}+\frac{1}{3}\mathrm{\Delta }^0)/(\pi ^++\mathrm{\Delta }^{++}+\frac{1}{3}\mathrm{\Delta }^+)`$ as a function of time for the central $`{}_{}{}^{132}Sn+^{124}Sn`$ reaction at beam energies from 200 to 1000 MeV/nucleon. This ratio naturally becomes the final $`\pi ^{}/\pi ^+`$ ratio at the freeze-out when the reaction time $`t`$ is much longer than the lifetime of the delta resonance $`\tau _\mathrm{\Delta }`$. The $`(\pi ^{}/\pi ^+)_{like}`$ ratio is rather high in the early stage of the reaction because of the large numbers of neutron-neutron scatterings near the surfaces where the neutron skins of the colliding nuclei overlap. It is seen that a variation of about 30% in the $`(n/p)_{\rho /\rho _0}`$ due to the different $`E_{sym}(\rho )`$ results in about 15% change in the final $`\pi ^{}/\pi ^+`$ ratio. It has thus an appreciable response factor of about 0.5 to the variation of n/p ratio and is approximately beam energy independent. Therefore, the $`(\pi ^{}/\pi ^+)_{like}`$ ratio is a sensitive probe of the high density behaviour of nuclear symmetry energy. Shown in Fig. 3 is the reaction system dependence of the $`\pi ^{}/\pi ^+`$ ratio. The symmetry energy effects increase with the isospin asymmetry energy as one expected.

In summary, within an isospin-dependent transport model we found that the $`\pi ^{}/\pi ^+`$ ratio is a promising probe of the EOS of dense neutron-rich matter. This work was supported by the National Science Foundation under grant No. PHY-0088934 and PHY-0243571.


text:
# Roughness of Sandpile Surfaces

## Abstract

We study the surface roughness of prototype models displaying self-organized criticality (SOC) and their noncritical variants in one dimension. For SOC systems, we find that two seemingly equivalent definitions of surface roughness yields different asymptotic scaling exponents. Using approximate analytical arguments and extensive numerical studies we conclude that this ambiguity is due to the special scaling properties of the nonlinear steady state surface. We also find that there is no such ambiguity for non-SOC models, although there may be intermediate crossovers to different roughness values. Such crossovers need to be distinguished from the true asymptotic behaviour, as in the case of a noncritical disordered sandpile model studied in .

Since the original proposal by Bak, Tang and Wisenfeld  there has been a large body of work directed towards understanding ubiquity of scale invariance in externally driven open dissipative systems using the concept of self-organized criticality (SOC). The sandpile is a prototype model system which has been extensively used as a paradigm of SOC . The principal aim has been to elucidate how slowly driven dissipative systems with fast relaxation mechanisms display long-tailed distributions of activity sizes. Recently, there have been resurge of interest in sandpile models in order to understand SOC in connection with better understood scale invariance in other nonequilibrium phenomena such as absorbing state phase transitions  and driven interfaces .

In , Krug, Socolar and Grinstein (KSG) studied the surface fluctuations in a prototype model of SOC, the limited local sandpile (LLS)  and its variations in one spatial dimension and concluded that the interfacial fluctuations, although nontrivial, are evidently unconnected to the criticality of the system. In fact, they argued that the evolution of height fluctuations $`\stackrel{~}{h}(x,t)`$ can be described by an extension of the Kardar-Parisi-Zhang (KPZ)  equation for an anchored interface

$$_t\stackrel{~}{h}=D_x^2\stackrel{~}{h}+c_x\stackrel{~}{h}+\lambda (_x\stackrel{~}{h})^2+\eta (x,t),$$
(1)
where $`\eta (x,t)`$ is a Gaussian white noise. The linear term $`c_x\stackrel{~}{h}`$ is more relevant than the diffusion term $`D_x^2\stackrel{~}{h}`$ and the KPZ nonlinearity $`\lambda (_x\stackrel{~}{h})^2`$ (in RG sense) and, since the interface is anchored, can not be eliminated by a Galilean shift. This term is responsible for transporting fluctuations up the pile thus relating the spatial roughness of the anchored interface to the dynamical roughening of a moving KPZ interface in $`d=1`$, $`\alpha _{LLS}=\beta _{KPZ(1d)}=1/3`$. This argument can be easily extended to higher dimensions.

In this paper we study interface roughness in LLS and related models and find that effect of the criticality of the system is somewhat subtle and is due to the nontrivial nature of the steady state surface. The value of surface roughness depends on whether it is measured about this nontrivial steady state profile or about the mean instantaneous surface, the latter being the correct choice for moving interfaces. The assertion made by KSG is valid only if $`\stackrel{~}{h}`$ in (1) is defined as the fluctuation around the nonlinear steady state surface. We also study modifications of the LLS which have nonlinear steady state profiles but which do not display SOC and find that there is no such ambiguity in the asymptotic roughness exponent. Thus, we conjecture that the roughness exponent is uniquely defined for systems not displaying SOC even though they may possess nontrivial steady state surfaces . For SOC systems the ambiguity exists and one needs to define the roughness appropriately. Our results have important bearings on studies such as the recent one by Barker and Mehta  who observed anomalously large roughness exponent in a sandpile model with structural disorder.

The limited local sandpile model in $`d=1`$ is defined as follows (see Fig. 1a): on a one dimensional lattice of sites $`i=0,1,\mathrm{},L`$ we define an integer height variable $`h_i`$. One grain of sand added at a randomly chosen site increases the height at that site by one: $`h_ih_i+1`$. The configuration is stable if the local slopes satisfy $`z_i=h_{i+1}h_iz_c`$ for all $`i`$, where we chose $`z_c=2`$ (our results are essentially unchanged for $`z_c2`$). An instability occurs when by addition of grains at site $`i+1`$, the local slope exceeds threshold $`z_i>z_c`$; in this case $`z_c`$ grains are transferred from column $`i+1`$ to column $`i`$. Subsequently, columns $`i`$ and $`i+2`$ may become unstable, setting off further topplings leading to an avalanche. The grains involved in an avalanche leave the system if they reach site $`i=0`$ ($`h_0(t)=0,t`$). A new grain is added after the system has attained a stable configuration. The final configuration reached is independent of the order in which the sites are updated in case more than one site become unstable. We count one unit of time (Monte Carlo step) for every $`L`$ grains added.

The standard way of quantifying avalanches in these systems is to count the number $`m`$ of grains that leave the pile after each deposition . The criticality of the model is reflected in a broad probability distribution of avalanche sizes, $`P(m)`$ (Fig. 1b), which has no other length scale except the system size $`L`$.

In , KSG also introduced and studied the inertial limited local sandpile (ILLS) which mimics the effect of inertia of the falling grains in a real sand pile. The instability condition setting off an avalanche is same as that in the LLS but the condition of stopping is changed. A cluster of grains when first destabilized is assigned a momentum $`p=0`$ and each time the front of the cluster reaches an unstable site ($`z_i>z_c`$) it gathers momentum, $`pp+1`$. If it comes across a stable site, $`p`$ decreases by an amount $`r`$: $`ppr`$. The cluster continues moving down as long as $`p>0`$ and leaves the pile if it reaches the site $`i=0`$. Clearly, $`r=\mathrm{}`$ corresponds to the LLS. It was noted by KSG that ILLS is not critical for any finite value of $`r`$ which is reflected in the corresponding avalanche size distribution (Fig. 1b).

To study the interfacial fluctuations in these models, we start with an initially flat pile ($`h_i(0)=0,i`$) and add grains till the pile reaches a steady state with the mean surface making a critical angle $`\psi `$ ($`\mathrm{tan}\psi _L\mathrm{}=3/2`$) with the horizontal. The closed boundary condition at $`i=L`$ ensures that the pile has only one surface with the critical slope . The width of the interface, in the steady state, can be measured in two ways

$`W_1^2(L)`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{i=1}{\overset{L}{}}}[h_i(t)h_i]^2\text{and,}`$ (2)
$`W_2^2(L)`$ $`=`$ $`{\displaystyle \frac{1}{L}}{\displaystyle \underset{i=1}{\overset{L}{}}}[h_i(t)s(t)i]^2,`$ (3)
where $`s(t)=2[_{i=1}^Lh_i(t)]/L(L+1)`$ is mean slope of the interface about which sum of instantaneous height fluctuations vanish $`_{i=1}^L[h_i(t)s(t)i]=0`$. In both definitions (3) above the ensemble average $`\mathrm{}`$ is identical to the time average in the steady state.

For large enough system sizes the scaling hypothesis is expected to be valid and the roughness exponent $`\alpha `$ is defined through

$$W_{1,2}(L)L^\alpha ,L\mathrm{}.$$
(4)
In Fig. 2, we plot the two widths $`W_1`$ and $`W_2`$ for a range of system sizes $`L=2^n`$ with $`6n14`$ (that is $`64L16384`$). It is seen that the two widths have different asymptotic behaviours and give two different values of $`\alpha `$. The roughness $`\alpha 0.33`$ computed from $`W_1`$ (and $`W_2`$ up to $`LL_c`$) is in accordance with the predictions of Eq.(1). The asymptotic roughness computed from $`W_2`$, beyond the crossover system size $`L_c1024`$, $`\alpha =0.65\pm 0.02`$ is larger than that for an interface generated by a simple random walk. This is surprising since, prima facie, the definitions in (3) are expected to be equivalent as far as asymptotic scaling is concerned (). In fact, for the ILLS, we have computed $`W_1`$ and $`W_2`$ for $`r=20`$ (Fig. 2, inset) and it is seen that they indeed have the same asymptotic behaviour and hence, a unique value of $`\alpha 0.33`$.

We now show that the crossover in $`W_2`$ in the LLS can be traced to the fact that the time averaged steady state profile is not linear, i.e., $`h_i(t)\overline{s}i`$ and this makes the two definitions nonequivalent. Although, this nonlinearity of $`h_i(t)`$ is a consequence of the boundary effects and is present for the ILLS as well, we demonstrate that the presence of SOC in LLS results in special scaling properties which is responsible for the observed differences between asymptotic scalings of $`W_1`$ and $`W_2`$. It can be shown that $`W_1`$ and $`W_2`$ are related by

$$W_2^2(L)+W_c^2(L)=W_1^2(L)+W_S^2(L),$$
(5)
where $`W_S(L)`$ is the root-mean-square (rms) wandering of the steady state interface profile $`h_i`$,

$$W_S^2(L)\frac{1}{L}\underset{i=0}{\overset{L}{}}\left[h_i\overline{s}i\right]^2,$$
(6)
with $`\overline{s}=s(t)`$ being the mean slope of the steady state profile satisfying $`_{i=0}^L[h_i\overline{s}i]=0`$. The second term on the left of (5), $`W_c^2(L)=[s(t)\overline{s}]^2(L+1)(2L+1)/6`$, is the contribution to interface width due to instantaneous slope fluctuations and is the analogue of the center of mass fluctuations in the case of moving interfaces . In fact, in the case of moving interfaces $`W_c`$ dominates for large length scales and hence $`W_2`$ is taken as the width. As it turns out, for the anchored interfaces we deal with here, $`W_c`$ is negligible compared to the other terms in (5) as $`L\mathrm{}`$. From Fig. 2, it is evident that the crossover in $`W_2(L)`$ for the LLS is due to the contribution from $`W_S`$. In the following we show that, for the LLS, the special form of the nonlinear steady state profile $`h_i(t)`$, which results in $`W_SL^{0.66}`$, is a consequence of the singular diffusion associated with the SOC state.

In , using a simple model, Carlson et. al. demonstrated that the SOC state of the system is associated with the vanishing of the density of troughs concommittant with the divergence of the corresponding diffusion coefficient. The troughs are defined as the sites for which $`z_i0`$ so that, in the LLS, an avalanche necessarily stops at a trough (or leave the system at $`i=0`$). On a coarse grained level we define a set of densities $`\{\rho _n(x,t);n=\mathrm{},\mathrm{},1,0,1,2\}`$ where $`\rho _n(x,t)`$ denotes the local density of sites with $`z_i=n`$. It follows that the coarse grained local slope $`z(x,t)`$, which is locally conserved by the dynamics, may be expressed in terms of the $`\rho _n`$’s as $`z(x,t)=_{n=\mathrm{}}^2n\rho _n(x,t)`$. The open boundary condition at $`i=0`$ implies $`h(0,t)=0`$ and the closed boundary condition at $`i=L`$ is modeled by setting $`z(L,t)=0`$. The trough density $`\rho (x,t)_{n=\mathrm{}}^0\rho _n(x,t)`$ is not strictly conserved. In the ’01’ model considered in  the slope has only two values $`z_i=0`$ (trough),$`1`$ and hence both $`z`$ as well as $`\rho `$ are strictly conserved and are related simply as $`z=1\rho `$. Although $`\rho `$ is not strictly conserved for LLS and ILLS we still approximate its dynamical evolution by the continuity equation, $`_t\rho (x,t)+_xJ(\rho (x,t))=0`$. Phenomenologically, the trough current $`J`$ consists of three parts: (i) current due to addition of grains $`J_0`$, (ii) avalanche current $`J_A`$, and (iii) a microscopic noise term $`\eta (x,t)`$. In analogy with driven diffusive systems  the phenomenological form of the avalanche current is written as $`J_A(\rho )=a\rho +b\rho ^2+\mathrm{}D(\rho )_x\rho `$. As in the 01 model, the critical state of the system ($`L\mathrm{},r=\mathrm{}`$) is associated with $`\rho ,_x\rho 0`$ and hence in order to balance the finite input flux $`J_0`$, the diffusion constant $`D(\rho 0)`$ has to diverge appropriately . Without loss of generality the leading divergence in $`D`$ is taken to be a simple pole and the evolution of $`\rho (x,t)`$ is thus ,

$$\frac{\rho }{t}=\frac{}{x}\left[D(\rho )\frac{\rho }{x}+\eta \right];D(\rho )=\frac{A(\rho )}{\rho ^\varphi }.$$
(7)
Here $`A(\rho )`$ is taken to be a smoothly varying function in the relevant interval $`0\rho 1`$.

In order to compute the steady state density profile $`\rho (x)=\rho (x,t\mathrm{})`$ of troughs, we note that, as $`L\mathrm{}`$, the system is arbitrarily close to criticality ($`\rho 0`$), and hence $`A(\rho )`$ may be set to a constant $`A_0=A(0)`$. One can then readily integrate the current equation $`D(\rho )d\rho /dx=J_0`$, subjected to the boundary condition $`\rho (L)=\rho _L`$, to obtain

$$\rho (x)\rho _L\left[1+\gamma (Lx)\right]^\theta $$
(8)
where $`\theta =1/(\varphi 1)`$ and $`\gamma =J_0\rho _L^{1/\theta }/(\theta A_0)`$. The average density of the troughs $`\overline{\rho }L^1_0^L\rho (x)𝑑x`$ thus scales with system size asymptotically as $`\overline{\rho }\rho _L(\gamma L)^\theta /(1\theta )`$. From Fig. 3 (inset), we obtain $`\theta 0.33`$ and thus $`\varphi 4`$ for the LLS, which are the same as the corresponding values obtained in  by direct measurement of $`D(\rho )`$ in a closed system.

In the following we find an approximate numerical relation between $`\rho (x)`$ to $`z(x)`$ which would enable us to find the steady state interface profile $`h_0(x)=_0^xz(x)𝑑x`$ form (8). In Fig. 3 (inset), we plot the average densities $`\overline{\rho }_nL^1_0^L𝑑x\rho _n(x,t)`$ in the steady state as functions of system size $`L`$. We note that, as $`L`$ increases, while $`\overline{\rho }_1\overline{\rho }_21/2`$ remain finite, densities of all the trough sites vanish algebraically with $`L`$: $`\overline{\rho }_0,\overline{\rho }_1L^{0.33},\overline{\rho }_2L^{0.6}`$, and $`\overline{\rho }_n`$’s with $`n3`$ are negligible. Hence, in the limit of large $`L`$, we can approximate the total trough density as $`\overline{\rho }\overline{\rho }_0+\overline{\rho }_1`$. From the normalization condition $`_{n=\mathrm{}}^2\rho _n(x,t)=1`$ it follows that $`\rho _1+\rho _2=1\rho `$ and we find numerically that $`\overline{\rho }_1\overline{\rho }_20.46\overline{\rho }`$. Thus, for large systems, we may write $`z(x)2\rho _2+\rho _1\rho _13/2\kappa \rho (x)`$, where numerically $`\kappa 2`$ . Thus, $`h_0(x)`$ is given approximately by

$$h_0(x)\frac{3}{2}x\frac{\kappa \rho _L(1+\gamma L)^{1\theta }}{\gamma (1\theta )}\left(1\left[\frac{1+\gamma (Lx)}{1+\gamma L}\right]^{1\theta }\right).$$
(9)
In Fig. 3, using the approximate expression for $`\rho (x)`$ from (8) with $`\varphi =4`$, we compare $`h_0(x)`$ with that obtained numerically and notice that the agreement is rather good given the nature of approximations involved .

The $`L`$ dependence of $`W_S`$ of the steady state profile can be computed from (9) and turns out to be

$$W_S(L)L^{1\theta }.$$
(10)
For the LLS with $`\theta 0.33`$, $`W_S(L)L^{0.66}`$ dominates $`W_1L^{1/3}`$ and $`W_cL^{0.18}`$ in (5) and hence, for large $`L`$, $`W_2(L)L^{0.66}`$ – in very good agreement with our numerical results (Fig. 2). It is interesting to note that $`W_1`$ and $`W_2`$ would have different asymptotic behaviours if $`\theta <2/3`$, i.e., if $`\varphi >5/2`$. For the ‘01’ model studied in  it is shown that $`\varphi =3`$ exactly and thus $`W_2(L)L^{1/2}`$ for the corresponding interface . Our conjecture that $`W_1`$ and $`W_2`$ would have different asymptotic behaviour for SOC systems would be invalid if one can devise a model with $`\varphi <5/2`$.

Next, we present a general argument as to why, for noncritical models such as the ILLS, we do not expect any ambiguity in roughness exponent. The essential difference between a critical and a noncritical model is the presence of an additional length scale (apart from system size and the microscopic cutoff) which shows up in the avalanche size distribution (e.g. Fig. 1 for ILLS) as well as governs the decay of the boundary effects into the bulk. In the ILLS this length scale is related to the parameter $`r`$. To see this we first note that as a cluster moves down its momentum $`p`$ makes a random walk . If the mean density of troughs $`\rho <1/r`$ then most avalanches leave the system resulting in net drainage and if $`\rho >1/r`$ then most avalanches stop on the pile leading to net growth of the pile. Thus, in the steady state one has a finite density of troughs $`\rho =1/r`$ and thus the mean spacing between the troughs $`1/\rho =r`$ appears as the additional length scale. As already discussed above, the current $`J(\rho )`$ now has a nonzero systematic part $`a\rho +b\rho ^2\mathrm{}`$, in addition to the finite diffusion term. Such a systematic current, e.g., a term such as $`a\rho `$, results in the effects of the boundaries decaying exponentially inside the bulk. This is in contrast to the long ranged power law decay in the LLS which is reflected in the forms of $`\rho (x)`$ and $`h_0(x)`$. Although the steady state surface is nonlinear (Fig. 3), its rms wandering $`W_S`$ is bounded and does not alter the asymptotic scaling of $`W_2`$.

Lastly, we briefly point out possible pitfalls of using $`W_2`$ naively without properly accounting for the underlying steady state profile, even in systems which do not show SOC. Recently, Barker and Mehta  studied a disordered version of the LLS where disorder was introduced by allowing the added grains to have an aspect ratio different from unity: grains are now rectangles and are deposited on the sandpile with fixed probabilities of landing on their larger or smaller edges. Thus, the height $`h_i`$ of column $`i`$, no longer an integer, is the sum of the vertical edges of all the grains in that column. In addition to the threshold dynamics of LLS, dynamical reorientation of cluster of grains were allowed. They found that the larger sandpiles cease to display SOC which is reflected in the emergence of a preferred size of the large avalanches in the associated drop number distributions. Their numerical studies of $`W_2`$ showed that while for very small systems ($`L100`$) the roughness exponent $`\alpha 0.34`$, it seems to crossover to a much larger value $`\alpha 0.72`$ for larger system sizes, $`100L400`$ (Fig. 4). Our preliminary numerical studies, using the same parameters as in , of yet larger systems ($`500L2048`$) show that in fact the crossover seen in  is evidently transient and the asymptotic roughness exponent goes back to $`\alpha 0.33`$ (see Fig. 4).

In fact, the crossover reported in  is less noticeable if one looks at $`W=\sqrt{W_2^2W_S^2}`$ instead of $`W_2`$ (Fig. 4) . This is in accordance with our conjecture since the disordered system is not critical, although the steady state profile is nonlinear (Fig. 4, inset), it does not affect the asymptotic roughness measured by $`W_2`$. The transient crossover seen in $`W_2`$ here is similar to what is seen for the ILLS in Fig. 2, inset.

In summary, we have studied the surface roughness of a prototype model of SOC and its modifications in one dimension. We find that one needs to be careful in defining quantities such as the interface width since special form of the steady state shape of the surface in systems with SOC can result in different asymptotic behaviours of otherwise equivalent definitions. Although there is no such ambiguity for noncritical models, still there may be crossovers at intermediate length scales which should not be taken as the true asymptotic behaviour.

JGO and GT acknowledge financial support of Centro de Física do Porto, University of Porto where part of this work was carried out. JGO also acknowledges financial support of Universidade de Aveiro. JFFM was partially supported by project POCTI/1999/FIS/33141 and POCTI/2002/MAT/46176.


text:
# A Unified Jet Model of X-Ray Flashes and Gamma-Ray Bursts

## 1 Introduction

Two-thirds of all HETE-2–localized bursts are either “X-ray-rich” or X-Ray Flashes (XRFs); of these, one-third are XRFs 111We define “X-ray-rich” GRBs and XRFs as those events for which $`\mathrm{log}[S_X(230\mathrm{kev})/S_\gamma (30400\mathrm{kev})]>0.5`$ and 0.0, respectively. (Sakamoto et al., 2003b). These events have received increasing attention in the past several years (Heise et al., 2000; Kippen et al., 2002), but their nature remains unknown.

XRFs have $`t_{90}`$ durations between 10 and 200 sec and their sky distribution is consistent with isotropy. In these respects, XRFs are similar to “classical” GRBs. A joint analysis of WFC/BATSE spectral data showed that the low-energy and high-energy photon indices of XRFs are $`1`$ and $`2.5`$, respectively, which are similar to those of GRBs, but that the XRFs had spectral peak energies $`E_{\mathrm{peak}}^{\mathrm{obs}}`$ that were much lower than those of GRBs (Kippen et al., 2002). The only difference between XRFs and GRBs therefore appears to be that XRFs have lower $`E_{\mathrm{peak}}^{\mathrm{obs}}`$ values. It has therefore been suggested that XRFs might represent an extension of the GRB population to bursts with low peak energies.

Clarifying the nature of XRFs and X-ray-rich GRBs, and their connection to GRBs, could provide a breakthrough in our understanding of the prompt emission of GRBs. Analyzing 42 X-ray-rich GRBs and XRFs seen by FREGATE and/or the WXM instruments on HETE-2, Sakamoto et al. (2003b) find that the XRFs, the X-ray-rich GRBs, and GRBs form a continuum in the \[$`S_\gamma (2400\mathrm{kev}),E_{\mathrm{peak}}^{\mathrm{obs}}`$\]-plane (see Figure 1, left-hand panel). This result strongly suggests that all of these events are the same phenomenon.

Furthermore, Lamb et al. (2003) have placed 9 HETE-2 GRBs with known redshifts and 2 XRFs with known redshifts or strong redshift constraints in the ($`E_{\mathrm{iso}},E_{\mathrm{peak}}`$)-plane (see Figure 1, right-hand panel). Here $`E_{\mathrm{iso}}`$ is the isotropic-equivalent burst energy and $`E_{\mathrm{peak}}`$ is the energy of the peak of the burst spectrum, measured in the source frame. The HETE-2 bursts confirm the relation between $`E_{\mathrm{iso}}`$ and $`E_{\mathrm{peak}}`$ found by Amati et al. (2002) (see also Lloyd, Petrosian & Mallozzi (2000)) for GRBs and extend it down in $`E_{\mathrm{iso}}`$ by a factor of 300. The fact that XRF 020903, one of the softest events localized by HETE-2 to date, and XRF 030723, the most recent XRF localized by HETE-2, lie squarely on this relation (Sakamoto et al., 2003a; Lamb et al., 2003) provides strong evidence that XRFs and GRBs are the same phenomenon. However, additional redshift determinations are clearly needed for XRFs with 1 keV $`<E_{\mathrm{peak}}<30`$ keV energy in order to confirm these results.

Figure 2 shows a simulation of the expected distribution of bursts in the ($`E_{\mathrm{iso}},E_{\mathrm{peak}}`$)-plane (left panel) and in the ($`F_N^{\mathrm{peak}},E_{\mathrm{peak}}`$)-plane (right panel), assuming that the (Amati et al., 2002) relation holds for XRFs as well as for GRBs (Lamb, Donaghy & Graziani, 2003), as is strongly suggested by the HETE-2 results. The SXC, WXM, and FREGATE instruments on HETE-2 have thresholds of $`16`$ keV and considerable effective areas in the X-ray energy range. Thus HETE-2 is ideally suited for detecting and studying XRFs. In contrast, BAT on Swift has a nominal threshold of 20 keV. This simulation shows that the WXM and SXC instruments on HETE-2 detect many times more bursts with $`E_{\mathrm{peak}}<10`$ keV than will BAT on Swift.

## 2 XRFs as a Probe of GRB Jet Structure, GRB Rate, and Core Collapse Supernovae

Most GRBs have a “standard” energy Frail et al. (2001); Panaitescu & Kumar (2001); Bloom, Frail & Kulkarni (2003); i.e, if their isotropic equivalent energy is corrected for the jet opening angle inferred from the jet break time, most GRBs have the same radiated energy, $`E_\gamma =1.3\times 10^{51}`$ ergs, to within a factor of $``$ 2-3.

Two models of GRB jets have received widespread attention:

* The “structured jet” model (see the left-hand panel of Figure 3). In this model, all GRBs produce jets with the same structure (Rossi, Lazzati, & Rees, 2002; Woosley, Zhang, & Heger, 2003; Zhang & Mészáros, 2002; Mészáros, Ramirez-Ruiz, Rees, & Zhang, 2002). The isotropic-equivalent energy and luminosity is assumed to decrease as the viewing angle $`\theta _v`$ as measured from the jet axis increases. The wide range in values of $`E_{\mathrm{iso}}`$ is attributed to differences in the viewing angle $`\theta _v`$. In order to recover the “standard energy” result (Frail et al., 2001), $`E_{\mathrm{iso}}(\theta _v)\theta _v^2`$ is required (Zhang & Mészáros, 2002).
* The “uniform jet” model (see the right-hand panel of Figure 3). In this model GRBs produce jets with very different jet opening angles $`\theta _{\mathrm{jet}}`$. For $`\theta <\theta _{\mathrm{jet}}`$, $`E_{\mathrm{iso}}(\theta _v)`$ = constant while for $`\theta >\theta _{\mathrm{jet}}`$, $`E_{\mathrm{iso}}(\theta _v)=0`$.

As we have seen, HETE-2 has provided strong evidence that the properties of XRFs, X-ray-rich GRBs, and GRBs form a continuum, and that these bursts are therefore the same phenomenon. If this is true, it immediately implies that the $`E_\gamma `$ inferred by (Frail et al., 2001) is too large by a factor of at least 100 (Lamb, Donaghy & Graziani, 2003). The reason is that the values of $`E_{\mathrm{iso}}`$ for XRF 020903 (Sakamoto et al., 2003a) and XRF 030723 (Lamb et al., 2003) are $``$ 100 times smaller than the value of $`E_\gamma `$ inferred by Frail et al. – an impossibility.

HETE-2 has also provided strong evidence that, in going from XRFs to GRBs, $`E_{\mathrm{iso}}`$ changes by a factor $`10^5`$ (see Figure 1, right-hand panel). If one tries to explain only the range in $`E_{\mathrm{iso}}`$ corresponding to GRBs, both the uniform jet model and the structured jet model work reasonably well. However, if one tries to explain the range in $`E_{\mathrm{iso}}`$ of a factor $`10^5`$ that is required in order to accommodate both XRFs and GRBs in a unified description, the uniform jet works reasonably well while the structured jet model does not.

The reason is the following: the observational implications of the structured jet model and the uniform jet model differ dramatically if they are required to explain XRFs and GRBs in a unified picture. In the structured jet model, most viewing angles $`\theta _v`$ are $`90^{}`$. This implies that the number of XRFs should exceed the number of GRBs by many orders of magnitude, something that HETE-2 does not observe (see Figures 1, 2, 4, and 5). On the other hand, by choosing $`N(\mathrm{\Omega }_{\mathrm{jet}})\mathrm{\Omega }_{\mathrm{jet}}^2`$, the uniform jet model predicts equal numbers of bursts per logarithmic decade in $`E_{\mathrm{iso}}`$ (and $`S_E`$), which is exactly what HETE-2 sees (again, see Figures 1, 2, 4, and 5) (Lamb, Donaghy & Graziani, 2003). Thus, if $`E_{\mathrm{iso}}`$ spans a range $`10^5`$, as the HETE-2 results strongly suggest, the uniform jet model can provide a unified picture of both XRFs and GRBs, whereas the structured jet model cannot. This means that XRFs provide a powerful probe of GRB jet structure.

A range in $`E_{\mathrm{iso}}`$ of $`10^5`$, which is what the HETE-2 results strongly suggest, requires a minimum range in $`\mathrm{\Delta }\mathrm{\Omega }_{\mathrm{jet}}`$ of $`10^410^5`$ in the uniform jet model. Thus the unified picture of XRFs and GRBs in the uniform jet model implies that there are $`10^410^5`$ more bursts with very small $`\mathrm{\Omega }_{\mathrm{jet}}`$’s for every such burst we see; i.e., the rate of GRBs may be $`100`$ times greater than has been thought.

In addition, since the observed ratio of the rate of Type Ic SNe to the rate of GRBs in the observable universe is $`R_{\mathrm{Type}\mathrm{Ic}}/R_{\mathrm{GRB}}10^5`$ (Lamb, 1999), a unified picture of XRFs and GRBs in the uniform jet model implies that the GRB rate is comparable to that of Type Ic SNe (Lamb, Donaghy & Graziani, 2003). More spherically symmetric jets yield XRFs and narrow jets produce GRBs. Thus XRFs and GRBs provide a combination of GRB/SN samples that would enable astronomers to study the relationship between the degree of jet-like behavior of the GRB and the properties of the supernova (brightness, polarization $``$ asphericity of the explosion, velocity of the explosion $``$ kinetic energy of the explosion, etc.). GRBs may therefore provide a unique laboratory for understanding Type Ic core collapse supernovae.

A unified picture of XRFs and GRBs in the uniform jet model also implies that many Type Ic SNe produce narrow jets, which may suggest that the collapsing cores of many Type Ic supernovae are rapidly rotating. Finally, such a unified picture implies that the total radiated energy in gamma rays $`E_\gamma `$ is $``$ 100 times smaller than has been thought Lamb, Donaghy & Graziani (2003).


text:
# Neutrino flux variations and solar activity

## 1 Introduction

A noticeable progress has been achieved in the so-called ”solar neutrino problem” during the last decade. This was caused, first of all, by the flux measurements taken with new neutrino experiments like Kamiokande, Super Kamiokande, GALLEX and SAGE. These experiments (together with the Homestake experiment) have, in particular, different energy thresholds that makes it possible to study different regions of the solar neutrino spectrum. Their results have confirmed that the observed neutrino flux is by a factor of 2-3 smaller than that calculated from the standard solar model (Suzuki Suzuki (1998), Abdurashitov Abdurashitov (2002)). The analysis of recent observational data from the Super Kamiokande and Sudbury Neutrino Observatory has shown that some fraction of electron neutrinos forming in the process of thermonuclear reactions in the Sun’s core are transformed on their transit to the Earth into $`\mu `$ and $`\tau `$ neutrinos which cannot be detected by radiochemical experiments (Mikheev & Smirnov Mikheev (1987)). In this concern, another problem has drawn more attention: whether the variations of the solar neutrino flux are real and correlate with solar activity.

In a series of investigations, made mainly with use of data of Homestake experiment, it was found that the flux of the solar neutrino varies in time and these variations are related with solar activity. (Subramanian Subramanian (1979), Sakurai Sakurai (1980), Bazilevskaya et al. Bazilevskaya (1982), Bieber et al. Bieber (1990), Gavryusev & Gavryuseva Gavryusev91 (1991, 1994), Massetti & Storini Massetti93 (1993, 1996), Oakley et al. Oakley (1994), Snodgrass & Oakley Snodgrass (1999), Rivin & Obridko Rivin (1997), Ikhsanov & Miletsky Ikhsanov99 (1999, 2000, 2002)). However, some authors regard the variations as insignificant. (Bahcall Bahcall (1989), Fukuda et al. Fukuda (1996), Walther Walther97 (1997, 1999), Cattaneo Cattaneo02 (2002)). Recently Sturrock et al. (Sturrock & Scargle Sturrock01 (2001), Sturrock & Weber Sturrock02 (2002)) basing on data of relatively short series of GALLEX-GNO and SAGE experiments found that the solar neutrino flux varies on time scale comparable with weeks. Besides, Milsztajn (Milsztajn Milsztajn (2003)) states that a period of 13.75 days is observed in Kamiokande data. The latter statement, however, is negated in paper of (Yoo et al. Yoo (2003)). Therefore, the question about variability of the solar neutrino flux is still open, especially as regards variations with period longer than one year, which is studied in the present work.

Searching for time variations caused by the solar activity cycle, most investigations focused on the 11- and 2-year periodicity. Our previous paper (Ikhsanov & Miletsky Ikhsanov99 (1999)) has shown that in the region of low frequencies a periodicity close to 5 years (hereafter, 5-year periodicity) plays the most important role in the variations of the neutrino flux. It is, therefore, not sufficient to analyze a connection between these variations and indices of solar activity, displaying predominantly 11-year period. Account for 5-year periodicity essentially changes our views on a probable character of interconnection between variations of the neutrino flux and solar activity.

We have made an attempt to find an index of solar activity variations of which show pronounced manifestation of the 5-year periodicity (Ikhsanov & Miletsky Ikhsanov99 (1999)). Among many solar indices under consideration, the number of polar coronal holes (PCH) is the best to meet this requirement. It is worthwhile to note that, in spite of high enough correlation coefficient between this index and the neutrino flux (0.68), it has no direct relation to the regions with large magnetic fields on the solar surface. Furthermore, the polar coronal holes are not situated on the way of neutrinos between the solar core and the Earth.

In the present paper we give further justification of necessity to account for 5-year periodicity in the variations of the solar neutrino flux. In order to clarify a connection between the neutrino flux and solar activity, we investigate, besides the activity indices relating to the solar surface (such as Wolf sunspot numbers, sunspot areas, PCH and some characteristics of the global magnetic fields), those indices, which characterize deeper solar layers (p-modes and solar radius), as well as the near-Earth parameters (the concentration of particles and cosmic rays).

## 2 Periodicities in variations of the solar neutrino flux

We use the data, obtained in 1970-1994 with the Cl-Ar detector Homestake (Cleveland et al.Cleveland (1998)). The neutrino flux was measured by detecting the number of 37Ar atoms generated in the reaction $`\nu `$+37Cl $``$ 37Ar+e-, with the counting rate being expressed in atoms per day. A series of measurements taken with the Homestake detector consists of consecutive runs of unequal duration and is contaminated with significant uncertainties. The duration of the runs had been changing around 1.2-2 months. The time series has a number of small gaps and one very continuous ($``$ 1.6 years). That is why we applied successively a number of methods to study the time series of the neutrino flux measurements.

We have analyzed two series of the neutrino flux measurements. The first one (hereafter, series I) contains the measurements of the neutrino flux assigned to the middle of the corresponding run. On the basis of this series, we reconstructed the second one (hereafter, series II), consisting of monthly averaged neutrino flux estimates. The series I contains 108 (the number of runs) unevenly spaced flux estimates, while the second one consists of 283 monthly means from October 1970 until April 1994 (without breaks). It should be noted that we have constructed the series II by merely practical reasons since the mathematical methods of data processing are mainly applied to the evenly spaced data. As will be shown below, the result is practically the same if only those periodicities are considered which essentially exceed the intervals between observations (runs).

In order to reveal peculiarities of temporal variations within these series, we calculated, at the first stage, the estimates of the power density spectrum (PDS), using a periodogram analysis. It is known (see, e.g., (Otnes & Enochson Otnes (1978), Marple Marple (1987)) that in the spectral analysis of such data spurious peaks and frequency shifts can arise. To avoid these disadvantages, a time series is usually transformed using a data-weighting window. We applied this procedure, having chosen a Blackman-Harris (BH) window (Otnes & Enochson Otnes (1978)). Then, after subtraction of the mean values for each series, we calculated the normalized periodogram estimates of PDS. Since the series I contains unevenly spaced data, we have adopted the Lomb-Scargle techniques for their analysis(Lomb Lomb (1976), Scargle Scargle (1982)). In order to compute the values of PDS for the series II, we applied a well-known algorithm (Otnes & Enochson Otnes (1978)) based on the Fast Fourier Transform (FFT) and allowing to estimate statistical significance of the obtained peaks.

Figure 1a presents the PDS for the series I (without using BH window). Comparison with PDS obtained with the use of BH window (Fig. 1b) shows essential differences in the region of 2-year periodicity. In the first case, one can see a set of peaks in the interval 1-3 years, with the maximum peak corresponding to the period 2.1 years (see Fig.1a), while in the second case (Fig.1b) their power is essentially smaller. The same picture is seen in the PDS for the series II (Fig.1c), with a high and stable peak, corresponding to the period 4.5-4.7 years, being seen in all three figures. This proves, in particular, that the series II can be used to reveal the periodicities in the low-frequency spectral region.

Tables I and II presents the values of significance levels $`\alpha `$ for the PDS peaks under consideration, determined with the use of shuffle test, randomly reassigning measurements among runs for the series I (Sturrock et al. Sturrock (1998)), with 10000 test spectra having been computed for each estimate. In the first case (Fig. 1a), two peaks (4.6 and 2.1 years) are significant at the level $`\alpha `$ $`<`$ 4%, which corresponds to the level of confidence probability P $`>`$ 96% (P=1-$`\alpha `$), while in the second case (Fig. 1b) - only one peak, which represents the 5-year periodicity (P=97%). Besides that, the peak corresponding to the period of 1.3 years can be noticed. Let us note that for the peak corresponding to the 11-year period, the level of confident probability is only 90%.

We have checked, to what extent the manifestations of the most significant periodicities, revealed in our analysis, depend on the time interval, covered by the series under investigation. For this purpose we calculated the PDS for two time intervals 1970-1982 and 1983-1994, containing data of 54 runs each. In the spectrum calculated over the first interval and shown in Fig. 1c, there is a noticeable peak near 2 years, while this peak is absent in the spectrum, computed for the second interval. This explains the fact that a number of authors have got significant amplitudes of 2-year oscillations, while the PDS, calculated over the whole series, does not contain significant signs of this periodicity. For the period around 5 years, a noticeable peak has been found in the spectra for both intervals, that indicates somewhat higher stability of the latter periodicity in comparison with 2-year oscillations.

On the above mentioned grounds, we have performed a more detailed study of the neutrino flux variations in the region of low frequencies by means of two different methods. The first method implied running-mean smoothing of the neutrino flux (Fig.2a), using a 15-month window with harmonic weights. A smoothed curve (a solid line in Fig.2b) clearly displays the presence of 5-year oscillations. The second method supposed implementation of the Butterworth filter (Otnes & Enochson Otnes (1978)) with the frequency interval of maximum transmission, corresponding to the period interval 4-14 years. The results are presented in Fig.2c. Both independent methods are in a good agreement regarding the shape of the 5-year wave. As seen in Fig.2b,c the average distance between maxima is 4.6 $`\pm `$ 0.7 years. Some uncertainty is connected with the gap in observations in 1985 – 1986 (1.6 year). The case under consideration corresponds to the linear interpolation of the flux values in this interval. Fig.2a presents another possible extreme variant of the flux behavior (shown by dotted line). In this case, the maximum flux value in the region of interest is somewhat lower, and its position is shifted to the right, that, however, has not changed a general picture of 5-year oscillations. In the next section will be shown that the first case seems to be more realistic. The appearance of six consecutive 5-year periods (taking into consideration the GALLEX data) leaves no doubts about the reality of these oscillations. Fig.6a demonstrates variability of the neutrino flux over the interval 1-3 years. It indicates that a periodicity of 2 years was in operation before 1980, but had afterwards disappeared and arose again after 1991.

A lot of attempts have been made in the recent years to discard the question about the neutrino flux variations because of large errors in the radiochemical measurements. However, this question is very important in the study of the Sun’s nature, especially with regard to the problems of solar magnetism, since the presence of variations in the solar neutrino flux (due to either the precession of the neutrino spin (VVO-effect) or the resonance spin-flavor precession (RSFP) (Akhmedov Akhmedov89 (1989), Akhmedov97 (1997)) could make it possible to investigate the magnetic field in the interiors of the Sun. Cattaneo (Cattaneo02 (2002)) has recently criticized the data from the Homestake experiment for insufficient background reduction. The main argument of the author was that the runs with large background value should be removed from further analysis because of large uncertainties. Thus he has increased the average production of the neutrino flux up to 0.566 at/day after removing one third of the runs with large background. In principle, this argument seems reasonable, but the question is whether it can influence the character of the neutrino flux variations revealed above. In order to clarify this question, we removed one third of the runs with large background value ( 0.027) from the Homestake data and calculated the power spectrum for the rest of the data. The result is shown in Fig. 3c. The appearance of the spectrum is even better than that of the complete series (Fig. 3a). Really, the periodicities of 10, 4.5 and 2.5 years are present at the same significance level (Table I). The peak corresponding to 2.5 years has increased while all the others decreased. In addition, the same figure shows the PDS of the background values for both the complete Homestake data and the series after reduction. A weak peak can be noticed at 8 years in the first case (Fig. 3b), and at 9 years in the second case (Fig. 3d). At the same time, there are no indications on the presence of 2- and 5-years periodicities in the PS of the background values. Thus, the above analysis of the neutrino flux variations in the time interval 1970-1994 shows that 5-year periodicity represents its most stable oscillation.

## 3 Connection between the neutrino flux and solar activity

As mentioned in the Introduction, variations of the neutrino flux can be caused by periodic changes of solar activity. The Wolf numbers (W) are the most frequently used to search for this interrelation. Together with the index of sunspots areas, they characterize the toroidal component of the Sun’s magnetic field. Fig.2e presents the coefficients of running correlation (using the running interval of 40 months with a step of 1 month), between the series II of the neutrino flux (Fig.2a) and the W series (Fig.2d). The correlation coefficient (CC) varies essentially between -0.8 and +0.6, with its changing from -0.77 up to +0.57 over the time interval 1987-1991. This fact explains the conclusion about an absence of significant correlation during this time interval, made by a number of authors.

The average value of the CC calculated over the whole series shown in Fig.2e is equal to -0.19. If we calculate the CC between the smoothed series of neutrino fluxes II and W presented in Fig.1b and 1d, respectively, we get the value -0.49. The same or close value has been obtained by other investigators (see, e.g., (Rivin Rivin (1997))). So week correlation can be connected with the fact that during 11-year cycle of solar activity we observe two maxima in the neutrino flux series, while in the Wolf numbers - only one (Fig.2b,d).

In the previous paper (Ikhsanov & Miletsky Ikhsanov99 (1999)) we have made an attempt to find an index of solar activity, in variations of which the 5-year periodicity were manifested the most clearly. This requirement turns out to be met by the PCH number. However, they are not situated on the transit from the Sun’s core to the Earth, as required within the mechanism suggested by Voloshin, Visotsky & Okun (Voloshin (1986)) (hereafter VVO).

It is, therefore, more interesting to consider the temporal changes of the Sun’s magnetic field on the transit of neutrinos to the Earth. If we take the radii of the Sun’s core regions, where 8B and pp-neutrinos are generated, as $`310^9cm`$ and $`1.810^{10}cm`$, respectively, then their projections on the Sun’s surface can be presented as the circles with diameter of $``$50 and 150. On the basis of Stanford data (1976-1998), we constructed the series of values of the magnetic field B and its module $`\left|B\right|`$ for both cases with account for annual change of the ecliptic’s inclination with respect to the plane of the Sun’s equator $`\left(\pm 7^015^{}\right)`$. As seen in Fig.4, the PDS of these two series are completely different. The PDS of the B series shows 22-, 5-, and 3-year periodicities, while that of the $`\left|B\right|`$ series displays only 10-year period. It is worthwhile to note that the PDS for 50\- and 150-circles differ only in some shift of the peaks, corresponding to the 3- and 5-year periods.

Figs. 5b,c present the curves for these series after removing all the components beyond the period interval 4-14 years (using a Butterworth filter). Comparison between the series obtained by means of the same transformation of the data on magnetic field B and neutrino fluxes reveals a low degree of similarity (CC=0.45). The correlation between the neutrino flux and $`\left|B\right|`$ series is also rather low, with the CC being negative (CC=-0.43), as in the case of W-series. A higher CC can be obtained by means of further smoothing, but some essential information would be lost.

Fig.5a presents the neutrino series (a dotted line), based on the data from Kamiokande (Gavryusev & Gavryuseva Gavryusev94 (1994)) and GALLEX (Kirsten Kirsten (1999)), with distinguishing the period interval of 4-14 years by means of the same method. The curve for the first of these series (covering a short time interval) resembles, in spite of a small phase shift, that for the Homestake measurements. Regarding the Gallex series, it demonstrates surprisingly strong similarity with the Homestake data for the period interval under consideration.

The same figure presents the plots for a few more series of the solar characteristics, mentioned in the Introduction and transformed using the same techniques. For the analysis of temporal variations of the solar radius (RA), we have chosen the data obtained in the experiment CERGA (Laclare et al.Laclare (1996)), which covers, unfortunately, a shorter period (1978-1994), than the neutrino flux estimates. The advantage of these measurements is that all of them had been carried out by the same observer. The frequency shifts of p-modes, obtained at the Observatorio del Teide (1980-1994), have been taken from (Regulo Regulo (1994)). All these data, together with those on the number of PCH (Ikhsanov & IvanovIkhsanov99a (1999)), the intensity of cosmic rays (CR), concentration of interplanetary particles (ION) and on the integral value of the sunspots areas (SA) have been used to form the series of monthly mean values. They are shown in Fig.5 (after filtration using the same procedure, as for the neutrino data). The highest CC=0.82 has been obtained between the neutrino flux measurements (Fig.5a) and the RA observations (Fig.5d). Somewhat lower correlation of the neutrino data (CC=0.62) has been found with the PCH numbers (fig.5e) and ION (Fig.5f), and essentially lower with the series PM (Fig.5g), CR and SA (Fig.5h), with the CC for the latter three series being negative.

Because of a small coverage of the radius observations, the value of its CC is not confident enough. A possible shift (dashed line) of the maximum position in the neutrino flux series in 1985-1986 leads to the essential decrease of the CC from 0.83 down to 0.33. The same can be seen for PM and SA, that is caused by a strong sensitivity of the CC to a displacement of even one maximum in case of short time series. However, examining Fig 5a, one can see that the maxima of the neutrino series at the phases of minima of the 11-year cycle are noticeably higher, than the neighboring maxima. This regularity as well as some other arguments favors the conclusion, that real variations of the neutrino flux in 1985-1986 are better represented by a solid line in Fig 5a.

Summarizing, we can infer that the neutrino flux series has a positive CC with those solar activity indices, which display a pronounced 5-year periodicity, with its value being in the interval 0.6-0.8 (except for the B series). As far as the indices SA, $`\left|B\right|`$, PM and CR are concerned, which show only 11-year period in the region of low frequencies, their CC with the neutrino data turn out to have a negative value in the interval between -0.4 and -0.5.

Fig.6 presents the curves for the series under investigation after removing (using the same procedure) all the components beyond the period interval 1-3 years. All the series demonstrate some changes in the periods and amplitudes of oscillations. As mentioned above, a periodicity close to 2 years in the neutrino flux measurements (Fig.6a) was in operation during 1970-1980, then it had disappeared, and was detected again only after 1992. The most clearly this 2-year periodicity can be seen in the RA series, but only until 1990, that is just when corresponding oscillations in the neutrino data were weak. The series of the SA and B measurements demonstrate 2-year oscillations, as a rule, near the maxima of 11-year cycles. Thus, all these indices show either very weak correlation with the neutrino flux data, or its absence, since the module of corresponding CC does not exceed 0.25 (except for RA, with CC = 0.31).

## 4 Conclusions

Thus, we have found, that the periodicity at approximately 5 years (4.6 $`\pm `$ 0.7 years) is the most stable in the variations of the neutrino flux. The 11-year cycle manifests itself only in the fact, that the maxima of these oscillations are somewhat higher at the phases of minima of the 11-year solar activity cycle. In 1970s and 1990s a 2-year periodicity was also observed. It is especially interesting to note that 2- and 5-year periods are present in the variations of the neutrino flux measured in both Homestake and GALLEX experiments. This fact supports the reality of these variations. The 5-year oscillations can be produced either in the interiors of the Sun (like the 11- and 2-year oscillations) or as a result of some near-Earth phenomena. If the former is true, the most probable seems the assertion that the 5-year oscillations arise on the transit of the electron neutrinos from the Sun’s core to the Earth. This, however, requires some conditions on the neutrino properties and interiors of the Sun to be satisfied (Voloshin Voloshin (1986), Akhmedov Akhmedov89 (1989), Akhmedov97 (1997)). Namely, within the VVO model, for the neutrino flux to be modulated, the neutrino magnetic moment should be significant as well as a cyclic toroidal component of the magnetic field, which should exist on their way inside the Sun.

As already mentioned in the Introduction, a number of authors have been looking for a connection between the solar neutrino flux and solar activity for more than twenty years. They concentrated, however, mainly on the attempts to find a correlation with 11-year cycle. The above analysis allows to conclude, that the series of solar activity indices displaying only 11-year periodicity (W, SA, $`\left|B\right|`$, PM and CR) show only weak correlation with the neutrino flux ($`\left|CC\right|`$ = 0.5), if not too smoothed data are analyzed. As shown above, this approach does not provide an appropriate choice of the indices, since it does not take into account an existence of oscillations at 4-5 years in the neutrino counting rates. The latter period manifests itself in the series of RA, PCH and ION indices. All of them correlate with the neutrino flux series with CC not less than 0.62. Let us note, however, that the data on RA and PM cover only short time intervals, and, therefore, some doubts may arise about the degree of their correlation with the neutrino time series. The first group of indices of solar activity, displaying 11-year periodicity, reflects mainly the properties of quadruple component of the Sun’s magnetic field. The second group (PCH, ION) relates to the dipole component of the global magnetic field (Ikhsanov & Miletsky Ikhsanov00 (2000)) and is determined by deep layers of the Sun. Finally, as seen in Fig.5a,b, near the minimum of 11-year cycle of solar activity, the magnetic field B is close to zero, while the neutrino flux reaches its maximum. This is in agreement with predictions of VVO model, that the maximum neutrino counting rates should be observed near the minimum of the toroidal magnetic field. Thus, there exists an essential interconnection between the neutrino flux and magnetic field of the Sun. In this case we can assume that a 5-year periodicity in the solar neutrino measurements can result from cyclic variations of the toroidal magnetic field, situated at the bottom of the convective zone or deeper.

###### Acknowledgements.

The authors are grateful to Dr. N. Beskrovnaya for her assistance in preparation of this manuscript.