text:
# Cosmological X-ray Flashes in the Off-Axis Jet Model

## 1 INTRODUCTION

The X-ray flash (XRF) is a class of X-ray transients (Heise et al. 2001, see also Barraud et al. 2003). Some properties of XRFs, such as the observed event rate, the duration and the isotropic distribution, are similar to that of Gamma-Ray Bursts (GRBs), while the spectral hardness of XRFs characterized by the peak flux ratio, the fluence ratio and the photon index is softer than that of GRBs. This class represents a large portion of the whole GRB population. Recently, the observational value of $`V/V_{\mathrm{max}}`$, which is the measure of the homogeneity of spatial distribution (Schmidt, Higdon, & Hueter 1988; see also Chang & Yi 2001; Kim, Chang, & Yi 2001), has been updated from $`0.56\pm 0.12`$ (J. Heise 2000, talk given in 2nd workshop Gamma-Ray Bursts in the Afterglow Era) to $`0.27\pm 0.16`$ (J. Heise 2002, talk given in 3rd workshop Gamma-Ray Bursts in the Afterglow Era). The updated value of $`V/V_{\mathrm{max}}`$ suggests that XRFs take place at a cosmological distance.

Various models accounting for the nature of the XRFs have been proposed (Yamazaki, Ioka, & Nakamura 2002; Heise et al. 2001; Dermer, Chiang, & Böttcher 1999; Huang, Dai, & Lu 2002; Mészáros et al. 2002; Mochkovitch et al. 2003; Daigne & Mochkovitch 2003). Heise et al. (2001) proposed that XRFs could be GRBs at high redshift. The redshifts of XRF011030 and XRF020427 have an upper limit of $`z3.5`$ (Bloom et al. 2003; J. Heise 2002, talk given in 3rd workshop Gamma-Ray Bursts in the Afterglow Era), which suggests that XRFs take place at not so high redshift but the same as that of GRBs. The photosphere-dominated fireball model may account the nature of the XRFs with peak energy $`E_p`$ more than $`20`$ keV (Mészáros et al. 2002, Ramirez-Ruiz & Lloyd-Ronning 2002). However, further considerations are needed to explain the event with $`E_p`$ of approximately a few keV, such as XRF020427 (L. Amati 2002, talk given in 3rd workshop Gamma-Ray Bursts in the Afterglow Era), XRF020903 and XRF010213 (N. Kawai 2002, talk given in 3rd workshop Gamma-Ray Bursts in the Afterglow Era). The models with small Lorentz factors, such as the dirty fireball model (Dermer, Chang, & Böttcher 1999; Huang, Dai, & Lu 2002) or the structured-jet model (Rossi, Lazzati, & Rees 2001; Woosley et al. 2002; Zhang & Mészáros 2002a) also have possibilities to explain the properties of the XRF, with the implicit assumption that the XRF arises not from internal shocks (Zhang & Mészáros 2002b). The models for internal shocks with small contrast of high Lorents factors might be the origin of XRFs (Mochkovitch et al. 2003; Daigne & Mochkovitch 2003).

For the other possibility, we have studied the off-axis jet model and proposed that if we observe the GRB jet with a large viewing angle, it looks like an XRF (Yamazaki, Ioka, & Nakamura 2002, Yamazaki, Ioka, & Nakamura 2003). In Yamazaki, Ioka, & Nakamura (2002) the value of the jet opening half-angle was adopted as $`\mathrm{\Delta }\theta =0.1`$. In this model the distance to the farthest XRF ever detected is about 2 Gpc ($`z0.4`$) so that the cosmological effect is small and $`V/V_{\mathrm{max}}0.5`$. Recent observations suggest that GRBs with relatively small opening angle exist, while the distribution of $`\mathrm{\Delta }\theta `$ is not yet clear (Panaitescu & Kumar 2002). If we assume the total emission energy to be constant as in the previous paper, the intrinsic luminosity is larger for the smaller opening half-angle. Such GRBs at the cosmological distance observed from off-axis viewing angle may be seen as XRFs and $`V/V_{\mathrm{max}}`$ is expected to be smaller than 0.5.

In this paper, we will show that our off-axis model has a possibility of accounting for the observational value of $`V/V_{\mathrm{max}}`$ if we change some of the model parameters in the previous paper (Yamazaki, Ioka, & Nakamura 2002). This paper is organized as follows. In § 2, we describe a simple jet model including the effect of cosmological expansion. We assume the uniform jet with sharp edges. Although one may consider the structured jet motivated by the simulation of the collapsar model, we cannot conclude, both observationally and theoretically, which model is preferable. The $`V/V_{\mathrm{max}}`$ for the XRFs detected by the Wide Field Cameras (WFCs) on BeppoSAX is calculated in § 3. Section 4 is devoted to a discussion. Throughout this paper, we adopt the flat universe with $`\mathrm{\Omega }_M=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ and $`h=0.65`$.

## 2 EMISSION MODEL OF X-RAY FLASHES

We consider a simple jet model of XRFs in the previous papers (Yamazaki, Ioka, & Nakamura 2002; Ioka & Nakamura 2001) taking into account the cosmological effect. A general formula to calculate the observed flux from an optically thin material is derived by Granot, Piran, & Sari (1999) and Woods & Loeb (1999). Here we adopt their formulations and notations. Let us use a spherical coordinate system $`𝒓=(r,\theta ,\varphi )`$ in the central engine frame, where the $`\theta =0`$ axis points toward the detector and the central engine is located at the origin. Consider a photon emitted at time $`t`$ and place $`𝒓`$ in the central engine frame. It will reach the detector at a time $`T`$ given by

$$T=(1+z)T_z=(1+z)(tr\mu /c),$$
(1)
where $`\mu \mathrm{cos}\theta `$ and $`z`$ is the cosmological redshift of the source; $`T=0`$ was chosen as the time of arrival at the detector of a photon emitted at the origin at $`t=0`$. Then the observed flux at the observed time $`T`$ and observed frequency $`\nu `$, measured in erg s-1 cm-2 Hz-1, is given by

$$F_\nu (T)=\frac{1+z}{d_L^2}_0^{2\pi }𝑑\varphi _1^1𝑑\mu _0^{\mathrm{}}r^2𝑑r\frac{j_\nu ^{}^{}(\mathrm{\Omega }_d^{},𝒓,T_z+r\mu /c)}{\gamma ^2(1\beta \mu )^2},$$
(2)
where $`d_L`$, $`\mathrm{\Omega }_d^{}`$ and $`j_\nu ^{}^{}`$ are the luminosity distance to the source, the direction towards the detector measured in the frame comoving with the jet (comoving frame), and the comoving frame emissivity in units of ergs s-1 cm-3 Hz-1 sr-1. The frequency $`\nu ^{}`$, which is measured in the comoving frame, is given by

$$\nu ^{}=\nu _z\gamma (1\beta \mu )=(1+z)\nu \gamma (1\beta \mu ).$$
(3)
We adopt an instantaneous emission of infinitesimally thin shell at $`t=t_0`$ and $`r=r_0`$. If the emission is isotropic in the comoving frame, the emissivity has a functional form of

$`j_\nu ^{}^{}(\mathrm{\Omega }_d^{},𝒓,t)`$ $`=`$ $`A_0f(\nu ^{})\delta (tt_0)\delta (rr_0)H(\mathrm{\Delta }\theta |\theta \theta _v|)`$ (4)
$`\times H\left[\mathrm{cos}\varphi \left({\displaystyle \frac{\mathrm{cos}\mathrm{\Delta }\theta \mathrm{cos}\theta _v\mathrm{cos}\theta }{\mathrm{sin}\theta _v\mathrm{sin}\theta }}\right)\right],`$
where $`f(\nu ^{})`$ and $`A_0`$ represent the spectral shape and the amplitude, respectively. $`\mathrm{\Delta }\theta `$ and $`\theta _v`$ are the jet opening half-angle and the viewing angle that the axis of the emission cone makes with the $`\theta =0`$ axis. The delta functions describe an instantaneous emission at $`t=t_0`$ and $`r=r_0`$, and $`H(x)`$ is the Heaviside step function which describes that the emission is inside the cone. Then the observed flux of a single pulse is given by

$`F_\nu (T)={\displaystyle \frac{2(1+z)r_0c\gamma ^2A_0}{d_L^2}}{\displaystyle \frac{\mathrm{\Delta }\varphi (T)f\left[\nu _z\gamma (1\beta \mathrm{cos}\theta (T))\right]}{\left[\gamma ^2(1\beta \mathrm{cos}\theta (T))\right]^2}},`$ (5)
where, $`1\beta \mathrm{cos}\theta (T)=(1+z)^1(c\beta /r_0)(TT_0)`$ and $`T_0=(1+z)(t_0r_0/c\beta )`$. For $`\mathrm{\Delta }\theta >\theta _v\mathrm{and}0<\theta (T)\mathrm{\Delta }\theta \theta _v`$, $`\mathrm{\Delta }\varphi (T)=\pi `$, otherwise $`\mathrm{\Delta }\varphi (T)=\mathrm{cos}^1\left[\frac{\mathrm{cos}\mathrm{\Delta }\theta \mathrm{cos}\theta (T)\mathrm{cos}\theta _v}{\mathrm{sin}\theta _v\mathrm{sin}\theta (T)}\right]`$. For $`\theta _v<\mathrm{\Delta }\theta `$, $`\theta (T)`$ varies from 0 to $`\theta _v+\mathrm{\Delta }\theta `$ while from $`\theta _v\mathrm{\Delta }\theta `$ to $`\theta _v+\mathrm{\Delta }\theta `$ for $`\theta _v>\mathrm{\Delta }\theta `$. In the latter case, $`\mathrm{\Delta }\varphi (T)=0`$ for $`\theta (T)=\theta _v\mathrm{\Delta }\theta `$. Pulse-starting and ending time are given by

$`T_{start}`$ $`=`$ $`T_0+(1+z)(r_0/c\beta )(1\beta \mathrm{cos}(\mathrm{max}[0,\theta _v\mathrm{\Delta }\theta ])),`$
$`T_{end}`$ $`=`$ $`T_0+(1+z)(r_0/c\beta )(1\beta \mathrm{cos}(\theta _v+\mathrm{\Delta }\theta )).`$ (6)
The observed spectrum of GRBs is well approximated by the Band spectrum (Band et al. 1993). In order to have a spectral shape similar to the Band spectrum, we adopt the following form of the spectrum in the comoving frame,

$$f(\nu ^{})=\{\begin{array}{cc}\hfill (\nu ^{}/\nu _0^{})^{1+\alpha _B},& \nu ^{}<\nu _0^{}\\ \hfill (\nu ^{}/\nu _0^{})^{1+\beta _B},& \nu ^{}>\nu _0^{}\end{array}$$
(7)
where $`\alpha _B`$ ($`\beta _B`$) is the low (high) energy power law index. Equations (5) and (7) are the basic equations to calculate the flux of a single pulse.

In order to study the dependence on the viewing angle $`\theta _v`$ and the jet opening half-angle $`\mathrm{\Delta }\theta `$, we fix the other parameters as $`\alpha _B=1`$, $`\beta _B=3`$, $`\gamma \nu _0^{}=200\mathrm{keV}`$, $`r_0/c\beta \gamma ^2=10\mathrm{s}`$, and $`\gamma =100`$ (Preece et al. 2000). We fix the amplitude $`A_0`$ so that the isotropic $`\gamma `$-ray energy $`E_{iso}=4\pi d_L^2(1+z)^1S(202000\mathrm{keV})`$ satisfies

$$\frac{(\mathrm{\Delta }\theta )^2}{2}E_{iso}=5\times 10^{50}\mathrm{ergs},$$
(8)
when $`\theta _v=0`$ and $`z=1`$ (Frail et al. 2001). Here $`S(\nu _1\nu _2)=_{T_{start}}^{T_{end}}F(T;\nu _1\nu _2)𝑑T`$ is the fluence in the energy range $`\nu _1\nu _2`$ and $`F(T;\nu _1\nu _2)=_{\nu _1}^{\nu _2}F_\nu (T)𝑑\nu `$ is the flux in the same energy range. The values of $`A_0`$ for different opening angles are summarized in Table 1. When the jet opening half-angle $`\mathrm{\Delta }\theta `$ becomes smaller, $`A_0`$ becomes larger.

## 3 CALCULATION OF $`V/V_{\mathrm{max}}`$

The formalism to calculate $`V/V_{\mathrm{max}}`$ is given by Mao & Paczyński (1992) and Piran (1992). In our case the absolute luminosity and spectrum depend on $`\theta _v`$ and $`\mathrm{\Delta }\theta `$. The observed-integrated number count from the source at redshift $`z`$ with the jet opening half-angle $`\mathrm{\Delta }\theta `$ and the viewing angle $`\theta _v`$ is given by

$$C(\mathrm{\Delta }\theta ,\theta _v,z)=_{T_{start}}^{T_{end}}𝑑T_{\nu _1}^{\nu _2}𝑑\nu \frac{F_\nu (T,\mathrm{\Delta }\theta ,\theta _v,z)}{h\nu }.$$
(9)
Let $`z_{min}`$ and $`z_{max}`$ be the minimum and the maximum redshift of the XRF for given $`\mathrm{\Delta }\theta `$ and $`\theta _v`$. In determining $`z_{min}`$ and $`z_{max}`$, we should note that the operational definition of the XRF detected by BeppoSAX is the fast X-ray transient with duration less than $`10^3`$ seconds which is detected by WFCs and not detected by the Gamma-Ray Burst Monitor (GRBM). Therefore, if the sources are nearby such that $`z<z_{min}`$, they are observed as GRBs because the observed fluence in the $`\gamma `$-ray band ($`40700`$ keV) becomes larger than the limiting sensitivity of GRBM ($`3\times 10^6\mathrm{ergs}\mathrm{cm}^2`$). If the sources are too far such that $`z>z_{max}`$, they cannot be observed by WFCs, which have an observation band of $`228`$ keV and a limiting sensitivity of about $`4\times 10^7\mathrm{ergs}\mathrm{cm}^2`$. Although the detection conditions of instruments vary with many factors of each event, such as the duration, the spectral index, or the peak photon energy (Band 2002), we adopt very simple criteria here. As shown in Figure 1, both $`z_{max}`$ and $`z_{min}`$ depend on $`\theta _v`$ and $`\mathrm{\Delta }\theta `$.

For given $`z_{max}`$, we can calculate the minimum integrated count $`C_{min}=C(\mathrm{\Delta }\theta ,\theta _v,z_{max})`$. Then, $`V/V_{\mathrm{max}}`$ for given $`\mathrm{\Delta }\theta `$ and $`\theta _v`$ is calculated as

$$\frac{V}{V_{max}}_{\mathrm{\Delta }\theta ,\theta _v}=\frac{_{z_{min}}^{z_{max}}[C(z)/C_{min}]^{3/2}𝒟(z)𝑑z}{_{z_{min}}^{z_{max}}𝒟(z)𝑑z},$$
(10)
where $`𝒟(z)`$ is given as

$$𝒟(z)=\frac{n(z)}{1+z}\mathrm{\hspace{0.17em}4}\pi \left(\frac{d_L}{1+z}\right)^2\frac{d}{dz}\left(\frac{d_L}{1+z}\right).$$
(11)
The luminosity distance is given by

$$d_L(z)=\frac{c}{H_0}(1+z)_0^z\frac{dz^{}}{\sqrt{\mathrm{\Omega }_\mathrm{\Lambda }+\mathrm{\Omega }_M(1+z^{})^3}}.$$
(12)
The function $`n(z)`$ is the comoving GRB rate density. We assume that $`n(z)`$ is proportional to the comoving rate densities of the star formation $`n_{SF}(z)`$ in the following form (Madau & Pozzetti 2000; Porciani & Madau 2001)

$$n_{SF}(z)=0.46h\frac{\mathrm{exp}(3.4z)}{\mathrm{exp}(3.8z)+45}\sqrt{\mathrm{\Omega }_M+\mathrm{\Omega }_\mathrm{\Lambda }(1+z)^3}M_{\mathrm{}}\mathrm{yr}^1\mathrm{Mpc}^3.$$
(13)
When we consider the case of $`z_{min}<z_{max}1`$, the cosmological effect can be neglected, so that $`𝒟z^2`$ and $`C/C_{min}(z/z_{max})^2`$. Then, Eq. (10) becomes

$$\frac{V}{V_{max}}_{\mathrm{\Delta }\theta ,\theta _v}=0.5\left[1+(z_{min}/z_{max})^3\right].$$
(14)
Therefore, $`z_{min}z_{max}1`$ implies $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta ,\theta _v}0.5`$.

Next, we integrate Eq. (10) over $`\theta _v`$ as

$$\frac{V}{V_{max}}_{\mathrm{\Delta }\theta }=\frac{V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta ,\theta _v}W(\theta _v)𝑑\theta _v}{W(\theta _v)𝑑\theta _v},$$
(15)
where $`W(\theta _v)`$ is the weight function which is the product of the solid angle factor and the volume factor:

$$W(\theta _v)=2\pi \mathrm{sin}\theta _v_{z_{min}}^{z_{max}}𝒟(z)𝑑z.$$
(16)
The results of the numerical integration are summarized in Table 1. For each $`\mathrm{\Delta }\theta `$, $`z_{max}(\theta _{v,p})`$ ($`z_{min}(\theta _{v,p})`$) means the maximum (minimum) redshift at which $`W(\theta _v)`$ takes the maximum value. If we take the jet opening half-angle as $`\mathrm{\Delta }\theta 0.03`$, $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta }`$ is smaller than $`0.4`$, which is consistent with the observational result at the 1 $`\sigma `$ level.

Let us consider the behavior of $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta ,\theta _v}`$ in the case of the fixed opening half-angle $`\mathrm{\Delta }\theta =0.03`$. Recent analysis of the GRB afterglows shows that some GRBs have a jet with a opening half-angle of less than 0.05 radian, i.e., the smallest value of $`\mathrm{\Delta }\theta `$ is about 0.03 (Panaitescu & Kumar, 2002). Let us decrease the viewing angle $`\theta _v`$ from a sufficiently large value ($`\gamma \theta _v20`$). As shown in Figure 1, both $`z_{max}`$ and $`z_{min}`$ increase monotonically, since the observed flux from the source increases due to the relativistic beaming effect. However, the behavior of $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta ,\theta _v}`$ is more complicated since the function $`𝒟(z)`$ in equation (11) has the maximum value at $`z=z_p1.5`$. We plot $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta ,\theta _v}`$ as a function of $`\gamma \theta _v`$ in Figure 2. If $`\theta _v`$ is large enough $`z_{max}`$ is smaller than $`z_p`$. Then the cosmological effect is small, so that $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta ,\theta _v}0.5`$. While if $`\theta _v`$ is small enough $`z_{min}`$ is larger than $`z_p`$. Then $`𝒟(z)`$ is a decreasing function in the range $`z_{min}<z<z_{max}`$, and the contribution of XRFs at smaller distance to $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta ,\theta _v}`$ is larger so that $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta ,\theta _v}`$ is small. We note that the above argument does not depend on $`\mathrm{\Delta }\theta `$ so much.

The behavior of the weight function $`W(\theta _v)`$ is shown in Figure 3. When $`\theta _v`$ is large enough for $`z_{max}`$ to be smaller than $`z_p`$ or when $`\theta _v`$ is small enough, $`W(\theta _v)`$ is small since $`𝒟(z)`$ or the solid angle factor are relatively small in the range $`z_{min}<z<z_{max}`$. One can see that in the case of $`\mathrm{\Delta }\theta =0.03`$, $`W(\theta _v)`$ takes the maximum value at $`\gamma \theta _{v,p}4.5`$, when $`z_{min}z_p`$.

Above discussions show that the sources with the viewing angle $`\theta _v0.05`$ at $`z1.5`$ are the most frequent class of the XRFs in the population for the opening half-angle $`\mathrm{\Delta }\theta =0.03`$. In other cases in which $`\mathrm{\Delta }\theta `$ is smaller than 0.03, the typical values are $`z1.5`$ and $`\theta _v\mathrm{\Delta }\theta +0.02`$ since $`W(\theta _v)`$ takes maximum at $`\theta _{v,p}\mathrm{\Delta }\theta +0.02`$ (see Table 1).

## 4 DISCUSSION

We have calculated $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta }`$ for the emission from a simple jet model, and shown that when the jet opening half-angle $`\mathrm{\Delta }\theta `$ is smaller than about 0.03, $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta }`$ for the XRFs detected by WFCs on BeppoSAX is smaller than 0.4. The value of $`\mathrm{\Delta }\theta 0.03`$ has been obtained from the fitting of the afterglow light curve (Panaitescu & Kumar 2002; Frail et al. 2001). Such a narrow jet which is inferred by afterglow observations is rare. However, following Eq. (1) of Frail et al. (2001), the jet break time is given by $`t_j13(\mathrm{\Delta }\theta /0.01)^{8/3}`$ minutes and so it requires fast localization to observe the jet break for a narrow jet. Therefore, at present, the small number of GRBs with small $`\mathrm{\Delta }\theta `$ may come from the observational selection effect. In the context of this scenario, we might be able to account for the fact that afterglows of XRFs have been rarely observed since the afterglow at a fixed time gets dimmer for an earlier break time. Furthermore, some “dark GRBs” might be such a small opening angle jet observed with an on-axis viewing angle for the same reason.

We briefly comment on how the results obtained in this paper will depend on the Lorentz factor of the jet $`\gamma `$. We see that when $`\gamma `$ becomes large, $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta }`$ becomes small. For example, when we fix $`\gamma =200`$, we obtain $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta =0.03}=0.39`$ and $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta =0.1}=0.43`$. This implies that the limitation on $`\mathrm{\Delta }\theta `$ can be relaxed.

We can estimate the typical observed photon energy as $`h\nu _{obs}(1+z)^1\delta \nu _0^{},`$ where $`\delta ^1=\gamma [1\beta \mathrm{cos}\stackrel{~}{\theta }]`$ and $`\stackrel{~}{\theta }=\mathrm{max}\{0,\theta _v\mathrm{\Delta }\theta \}`$ (Yamazaki, Ioka, & Nakamura 2002). Since $`\gamma 1`$ and $`\theta _v,\mathrm{\Delta }\theta 1`$, we obtain

$$h\nu _{obs}\frac{2\gamma \nu _0^{}}{(1+z)[1+(\gamma \stackrel{~}{\theta })^2]}.$$
(17)
In § 3, we have shown that for fixed $`\mathrm{\Delta }\theta 0.03`$, the typical value of $`\theta _v`$ is $`\mathrm{\Delta }\theta +0.02`$. Therefore, for the adopted parameters $`\gamma \nu _0^{}=200`$keV and the typical redshift $`z=z_p1.5`$, one can derive $`h\nu _{obs}30`$ keV, which is the typical observed peak energy of the XRFs (Kippen et al. 2002). We can propose from our argument that the emissions from the jets with a small opening half-angle such as $`\mathrm{\Delta }\theta 0.03`$ are observed as XRFs when they are seen from off-axis viewing angle. If one can detect the afterglow of the XRF, which has the maximum flux at about several hours after the XRF, the fitting of light curve may give us the key information about the jet opening angle (Granot et al. 2002). Therefore, our theoretical model can be tested by the near-future observations.

We can estimate the observed event rate of the XRF for fixed $`\mathrm{\Delta }\theta `$ as,

$$R_{\mathrm{XRF},\mathrm{\Delta }\theta }=\frac{1}{4\pi }W(\theta _v)𝑑\theta _v.$$
(18)
In order to calculate $`R_{\mathrm{XRF},\mathrm{\Delta }\theta }`$, we consider the proportionality constant $`=n(z)/n_{SF}(z)`$. One can write approximately as $`=r_{>8M_{\mathrm{}}}k`$, where $`r_{>8M_{\mathrm{}}}`$ is the number of stars with masses $`M>8M_{\mathrm{}}`$ per unit mass. Since we assume all stars with masses $`M>8M_{\mathrm{}}`$ explode as core-collapse SNe, $`k`$ represents the ratio of the number of XRF sources to that of core-collapse SNe. We adopt the value $`k=1\times 10^3`$, which is derived by the result of Porciani & Madau (2001) combined with the effect of the solid angle factor $`(\mathrm{\Delta }\theta )^2/2`$. Using a Salpeter initial mass function $`\varphi (M)`$, we obtain $`r_{>8M_{\mathrm{}}}=(_{8M_{\mathrm{}}}^{125M_{\mathrm{}}}\varphi (M)𝑑M)/(_0^{125M_{\mathrm{}}}M\varphi (M)𝑑M)=1.2\times 10^2M_{\mathrm{}}^1`$. Then, in the case of $`\mathrm{\Delta }\theta =0.03`$, we derive $`R_{\mathrm{XRF},\mathrm{\Delta }\theta }=1\times 10^2\mathrm{events}\mathrm{yr}^1(/1\times 10^5M_{\mathrm{}}^1)`$, which is comparable to the observed event rate of the XRF (Heise et al. 2001). Note that the value of $`R_{\mathrm{XRF},\mathrm{\Delta }\theta }`$ remains unchanged within a factor of 2 when we vary $`\mathrm{\Delta }\theta `$ from 0.01 to 0.07.

When the jet opening half-angle $`\mathrm{\Delta }\theta `$ has a distribution $`f_{\mathrm{\Delta }\theta }`$, we integrate $`V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta }`$ and $`R_{\mathrm{XRF},\mathrm{\Delta }\theta }`$ over the distribution of $`\mathrm{\Delta }\theta `$ as

$$V/V_{\mathrm{max}}=\frac{d(\mathrm{\Delta }\theta )f_{\mathrm{\Delta }\theta }R_{\mathrm{XRF},\mathrm{\Delta }\theta }V/V_{\mathrm{max}}_{\mathrm{\Delta }\theta }}{d(\mathrm{\Delta }\theta )f_{\mathrm{\Delta }\theta }R_{\mathrm{XRF},\mathrm{\Delta }\theta }},$$
(19)
$$R_{\mathrm{XRF}}=\frac{d(\mathrm{\Delta }\theta )f_{\mathrm{\Delta }\theta }R_{\mathrm{XRF},\mathrm{\Delta }\theta }}{d(\mathrm{\Delta }\theta )f_{\mathrm{\Delta }\theta }},$$
(20)
respectively. We assume a power-low distribution as $`f_{\mathrm{\Delta }\theta }(\mathrm{\Delta }\theta )^q`$. When we adopt $`q=4.54`$ (Frail et al. 2001), and integrate over $`\mathrm{\Delta }\theta `$ from 0.01 to 0.2 rad, we find $`V/V_{\mathrm{max}}=0.36`$ and $`R_{\mathrm{XRF}}=1\times 10^2`$ events yr-1. These values mainly depend on the lower bound of the integration. For example, we obtain $`V/V_{\mathrm{max}}=0.43`$ and $`R_{\mathrm{XRF}}=3`$ events yr-1 if the integration is done over $`\mathrm{\Delta }\theta `$ from 0.03 to 0.2 rad (Note that we may let the value of $`R_{\mathrm{XRF}}`$ be consistent with observed value by adjusting $``$). Since the statistics of the observational data will increase in the near future owing to instruments such as HETE-2 and Swift, we will be able to say more than above discussion, including more accurate functional form of $`f_{\mathrm{\Delta }\theta }`$ than that we have considered above, as well as the relation to the GRB event rate.

We would like to thank the referee for useful comments and suggestions. We are grateful to L. Amati for valuable comments. Numerical computation in this work was carried out at the Yukawa Institute Computer Facility. This work was supported in part by Grant-in-Aid for Scientific Research of the Japanese Ministry of Education, Culture, Sports, Science and Technology, No.05008 (R.Y.), No.00660 (K.I.), No.14047212 (T.N.), and No.14204024 (T.N.).


format: text

text:
# General Stability Analysis of Synchronized Dynamics in Coupled Systems

## Figure Caption

Figure 1: Schematic illustrations of the stability zone. (a) unbounded area, (b) bounded area.


text:
# Matrix controlled channel diffusion of sodium in amorphous silica

## Abstract

To find the origin of the diffusion channels observed in sodium-silicate glasses, we have performed classical molecular dynamics simulations of Na2O–4SiO2 during which the mass of the Si and O atoms has been multiplied by a tuning coefficient. We observe that the channels disappear and that the diffusive motion of the sodium atoms vanishes if this coefficient is larger than a threshold value. Above this threshold the vibrational states of the matrix are not compatible with those of the sodium ions. We interpret hence the decrease of the diffusion by the absence of resonance conditions.

The mechanism of ionic transport in amorphous materials  at the atomic scale is not completely elucidated so far and therefore it is the subject of numerous experimental  and numerical  studies. Among the most studied systems are the sodium-silicate glasses since they contain the essential ingredients, namely the amorphous matrix (silica) and the mobile ions (sodium), as a first step in the simulation of more complex glasses of higher practical interest.

In previous studies , we have shown by means of classical molecular dynamics simulations of Na2O–4SiO2 (NS4), that the sodium atoms diffuse through a well connected network of pockets (which represents only a limited fraction of the entire available space) that we have called “channels” to be coherent with the literature. The existence of the channels, which are not due to micro-segregation effects , has been confirmed by the existence of a pre-peak in the partial Na-Na structure factor at a wave-vector $`q=0.95`$ Å-1. This pre-peak has also been observed experimentally  and numerically  in another study. We have also shown that the location of the channels is strongly correlated to the positions of the non-bridging oxygens  and Horbach et al. have shown that the sodium dynamics should be related to that of the underlying silica network . This suggests that the origin of the channels could be related to the dynamical properties of the matrix. To check this idea we present in this letter classical molecular dynamics simulations on a series of “toy” systems in which the atomic masses of both the oxygen and silicon atoms have been systematically changed after artificially multiplying their experimental values by a common factor $`\mu `$ varying from 0.5 to infinity, the usual NS4 system being recovered for $`\mu =1`$.

By studying both the mean square displacement of the sodium atoms and the characteristics of the channels, we find a change in the sodium diffusion properties when the parameter $`\mu `$ is increased above a value of about 30. Above this threshold, the sodium diffusion decreases and the channels can no more be clearly defined. Guided by the concomitant change in the short time characteristics of the velocity autocorrelation function, we have calculated the vibrational density of states (VDOS) for both the sodium atoms and the atoms of the matrix for various values of $`\mu `$. We observe that the threshold corresponds to the value of $`\mu `$ above which the sharp VDOS of the sodium atoms starts to separate from the band of vibrational states of the matrix. Therefore we propose that the diffusive motion of the sodium atoms, i.e. their ability to escape from their local cage, is facilitated when their vibrational frequencies are compatible with the ones of the matrix. We argue that this mechanism is responsible for the channel diffusion of sodium in the silica matrix.

In this letter we present classical molecular dynamics calculations of a system of 648 particles (86 sodium, 173 silicon and 389 oxygen atoms) confined in a cubic box of edge length 20.88 Å with periodic boundary conditions. The density is thus the experimental density of glassy NS4, i.e. 2.38 g.cm-3 . The interactions between the particles are given by a modified version of the so-called “BKS” potential  which is able to reproduce the structure as well as the dynamics of several sodium-silicate systems  (for more details see ). In the present study we have generated three independent samples (in the following, all the results are averaged over these samples in order to improve the statistics) at a temperature of $`1900`$ K for which the channels have been evidenced in our previous studies. For each sample we have performed ten different simulations in which the mass of the atoms of the matrix has been artificially multiplied by a factor $`\mu `$ of 0.5, 1, 2, 5, 10, 30, 102, 103, 104, and 106. We have also carried out the limiting case $`\mu =\mathrm{}`$ by performing simulations in which only the sodium atoms move, the atoms of the matrix being kept fixed (frozen matrix approximation). Each of these 33 samples has first been relaxed at $`1900`$ K for 106 steps (1 step = 1.4 fs) and the following 106 steps were used to produce 1000 configurations recorded every 1000 time steps. These configurations have then been used to analyze the trajectories of the sodium atoms during these 1.4 ns.

In Fig. 1 we show $`R^2(t)`$, the mean square displacement (MSD) of the sodium atoms, for various values of $`\mu `$. While these curves have the characteristic long-time shape of strong diffusion for small values of $`\mu `$, they flatten out as $`\mu `$ increases. In the limiting case $`\mu =\mathrm{}`$ (frozen matrix) the curve becomes even completely flat, at least within the period of time covered by our simulations. In order to give a quantitative idea of how the ionic diffusion properties decrease with increasing $`\mu `$, we have determined a characteristic time $`\tau _{\mathrm{MSD}}`$ necessary for $`R^2(t)`$ to be equal to 4 Å2. This corresponds to an average travel distance of 2 Å which is the minimum distance between two Na neighbors (as determined from the Na–Na radial pair distribution function ). In the inset of Fig. 1 we give the variation of $`\tau _{\mathrm{MSD}}`$ with $`\mu `$. While $`\tau _{\mathrm{MSD}}`$ is almost independent of $`\mu `$ for small $`\mu `$, it starts to increase for values of $`\mu `$ larger than a value of about 30.

Since we have previously shown that the diffusive motion of the sodium atoms takes place within channels , it is worth studying how these channels are modified when the mass of the atoms of the matrix is changed. Therefore, for all our values of $`\mu `$, we have determined the channels, in the same manner as previously : using a 3-dimensional mesh we determine the number of different sodium atoms that have visited each cube of the mesh during the total simulation time. Then we define $`\xi `$ which is the minimal occupation number such that the cubes visited more than $`\xi `$ times represent the upper 10 % of all the visited cubes (for more details see ). In previous studies we have shown that at $`T1900`$ K a cube needed to be visited by at least $`\xi =8`$ different sodium atoms during a 1.4 ns simulation, in order to be part of the channel structure.

Here we study how $`\xi `$ changes with $`\mu `$ and the results (averaged over three samples) are reported in Fig. 2. While for $`\mu <30`$, $`\xi `$ remains almost constant between 7 and 8 within the error bars, it starts to decrease dramatically for $`\mu >30`$. Obviously, as explained in refs. , one can no more speak of “channels” if $`\xi `$ becomes as small as unity. Therefore one can interpret the results in Fig. 2 by assuming that the channels disappear progressively as soon as $`\mu >30`$. Then one can argue that the increase of the characteristic time for the diffusion observed in Fig. 1 is intimately correlated with the disappearance of the channels shown in Fig. 2.

In order to know if the observed change in the long time diffusive behavior of the sodium atoms is accompanied by a change in their short time dynamical properties, we have calculated the Na velocity autocorrelation function $`\vartheta (t)=\stackrel{}{v}(t_0)\stackrel{}{v}(t_0+t)/v^2(t_0)`$ and studied its short time behavior. As seen in Fig. 3 the short time part of $`\vartheta (t)`$ is characterized by a first minimum at $`\tau _\vartheta `$, typical of the time for a sodium atom to bounce back and forth against the internal boundaries of the cage in which it vibrates (well before being eventually able to jump towards another cage). In the inset of Fig. 3 the intensity of the minimum $`\vartheta _\mathrm{m}=\vartheta (\tau _\vartheta )`$ has been plotted as a function of $`\mu `$. $`\vartheta _\mathrm{m}0.28`$ for $`\mu <10`$ while $`\vartheta _\mathrm{m}`$ drops down to $`0.5`$ for $`\mu >30`$. This increasing depth can be interpreted as an increasing stiffness of the internal cage boundaries when $`\mu `$ is increased. In a classical mechanics picture, due to this increased stiffness, the probability for a sodium atom to escape from its cage becomes weaker and this is consistent with the vanishing diffusion observed in Fig. 1.

To go further in our microscopical analysis, and to understand the changes occurring around $`\mu 30`$, we have calculated the VDOS of the sodium atoms and the atoms of the matrix, by Fourier transforming the corresponding velocity autocorrelation functions. It is known that such a method reproduces only approximately the VDOS, since at finite temperature (here $`T1900`$ K), the harmonic hypothesis does not fully hold but this approximation will not be crucial for our arguments. In Fig. 4 (a) the VDOS of the oxygen atoms for different values of $`\mu `$ are represented (a similar picture could be drawn for the Si atoms). The VDOS for $`\mu =1`$ is of course close to the one of amorphous silica obtained experimentally except that it is less structured, which is a known drawback of the BKS potential . Such a VDOS is typical of a broad band of vibrational states, extending from $`\nu =0`$ up to $`\nu _1`$ 35 THz, as a result of the coupling between neighboring oxygen and silicon atoms forming a strong random covalent network. When increasing the mass of the atoms of the matrix, one observes a shrinkage of the VDOS towards low frequencies. The top of the oxygen band, $`\nu _{\mathrm{max}}^\mathrm{O}`$, has been plotted versus $`\mu `$ in Fig 4 (c): as expected from standard solid state physics $`\nu _{\mathrm{max}}^\mathrm{O}`$ varies like $`\nu _1\mu ^{1/2}`$. The situation is quite different for the Na VDOS represented in Fig. 4 (b) for three typical values of $`\mu `$: it looks like a well defined peak centered at a frequency $`\nu _1^{\mathrm{Na}}`$ close to 5 THz for small values of $`\mu `$. Such a peak can be interpreted as being mostly due to the vibrations of the sodium atoms in their cage. Using the approximate formula for such a motion, $`R^1\sqrt{kT/m_{\mathrm{Na}}}`$, where $`R=1`$ Å is a typical size of the cage deduced from the MSD and $`m_{\mathrm{Na}}`$ is the mass of a sodium atom, one finds a frequency of about 8 THz in reasonable agreement with $`\nu _1^{\mathrm{Na}}`$. The broadening of the peak is then certainly due to the polydispersity of the $`R`$ values as well as to the coupling of the sodium atoms with other species. With increasing $`\mu `$ a second peak at $`\nu _2^{\mathrm{Na}}`$ grows out of the high frequency part and increases while the peak at $`\nu _1^{\mathrm{Na}}`$ decreases and shifts towards low frequencies. The variation of $`\nu _1^{\mathrm{Na}}`$ and $`\nu _2^{\mathrm{Na}}`$ with $`\mu `$ is represented in Fig. 4 (c) where the main peak of the Na VDOS is represented by the filled symbols. One sees that up to $`\mu =30`$ the peak at $`\nu _1^{\mathrm{Na}}`$ is the main peak while for higher values of $`\mu `$ the peak at $`\nu _2^{\mathrm{Na}}`$ becomes the main peak. Once the transition is fulfilled, the secondary Na peak follows the variation of $`\nu _{\mathrm{max}}^\mathrm{O}`$ with $`\mu `$ while the principal peak remains at a constant frequency. Such a behavior is typical of hybridization effects as commonly seen in electronic and vibrational systems. It is due to the coupling between the sodium atoms and the atoms of the matrix. The characteristic value $`\mu _0`$ for which the principal peak in the Na VDOS starts to escape from the broad oxygen band can be estimated by equating the top of the band $`\nu _1\mu _0^{1/2}`$ and the typical sodium frequency, 5 THz. This gives $`\mu _050`$ (indicated by the arrow in Fig. 4 (c)), which is remarkably close to the value of 30 above which we have observed the change in the diffusive properties of the Na atoms.

Therefore the above analysis of the vibrational density of states provides a very simple explanation of the predominance of channel diffusion for $`\mu `$ smaller than $`\mu _0`$. In that case the sodium frequency peak lies within the limits of the vibrational band of the matrix as shown in Fig. 4 (c). Therefore there always exists a vibrational mode of the matrix with the same frequency as the one of the sodium atoms. Using a simple picture in direct space, the vibrational amplitudes of the sodium atoms can become very large due to “resonance conditions”, giving them the possibility to escape from their cage and to jump to a neighboring cage, as it is expected in a local picture of the diffusive motion. Moreover, at the location where a sodium atom is connected to an oxygen atom via a covalent-like bond, it is clear that the local hybridization is important and this explains why the location of the channels, i.e. the preferential pathways for the sodium diffusive motion, is strongly correlated to the position of the non-bridging oxygen atoms. On the contrary, for $`\mu >\mu _0`$, the vibrational peak of the sodium atoms is located outside the VDOS of the matrix (Fig. 4 (c)) and the resonance conditions are more difficult to be fulfilled. As a consequence, the sodium atoms stay much longer confined in their cages. This picture is consistent with our interpretation of the increasing depth of the first minimum of the Na velocity autocorrelation function, as being due to an increasing stiffness of the internal cage boundaries when $`\mu `$ is increased above $`\mu _0`$.

In conclusion, we have demonstrated the essential role played by the vibrations of the atoms of the matrix in the existence of the channel diffusion of sodium inside an amorphous silica matrix. In a sense this is close to the concept of “matrix-mediated-coupling” recently used to interpret the mixed cation effect in glasses  except that we show here that the direct hybridization between the ionic modes and the modes of the matrix is necessary to insure a fast diffusion process of the ions. This result is also coherent with previous studies showing indirectly that the sodium dynamics should be intimately linked to the one of the matrix . Of course it would be extremely interesting to extend the present study to other kinds of ions, such as Li, H, etc., to test the generality of our interpretation (according to preliminary results it appears indeed that a similar behavior is observed for Li in SiO2 ). This would constitute a real improvement in the understanding of the mechanisms of ionic transport in random media.

We thank Pr. K. Funke and Pr. A. Heuer for interesting discussions. Part of the numerical calculations were done at the “Centre Informatique National de l’Enseignement Supérieur” (CINES) in Montpellier.


text:
# Flavor Symmetry Breaking and Strangeness in the Nucleon

## 1 Introduction and Motivation

There has been much interest in the strangeness content of the nucleon ever since the analysis of the DIS data EMC  suggested a large (negative) polarization of strange quarks in the nucleon Br88 ; El95 , $`\mathrm{}S_N0.15`$. This surprising result particularly relies on the assumption of flavor covariance for the axial current matrix elements of the $`\frac{1}{2}^+`$ baryons. This assumption originates from the feature that the Cabibbo scheme Ca63 , that utilizes the $`F`$&$`D`$ parameterization for the flavor changing axial charges, works unexpectedly well Fl98  as the comparison in table 1 shows.

Here we will investigate in how far this agreement justifies to carry over flavor covariance to strangeness conserving axial current matrix elements in order to disentangle the various quark flavor components of the nucleon axial current matrix element. This investigation requires baryon axial current matrix elements as functions of the (effective) strength of flavor symmetry breaking. This can be achieved within the three flavor version of the Skyrme model (and generalizations thereof) in which baryons emerge as solitons. In such models baryon states are constructed by quantizing the large amplitude fluctuations about the soliton and constructing exact eigenstates in the presence of symmetry breaking. We focus on a picture where symmetry breaking mainly resides in the baryon wave–functions, including important contributions which would be missed in a first order treatment. In contrast, we assume that the current operators, from which the charges are computed, are dominated by flavor covariant components. In a first step we do not specify the model Lagrangian but adjust the prefactors of the few possible flavor covariant components of the axial current operator to observables in hyperon beta–decay and analyze their matrix elements as functions of flavor symmetry breaking. We also present results obtained from a realistic vector meson soliton model that supports the suggested picture. Details omitted here may be traced from ref We00 .

## 2 Symmetry Breaking in the Baryon Wave–Functions

The collective coordinates $`A`$ that parameterize the large amplitude fluctuations off the soliton are introduced via

$$U(\stackrel{}{r},t)=A(t)U_0(\stackrel{}{r})A^{}(t),A(t)SU(3).$$
(1)
$`U_0(\stackrel{}{r})`$ describes the soliton embedded in the isospin subgroup. A prototype model Lagrangian for $`U(\stackrel{}{r},t)`$ consists of the Skyrme model supplemented by the Wess–Zumino–Witten term and suitable symmetry breaking pieces. We parameterize the collective coordinates by “Euler–angles”

$$A=D_2(\widehat{I})\mathrm{e}^{i\nu \lambda _4}D_2(\widehat{R})\mathrm{e}^{i(\rho /\sqrt{3})\lambda _8}.$$
(2)
Here $`D_2`$ denote the rotation matrices for rotations in isospace ($`\widehat{I}`$) and coordinate–space ($`\widehat{R}`$). Substituting (1) into the model Lagrangian yields upon canonical quantization the Hamiltonian for the collective coordinates $`A`$:

$$H=H_\mathrm{s}+\frac{3}{4}\gamma \mathrm{sin}^2\nu .$$
(3)
The symmetric piece of this collective Hamiltonian only contains Casimir operators and may be expressed in terms of the $`SU(3)`$–right generators $`R_a(a=1,\mathrm{},8)`$:

$$H_\mathrm{s}=M_{\mathrm{cl}}+\frac{1}{2\alpha ^2}\underset{i=1}{\overset{3}{}}R_i^2+\frac{1}{2\beta ^2}\underset{\alpha =4}{\overset{7}{}}R_\alpha ^2.$$
(4)
$`M_{\mathrm{cl}},\alpha ^2,\beta ^2`$ and $`\gamma `$ are functionals of the soliton, $`U_0(\stackrel{}{r})`$. The generators $`R_a`$ can be expressed in terms of derivatives with respect to the ‘Euler–angles’. The essential feature of the parameterization (2) is that the flavor symmetry breaking part of the full Hamiltonian (3) only depends on the flavor changing angle $`\nu `$. Therefore the eigenvalue problem $`H\mathrm{\Psi }=ϵ\mathrm{\Psi }`$ reduces to ordinary second order differential equations for isoscalar functions which only depend on $`\nu `$ Ya88 . Solely the product $`\omega ^2=\frac{3}{2}\gamma \beta ^2`$ appears in these differential equations as the effective strength of symmetry breaking on which the eigenfunctions of $`H`$ depend parametrically. A value in the range $`5\stackrel{<}{}\omega ^2\stackrel{<}{}8`$ is required to obtain reasonable agreement with the empirical mass differences for the $`\frac{1}{2}^+`$ and $`\frac{3}{2}^+`$ baryons We96 . Such large a value for $`\omega ^2`$ is without reach of a perturbation expansion as the resulting baryon wave–functions exhibit strong distortion from flavor covariance.

## 3 Charge Operators

In the soliton description the effect of the derivative type symmetry breaking terms is mainly indirect. They provide the splitting between the various decay constants and thus increase $`\gamma `$ since it is proportional to $`f_K^2m_K^2f_\pi ^2m_\pi ^21.5f_\pi ^2(m_K^2m_\pi ^2)`$. Otherwise the derivative type symmetry breaking terms are negligible. Whence symmetry breaking terms can be omitted in the current operators and the non–singlet axial charge operator is parameterized as ($`a=1,\mathrm{},8,i=1,2,3`$)

$$d^3rA_i^{(a)}=c_1D_{ai}c_2D_{a8}R_i+c_3\underset{\alpha ,\beta =4}{\overset{7}{}}d_{i\alpha \beta }D_{a\alpha }R_\beta ,$$
(5)
where $`D_{ab}=\frac{1}{2}\mathrm{tr}\left(\lambda _aA\lambda _bA^{}\right)`$. For $`\omega ^2\mathrm{}`$ (infinitely heavy strange degrees of freedom) the strangeness contribution to the nucleon axial charge should vanish. Noting that $`N|D_{83}|N0`$ and $`N|_{\alpha ,\beta =4}^7d_{3\alpha \beta }D_{8\alpha }R_\beta |N0`$ while $`N|D_{88}|N1`$ for $`\omega ^2\mathrm{}`$, we demand

$$d^3rA_i^{(0)}=2\sqrt{3}c_2R_ii=1,2,3$$
(6)
for the axial singlet current because it leads to the strangeness projection, $`A_i^{(s)}=(A_i^{(0)}2\sqrt{3}A_i^{(8)})/3`$ that vanishes for $`\omega ^2\mathrm{}`$. Actually all model calculations in the literature Pa92 ; Bl93  are consistent with this relation between singlet and octet currents. The singlet current matrix element, $`\mathrm{}\mathrm{\Sigma }_B=\sqrt{3}c_2`$, is the quark spin contribution to the spin of the considered baryon, $`B`$. It is well known that the empirical value for the nucleon matrix element, $`\mathrm{}\mathrm{\Sigma }_N0.20\pm 0.10`$ El95  is insensitive to the strength of flavor symmetry breaking Jo90 . This suggests to adjust $`c_2`$ accordingly. In order to completely describe the hyperon beta–decays we also demand matrix elements of the vector charges. These are obtained from the operator

$$d^3rV_0^{(a)}=\underset{b=1}{\overset{8}{}}D_{ab}R_b=L_a,$$
(7)
which introduces the $`SU(3)`$–left generators $`L_a`$.

The values for $`g_A`$ and $`g_V`$ (only $`g_A`$ for $`\mathrm{\Sigma }^+\mathrm{\Lambda }e^+\nu _e`$) are obtained from the matrix elements of the operators in eqs (5) and (7), respectively, sandwiched between the eigenstates of the full Hamiltonian (3). We still have to specify $`c_1`$ and $`c_3`$. We determine these two parameters such that nucleon axial charge, $`g_A`$ and the $`g_A/g_V`$ ratio for $`\mathrm{\Lambda }pe^{}\overline{\nu }_e`$ are reproduced111In this section we will not address the problem of the too small model prediction for $`g_A`$. at a prescribed strength of flavor symmetry breaking, $`\omega _{\mathrm{fix}}^2=6.0`$. Then we are not only left with predictions for the other decay parameters but we can in particular study the variation with symmetry breaking. This is shown in figure 1.

The dependence on flavor symmetry breaking is very moderate222However, the individual matrix elements entering the ratios $`g_A/g_V`$ vary strongly with $`\omega ^2`$ We00 . and the results can be viewed as reasonably agreeing with the empirical data, cf. table 1. The observed independence of $`\omega ^2`$ shows that these predictions are not sensitive to the choice of $`\omega _{\mathrm{fix}}^2`$. The two transitions, $`np`$ and $`\mathrm{\Lambda }p`$, which are not shown in figure 1, exhibit a similar negligible dependence on $`\omega ^2`$. We therefore have a two parameter ($`c_1`$ and $`c_3`$, $`c_2`$ is fixed from $`\mathrm{\Delta }\mathrm{\Sigma }_N`$) fit of the hyperon beta–decays. Comparing the results in figure 1 with the data in table 1 we see that the present calculation using the strongly distorted wave–functions agrees equally well with the empirical data as the flavor symmetric $`F`$&$`D`$ fit. On the other hand, the strangeness contribution to the nucleon axial current matrix element reduces from $`\mathrm{}S_N0.13`$ in the symmetric treatment to $`\mathrm{}S_N0.07`$ in the realistic case.

## 4 Model Calculation

We consider a realistic soliton model containing pseudoscalar and vector meson fields. It has been established for two flavors in ref Ja88  and been extended to three flavors in ref Pa92  where it has been shown to fairly well describe the parameters of hyperon beta–decay (cf. table 4 in ref Pa92 ). The model Lagrangian contains terms which involve the Levi–Cevita tensor $`ϵ_{\mu \nu \rho \sigma }`$, to accommodate processes like $`\omega 3\pi `$ Ka84 . Such terms contribute to $`c_2`$ and $`c_3`$. A minimal set of symmetry breaking terms is included Ja89  to account for different masses and decay constants. They add symmetry breaking pieces to the axial charge operator,

$`\delta A_i^{(a)}`$ $`=`$ $`c_4D_{a8}D_{8i}+c_5{\displaystyle \underset{\alpha ,\beta =4}{\overset{7}{}}}d_{i\alpha \beta }D_{a\alpha }D_{8\beta }`$
$`+c_6D_{ai}(D_{88}1),`$
$`\delta A_i^{(0)}`$ $`=`$ $`2\sqrt{3}c_4D_{8i}.`$
Unfortunately the model parameters cannot be completely determined in the meson sector Ja88 . We use the remaining freedom to accommodate baryon properties in three different ways as shown in table 2. The set denoted by ‘b.f.’ refers to a best fit to the baryon spectrum. It predicts the axial charge somewhat on the low side, $`g_A=0.88`$. The entry ‘mag.mom.’ labels parameters that yield magnetic moments close to the empirical data (with $`g_A=0.98`$) and finally the set labeled ‘$`g_A`$’ reproduces the axial charge of the nucleon Pa92 .

We observe that in particular the strangeness projection of the nucleon axial current is very small and depends only mildly on the model parameters. This confirms the above conclusion from the general structure of the axial current matrix elements that the strangeness admixture in the nucleon is significantly smaller than an analysis based on flavor covariance suggests. Also the predictions for the axial properties of the $`\mathrm{\Lambda }`$ hyperon are quite insensitive to the model parameters. Sizable polarizations of the up and down quarks in the $`\mathrm{\Lambda }`$ are predicted; comparable to those obtained from the SU(3) analysis Ja96  of the available data.

## 5 Conclusions

We have suggested a picture for the axial charges of the low–lying $`\frac{1}{2}^+`$ baryons which manages to reasonably reproduce the empirical data without introducing (significant) flavor symmetry breaking components in the corresponding operators. Rather, a sizable symmetry breaking, as demanded by the baryon spectrum, resides almost completely in the baryon wave–functions. In this picture the empirical data for hyperon beta–decay are as reasonably reproduced as in the Cabibbo scheme. We emphasize that the present picture is not a re–application of the Cabibbo scheme since here the ‘octet’ baryon wave–functions have significant admixture of higher dimensional representations. Especially, when compared with the flavor covariant treatment, the present approach predicts a sizable suppression of strangeness in the nucleon.

This work is supported by the Deutsche Forschungsgemeinschaft (DFG) under contract We 1254/3-2.


text:
# Modularity and Extreme Edges of the Internet

## Abstract

We study the spectral properties of a diffusion process taking place on the Internet network focusing on the slowest decaying modes. These modes allow us to identify an underlying modular structure of the Internet roughly corresponding to individual countries. For instance in the slowest decaying mode the diffusion current flows from Russia towards US military sites. These two regions thus constitute the extreme edges of the Internet. Quantitatively the modular structure of the Internet manifests itself in approximately 10 times larger participation ratio of its slow decaying modes compared to the null model – a random scale-free network. We propose to use the fraction of nodes participating in slow decaying modes as a general measure of the modularity of a network. For the 100 slowest decaying modes of the Internet we measured this fraction to be around 30%. Finally we suggest, that the degree of isolation of an individual module can be assessed by comparing its participation in different diffusion modes. Using the proportionality of response as a criterion we find that the independent module approximation works well for the Internet.

Virtually any complex system has an underlying network that defines the backbone of interactions among its components. Examples of such networks include the Internet and the World Wide Web, molecular networks of living cells, food webs in ecosystems, etc. An important question is whether nodes of such a network can be divided into smaller sub-networks (modules), which interact with each other relatively weakly Girvan2002 . Estimating the strength of inter-modular interactions, localizing crucial links connecting these modules to each other, and finding pairs of modules which are the most distant from each other is important for several reasons. First of all, it serves as a test of stability of the system with respect to breaking it up into truly isolated components. Such a break-up would be undesirable in, for example, the Internet. Creation of extra connections between the most distant modules in the network and reinforcement of crucial links is an efficient way to increase its stability. Secondly, by measuring the relative strength of inter- and intra-modular connections one directly assesses the quality of the independent module approximation, which may turn out to be important in modeling the actual dynamics of a given complex system.

In this work we explore the modular structure present in the physical layout of the Internet. To this end we study an auxiliary diffusion process taking place on this network. The slowest modes of diffusion, easily identifiable from the spectrum of its transfer matrix, allow us to detect the weakly interacting modules of the Internet. These modules turn out to roughly correspond to individual countries or for large countries to cultural or geographical regions within the country. Of course, the diffusion process studied in this work does not describe the real dynamics of the information flow over the Internet. However, the detected modular features play an important role in any local dynamical process taking place on this network including the real Internet traffic.

Analysis of spectral properties of a similar diffusion process lies at the heart of the popular search engine www.google.com google . Its variants have also been applied to social networks (the correspondence analysis) wasserman\_book , random and small-world networks (the Laplace equation analysis) monasson , artificial scale-free networks, FarkasBarabasiSpectra ; GohKahngKimSpectra  and the community structure of the World Wide Web Kleinberg .

In this work we explore the physical layout of the Internet on a coarse-grained level of the so-called Autonomous Systems (AS), which are large groups of routers and servers belonging to one organization such as a university or a business enterprise (e.g. an Internet Service Provider). To this end we use the January 3, 2000 dataset when the Internet consisted of 6474 Autonomous Systems exchanging information via 12572 undirected links AS . As expected for the Internet, any pair of Autonomous Systems is connected to each other by at least one path, so that topologically the network consists of just one large cluster. The diffusion process we analyze here describes the dynamics of a large number of random walkers moving on the network at discrete time steps. Statistical properties of returns to the origin of such random walks have recently been used to measure the effective dimensionality of several complex networks Bilke2002 . At each time-step every walker moves from its current node to one of the neighboring nodes along a randomly selected link. The average dynamics of this process is described by

$$\rho _i(t+1)=\underset{j}{}T_{ij}\rho _j(t),$$
(1)
where $`\rho _i(t)`$ is the expectation value of the number of random walkers at site $`i`$ and time $`t`$. The elements $`T_{ij}`$ of the transfer matrix are equal to $`1/K_j`$ for neighboring nodes $`i`$ and $`j`$ and zero otherwise. Here $`K_j`$ is the connectivity (the number of immediate neighbors) of the node $`j`$ from which a walker steps to the node $`i`$. Note that $`_iT_{ij}=1`$, so that the total number of walkers is conserved at all times. Eq. (1) can also be rewritten as a discrete time diffusion equation $`\rho _i(t+1)\rho _i(t)=_j(T_{ij}\delta _{ij})\rho _j(t)`$. Hence the diffusion matrix $`𝐃`$ is related to the transfer matrix $`𝐓`$ simply as

$$𝐃=𝐓\mathrm{𝟏}.$$
(2)
As time advances the distribution of random walkers approaches a steady state $`\rho _i(\mathrm{})`$ in which the diffusion current flowing from a node $`i`$ to a node $`j`$ is exactly balanced by that flowing from $`j`$ to $`i`$. This is satisfied when the average number of walkers $`\rho _i(\mathrm{})`$ on every node $`i`$ is proportional to its connectivity $`K_i`$.

The relaxation of any initial distribution of random walkers among nodes, $`\rho _i(0)`$, towards the steady state configuration $`\rho _i(\mathrm{})`$ is determined by the spectral properties of the matrix $`𝐓`$ CommentSymT  (or alternatively $`𝐃`$). For instance, the steady state configuration $`\rho _i(\mathrm{})`$ itself is proportional to the principal eigenvector $`\rho _i^{(1)}`$ of $`𝐓`$ corresponding to its largest eigenvalue $`\lambda ^{(1)}=1`$, which is unique for single component networks such as the Internet. The remaining eigenvectors $`\rho ^{(\alpha )}`$ describe the decay of the initial configuration towards the steady state with a characteristic decay time $`\tau ^{(\alpha )}`$ related to the corresponding eigenvalue $`\lambda ^{(\alpha )}`$ through $`\mathrm{exp}(1/\tau ^{(\alpha )})=|\lambda ^{(\alpha )}|`$. Note that in general there exist both non-oscillatory ($`\lambda ^{(\alpha )}1`$) and oscillatory ($`\lambda ^{(\alpha )}1`$) slowly decaying modes.

The modularity of a given complex network reflects itself in statistical properties of its diffusion eigenvectors $`\rho _i^{(\alpha )}`$. One such property is the Participation Ratio (PR), which quantifies the effective number of nodes participating in a given eigenvector with a significant weight. In the Internet the components of the principal eigenvector $`\rho _i^{(1)}\rho _i(\mathrm{})K_i`$ as well as those of other slow decaying modes are broadly distributed (scale-free) Faloustos  and as such tend to be localized on just a few highly connected nodes. In this case participation ratios are best calculated using the normalized eigenvector

$$c_i^{(\alpha )}=\rho _i^{(\alpha )}/K_i$$
(3)
of outgoing currents flowing from node $`i`$ along each of its links. More formally $`c_i^{(\alpha )}`$ is also the eigenvector of the transposed transfer matrix $`𝐓^{}`$ with the same eigenvalue $`\lambda ^{(\alpha )}`$. For such a vector normalized by $`c_i^2=1`$ the participation ratio is defined as: $`PR=\left(_{i=1}^Nc_i^4\right)^1`$.

In Fig. 1 the participation ratio of eigenvectors $`c_i^{(\alpha )}`$ (top) and the eigenvalue density (i.e. the spectrum of the matrix) (bottom) is plotted as a function of the corresponding eigenvalue $`1<\lambda _i^{(\alpha )}<1`$.

The data for the Internet (filled circles) is displayed together with the data for its randomized counterpart (open squares). The randomization of the Internet was performed in such a way that the connectivity of every node is strictly preserved maslov\_sneppen\_internet . It was argued maslov\_sneppen\_science ; maslov\_sneppen\_internet  that such a network constitutes a proper null model of a given complex network. Since this random network has the same scale-free distribution of connectivities as the Internet Faloustos  it will be referred to as a Random Scale-Free Network (RSFN).

Comparing the data for the Internet and an RSFN we note that while the density of states is rather similar in these two networks (Fig. 1 B ), the participation ratio of the slowly decaying modes (especially for the non-oscillatory ones with $`\lambda `$ close to 1) is markedly higher in the Internet than in an RSFN (Fig. 1 A). In these non-oscillatory modes the diffusion current flows from relatively isolated regions (modules) along the few links connecting them to the rest of the network. If for such a module these links would be hypothetically cut one by one, the corresponding eigenvalue would gradually increase towards unity, $`\lambda ^{(\alpha )}1`$, while the eigenvector would become more and more localized on the module. When finally the module is completely disconnected from the network the eigenvector has evolved to the steady state solution on the module, which has the participation ratio equal to its size. Thus the PR of slowly decaying eigenmodes serves as a good quantitative estimate of the size of modules in the network. In an RSFN these modules are small consisting of just a handful of nodes that accidentally happen to be loosely connected to the rest of the network. The fact that the participation ratios of slow decaying modes on the Internet significantly exceed those in an RSFN indicates that the corresponding modules are real and not accidental. The average participation ratio of slowly decaying modes can be quantified by $`_\alpha \text{PR}^{(\alpha )}|\lambda ^{(\alpha )}|^k/_\alpha |\lambda ^{(\alpha )}|^k`$. For $`5k10`$ this average changes only slowly in both the Internet and an RSFN and equals approximately $`60`$ and $`5`$, respectively. A rough estimate for the number of different modules is given by the number of slowly decaying non-oscillatory states in Fig. 1 A that have a participation ratio significantly exceeding that of an RSFN. For the Internet the number is around $`100`$. The sum of the participation ratios for these first 100 modes, $`5400`$ , is a rough estimate of the total number of nodes in the modular part of the network. This should be compared to the same sum being approximately equal to $`520`$ in an RSFN. If one takes special care to avoid double counting nodes that appear more than once among the set of $`PR^{(\alpha )}`$ nodes with the largest $`|c_i^{(\alpha )}|`$ taken for each eigenmode $`1<\alpha 100`$, this number gets reduced to $`1800`$. Thus the overall modularity of the Internet network is at least $`1800/650030\%`$.

To determine the organizing principle behind these Internet modules in Fig. 2 we plot the outgoing current $`c_i^{(2)}`$ in the slowest decaying diffusion mode ($`\lambda ^{(2)}=0.9626`$) as a function of the AS number (note that some AS numbers are not yet in use).

Autonomous Systems known to be located in Russia are marked with a circle. The PR for this eigenmode is $`107`$, while the total number of Russian AS in our dataset is 174. In Fig. 2 one can see that almost all the Autonomous Systems that significantly participate in this mode (large positive $`c_i^{(2)}`$) are Russian. We have checked that the few exceptions to this rule are in fact Autonomous Systems closely related to Russia. Thus in the slowest decaying mode the diffusion current flows from a module that may be identified with Russia towards the rest of the Internet. Curiously enough the set of Autonomous Systems furthest away from Russia (the most negative $`c_i^{(2)}`$) are located in the US and belong to the US Military. This possible legacy of the cold war makes Russia and the US Military the extreme edges of the Internet. Performing a similar analysis for other slowly decaying modes we get a similar picture, just with other pairs of countries being pulled out. For the Internet the modules thus correspond to individual countries, or for large countries – to organizational or geographical features within the country.

It is interesting to note that these country-modules cannot be detected using the spectral analysis of the adjacency matrix of the network SpectraOfGraphs ; FarkasBarabasiSpectra ; GohKahngKimSpectra . The elements of this matrix, closely related to $`𝐓`$, are equal to 1 for a pair of neighboring nodes and 0 otherwise. The largest eigenvectors of the adjacency matrix are known to be localized primarily on the highest connected hubs and their neighbors FarkasBarabasiSpectra ; GohKahngKimSpectra . However, unlike in the case of $`𝐓`$, this undesirable localization cannot be properly eliminated simply by dividing the components of the eigenvectors by the connectivity $`K_i`$. Hence eigenvectors of the adjacency matrix do not properly reflect the country-based modular structure uncovered in this work.

Having established that the Internet is indeed modular we now address the question of how good these individual modules are. To this end we compare different eigenmodes $`c_i^{(\alpha )}`$ to each other. Although the primary feature in a slowly decaying eigenmode $`c_i^{(\alpha )}`$ is the flow between a dominant pair of country-modules (such as between Russia and the US Military in $`c_i^{(2)}`$), other modules may also participate in it but to a smaller extent. This gives rise to a fine structure within slow decaying modes that is not captured by the participation ratio. The hallmark of a good module is that even though it participates in different eigenmodes $`\alpha `$ to a different extent, the relative distribution of $`\alpha `$-currents within the module stays approximately the same. In other words it enters different eigenmodes as just one degree of freedom. In this case the ratio $`c_i^{(\alpha )}/c_j^{(\alpha )}`$ is approximately independent of the eigenmode $`\alpha `$ for any pair of nodes $`i`$ and $`j`$ within the module. This is equivalent to the condition that for any two different eigenmodes $`\alpha `$ and $`\beta `$, $`c_i^{(\alpha )}/c_i^{(\beta )}=\text{const }`$ for every node $`i`$ belonging to the given module. In Fig. 3 we plot the outgoing currents in the two slowest decaying non-oscillatory eigenmodes – $`c_i^{(2)}`$, and $`c_i^{(3)}`$ – as a function of each other. Similar plots can be made for other pairs of slow decaying eigenmodes.

The principal feature in this kind of plot is a star like shape, where different rays of the star correspond to individual country-modules. This type of plot is more powerful in identifying individual modules than the participation ratio alone. Indeed, in Fig. 3 one can easily detect not only the most excited modules like Russia, France, and US (red squares, green circles, and blue crosses), but also less excited ones like Korea (orange triangles). We believe that the idea of measuring the quality of individual modules by how proportionally their nodes participate in different slowly decaying modes, can be easily generalized to other dynamical processes taking place on the network such as e.g. spin dynamics, vibrational modes, etc.

Finally, we would like to point out another interesting feature of Fig. 1. Both the density of states and the participation ratio are nearly symmetric around $`\lambda =0`$ for both the Internet and an RSFN. This near symmetry indicates that both these networks are almost bipartite Comment1 , a feature also observed in citation networks, but not in metabolic networks Bilke2002 . In fact, a more detailed analysis shows that while almost every slow oscillatory mode ($`\lambda ^{(\mathrm{osc})}1`$) is related to the corresponding non-oscillatory mode with $`\lambda |\lambda ^{(\mathrm{osc})}|`$, the reverse is not true as there are roughly 30% more modes near $`\lambda =1`$ than near $`\lambda =1`$. Those country-modules that are present in both the oscillatory and non-oscillatory parts of the spectrum are internally almost bipartite. The simplest bipartite graph is a tree and this seems to be the dominant structure within the Internet modules.

In summa, we have demonstrated how a diffusion process taking place on the Internet network allows one to extract information about its modules and extreme edges. For many “real-world” complex networks the local context of a node (in terms of the linkage pattern) reflects, or, perhaps, even determines the importance and function of the given node. For instance, in biology one can successfully assign putative functions to unclassified proteins based on the function of their interaction partners SchwikowskiUetz . In general the diffusion process introduced in this work can be seen as a systematic way to explore the local linkage structure of a network beyond just the nearest neighbors. The detection of the modular structure of a network is just one possible application of such a process google ; Bilke2002 .

Work at Brookhaven National Laboratory was carried out under Contract No. DE-AC02-98CH10886, Division of Material Science, U.S. Department of Energy. Two of us (K.E and K.S.) thank the Institute for Strongly Correlated and Complex Systems at Brookhaven National Laboratory for financial support during visits when part of this work was completed.