text:
# 1 INTRODUCTION

## 1 INTRODUCTION

Galactic isolated neutron stars (NSs) are expected to be as many as $`10^8`$$`10^9`$, a non–negligible fraction of the total stellar content of the Galaxy. In the last 10 years they received special attention. The idea of observing such objects in the X-ray range has emerged in soon after discovery of the first X-ray sources (Ostriker et al. 1970). Now we are sure that few dim ROSAT sources are nearby isolated neutron stars (see Neuhäuser 2001, Walter 2001).

All these sources emit a thermal spectrum at $`100`$ eV and the derived column densities place them at relatively close distances ($`150`$ pc) with luminosities $`L10^{30}`$$`10^{31}`$ erg $`\mathrm{s}^1`$. Several of them can be visible due to accretion from the interstellar medium (ISM), others — due to thermal emission of cooling relatively young, $`1`$ Myr, NSs (see Treves et al. 2000 for a recent review of accreting isolated NSs, and Popov et al. 2000b for discussion on relative fraction of cooling and accreting NSs in the solar vicinity).

X-ray observations of globular clusters (GCs) show large population, $`>30`$, of dim X-ray sources (hereafter dim GC sources) with $`L10^{31}`$$`10^{34}`$ erg $`\mathrm{s}^1`$ (Verbunt 2001). Pfahl & Rappaport (2001) suggested, that part of dim GC sources can be old isolated accreting NSs. In this short note we try to analyze this idea from the point of view of population synthesis (see Lipunov et al. 1996 for detailed description of population synthesis method).

Population synthesis of isolated NSs in the Galactic disk was made previously by several authors (see Popov et al. 2000a,b for recent results and references). Evolution of isolated NSs in the Galactic disk and in GCs must be significantly different. Here we try to explore with a simple model evolution of progenitors of isolated accreting NSs in GCs. In the next section we briefly discuss our models. In section 3 we present our results, and in the last section we discuss properties of old accreting NSs in GCs.

## 2 THE MODEL

An isolated NS can pass through four main stages of evolution: Ejector, Propeller, Accretor and Georotator (see Lipunov 1992 for detailed description of each stage). We consider, that all NSs are born as Ejectors, and then while they slow down (or/and their magnetic field decays) they pass through the Propeller and Accretor stages. If spatial velocity or/and magnetic field of a NS are very high, then this NS can appear as Georotator after the end of the stage of ejection.

For constant magnetic field we can estimate the period of time, which a NS spends on the stage of Ejector, $`t_E`$, as (see for example Popov & Prokhorov 2000):

$$t_E10^9\mathrm{yrs}\mu _{30}^1n^{1/2}v_{10},$$
(1)
where $`\mu _{30}=\mu /10^{30}`$ G cm3 — magnetic moment, $`v_{10}=v/10`$ km s-1 — velocity relative to ISM, $`n`$ — ISM number density.

For the Propeller stage one can find in the literature a lot of different spin-down regimes, but in the case of constant magnetic field the Propeller stage is always considered to be much shorter than the Ejector stage (see Lipunov & Popov 1995). Taking this into account we neglect the time spent by NSs on the Propeller stage in our simple estimates (for field decay this stage become very important! see Livio et al. 1998, Colpi et al. 1998, Popov et al. 2000a).

We calculate relative fraction of Accretors (i.e. objects with $`t_E<t_{Hubble}=10^{10}\mathrm{yrs}`$ which cannot appear as Georotators due to relatively low magnetic field or/and spatial velocity) among other isolated NSs for typical disk and GC conditions after $`10^{10}`$ yrs of evolution. In the case of the Galactic disk we use constant starformation rate, and in the case of GC starformation rate is considered as a $`\delta `$-function (i.e. a starformation burst at $`t=0`$).

For both cases we use Maxwellian velocity distribution with $`\sigma =140`$ km s-1, and log-normal distribution for magnetic fields of NSs typical for radiopulsars (see Popov et al. 2000a,b for discussion of these choices). For GC velocity was truncated at 30 km s-1, which is a typical escaping velocity from the cluster. For the Galactic disk we use typical value of ISM number density $`n=1`$ cm-3, for GCs — $`n=100`$ cm-3 (see Pfahl & Rappaport 2001 and short discussion below). In our simple model $`n`$, $`\mu `$ and $`v`$ are considered to be constant during the whole evolution.

For both populations (disk and GCs NSs) we assume the same parameters of NSs (initial spin periods, moments of inertia, masses etc.). Initial spin periods are assumed to be equal to 0.020 sec for all NSs in our calculations. Our previous estimates (Popov et al. 2000a) show that small variations of this parameter do not change final results for Accretors.

## 3 CALCULATIONS AND RESULTS

Our aim is to calculate the number of accreting isolated NSs per typical GC with the total mass $`10^5M_{}`$ (we neglect the fact, that the GC mass slightly decreases during GC evolution).

As the first step we compare fractions of Accretors in GCs and in the Galactic disk. To do it we just run our simple models (described in the previous section) for the same numbers of NSs for each case. An important result is the following: for GCs fraction of Accretors relative to other stages is about 26 times lower than in the Galactic disk. Of course it is just a rough estimate, but it shows, that even for high ISM density ($`n=100`$ cm-3) due to escaping of the most part of population of isolated NSs number of isolated accreting NSs in GCs is small. Escaped NSs are considered as Ejectors as far as they have high spatial velocities and move in low-density surrounding. One can compare typical escaping velocity, 30 km s-1, with typical velocity of Accretors in the Galactic disk (Popov et al. 2000b), $`>50`$ km s-1.

Then we try to estimate the number of Accretors for a typical GC. In the Galactic disk Accretors form about 1% of the total NS population (Popov et al. 2000a). Mass of the disk population of stars is about $`710^{10}M_{}`$. Typical mass of a GC is assumed to be $`10^5M_{}`$. We note, that for the Galactic disk fraction of 1% was calculated (Popov et al. 2000a) for constant starformation rate, for GCs we calculate it for a starformation burst (age of the burst is equal to the age of the disk).

In the disk we expect about $`N=10^9`$ isolated NS. I.e. about $`10^7`$ Accretors. After simple calculations we obtain about $`(10^90.0110^5)/(26710^{10})0.55`$ Accretors per typical GC ($`M=10^5M_{}`$). If total number of NSs in the disk is lower, than our estimate should be lower too.

We can also calculate the number of Accretors per GC in a more simple way. Let us take the Salpeter mass function, and assume, that NSs are born from mass interval 10–40 solar masses (Lipunov et al. 1996). In a simple model described above (starformation burst, $`10^{10}`$ years of evolution, constant parameters of the ISM and NSs during their evolution) we can calculate that about 0.25% NSs become Accretors in GCs (mainly because most part of NSs leave GCs due to high kick velocity, but nearly all NSs which stay bounded with the GC become Accretors). It gives us about 1.1 Accretor per GC. It can be considered as an upper limit on the number of Accretors, as far as for GCs number of massive stars can be lower (initially) than it follows from the Salpeter mass function, and evolution of NSs can proceeds with an average ISM density lower than the value $`n=100`$ cm-3 which we used in our calculations: all these factors decrease the number of Accretors.

So, we can roughly estimate the number of Accretors per GC with the total mass $`10^5`$ solar masses as 0.5-1. We also note, that coincidence of the two estimates can be considered as an additional (but not very strong) argument for the number of Galactic NSs equal to $`10^9`$.

To summaries, simple evolutionary estimates do not contradict the idea, that at least part of dim X-ray sources in GCs can be old accreting isolated NSs. It is important to repeat such calculations in more details (spin evolution, dynamical evolution, realistic ISM distribution etc.) and with realistic models of field decay.

## 4 DISCUSSION

It is reasonable to discuss in brief the following important topics:

1. Interstellar gas density in GCs.

2. Observed temperatures of dim sources in GCs.

3. Periods of accreting NSs in GCs (no pulsations were detected from dim GC sources) and rotational equilibrium (period changes).

4. Magnetic field distribution for Accretors in GCs.

Pfahl & Rappaport in their paper suggested and used high value of the gas number density in GCs: $`100`$ cm-3. We used the same value in our calculations. Actually, nobody observed gas in GCs at such high density. But it is reasonable to expect these high values and observations provide high upper limits: $`(50`$$`100)`$ cm-3 (Knapp et al. 1996).

Average density along the track during a NS evolution can be lower than the high value we used. This effect will increase the duration of the Ejector stage (see eq.1), and correspondently decrease the number of Accretors. But if we a priori consider dim X-ray sources to be Accretors, we have to take high number density in order to explain observed luminosities (Pfahl & Rappaport 2001). The same high density comes from estimates of gas accumulation in GCs due to mass loss from the stars (Knapp et al. 1996, Pfahl & Rappaport 2001).

That’s why we argue, that the effect of decreasing of the number of Accretors is not very significant. Otherwise, in the case of low density the hypothesis of accreting isolated NSs should be rejected due to difficulties with the explanation of X-ray luminosity of the sources under discussion.

Relatively high temperatures of dim GC sources (0.1-0.5 keV), which is higher than for dim ROSAT INS candidates ($`<200`$ eV with typical values 50-100 eV, see Treves et al. 2000, Walter 2001, Neuhäuser 2001) are in favour of the accreting hypothesis (both for single and binary NSs). For example, young cooling NSs, age $`1`$ Myr, should be cooler than Accretors in the case of polar cap accretion. In Popovet al. 2000b the authors plot temperature distribution for Accretors (without magnetic field decay) for the disk population. Maximum of the distribution corresponds to $`400`$$`500`$ eV, the total range is about 100–1000 eV (depending on the magnetic field and accretion rate). But as far as luminosities of dim GC sources are higher than luminosities of the disk accreting sources (the ISM density is higher), they should be also slightly warmer. May be accretion in that sources proceeds onto very large polar cap due to magnetic field decay (objects are very old).

Old Accretors in disk should have relatively long periods if their magnetic field does not decay significantly (Lipunov & Popov 1995). But in GCs spin periods for the same magnetic field should be much shorter, as far as the ISM density is higher:

$$p_A42\mu _{30}^{6/7}v_{10}^{9/7}n_{100}^{3/7}\mathrm{s}.$$
(2)
Here $`p_A`$ is a critical period for transition from the stage of Propeller to the stage of accretion (Lipunov 1992). On that stage a NS continues to spin down till it comes to rough equilibrium with turbulized ISM. So, typical periods should be close to some value,

$$p_{turb}(100300)v_{10}^{13/3}n_{100}^{2/3}\mu _{30}^{2/3}\mathrm{s},$$
(3)
for $`n=100`$ cm-3 (Lipunov & Popov 1995, Konenkov & Popov 1997), slightly longer than $`p_A`$. Realistic period distribution for accreting isolated NSs will be presented elsewhere (Popov, Prokhorov & Khoperskov 2001). It is characterized by broad maximum near $`p_{turb}`$, slightly shifted to longer periods.

In the picture described above spin periods of isolated accreting NSs in GCs can be observed, and lack of period detection is not in favour of Pfahl & Rappaport (2001) hypothesis. To solve this problem it is necessary to introduce something non-standard for NSs in GCs. For example significant field decay, so that Alfven radius will be about the size of a NS, and no pulsations can be observed in that case due to low level of modulation. As we note above it can also explain relation between temperature and luminosity of that sources (due to large polar cap area).

Periods and their derivatives (if detected) should show significant fluctuations on the time scale $`R_G/vv^35.9\mathrm{yrs}v_{10}^3`$. On the same time scale one expects variations of luminosity, and they are observed (Pfahl & Rappaport 2001) on the time scale $`1`$ year. Detection of spin periods and analysis of their correlation with luminosity changes are very desirable.

Magnetic fields of Accretors in GCs should be, on average, smaller than in disk, because it is much easier to become Accretor in a GC (density is higher, and typical spatial velocity is lower: in disk it is about 50 km s-1, and in GCs about 20 km s-1). For lower magnetic fields it is more difficult to reach $`p_A`$ (see eq.2), so there should be low values of the magnetic field for which only in high density environment of GC it is possible for a NS to become an Accretor.

We also note, that field decay can both: decrease and increase number of Accretors (Colpi et al. 1998, Livio et al. 1998, Popov et al. 2000a, Popov & Prokhorov 2000). So, significant field decay, which make it easier for a NS to reach the stage of accretion, can compensate the effect of lower average ISM density along the evolutionary path of the NS.

In principle some of dim GC sources can also be Georotators accreting warm ISM (aka MAGACs, see Rutledge 2001). These objects should be relatively hard sources, without any periodicity on the timescale $`<10^5`$ s. But for GC conditions very large magnetic moments are required to become a Georotator, because in dense ISM ($`n100`$ cm-3) among low velocity isolated NSs, which can stay in the GC, only most magnetized objects can become Georotators instead of becoming Propellers (and then — Accretors). It makes this hypothesis unlikely. If one try to explain dim GC sources as young cooling NSs, then it is necessary to accept, that they are really very young to be so hot and luminous. It is difficult to explain why in several GCs young NSs are observed (Andrei Zakharov noted to us, that there is a possibility of NS formation in GCs due to recent star formation from the accumulated matter, which was lost by cluster members on their late stages of evolution). Most probably dim GC sources are accreting isolated NSs with decayed magnetic fields in correspondence with Pfahl & Rappaport (2001) hypothesis.

### Acknowledgements

We thank Andrei Zakharov for discussions and unknown referees for useful comments. This work was supported by the RFBR (01-02-06265), and NTP “Astronomy” (1.4.4.1; 1.4.2.3) grants.

### References

Colpi, M., Turolla, R., Zane, S., and Treves, A. (1998) ApJ 501, 252.

Knapp, G.R., Gunn, J.E., Bowers, P.F., and Vasques Poritz, J.F. (1996) ApJ 462, 231.

Konenkov, D.Yu., and Popov, S.B. (1997) Astronomy Letters 23, 498.

Lipunov, V.M. (1992) Astrophysics of Neutron Stars Berlin: Springer & Verlag.

Lipunov, V.M., and Popov, S.B. (1995) AZh 71, 711.

Lipunov, V.M., Postnov, K.A. and Prokhorov, M.E. (1996) Astron. and Space Phys. Rev. 9, part 4.

Livio, M., Xu, C., and Frank, J. (1998) ApJ 492, 298.

Neuhäuser, R. (2001) Astron. Nachr. 322, 3 (astro-ph/0102004).

Ostriker, J.P., Rees, M.J. and Silk J. (1970) Astrophys. Letters 6, 179.

Pfahl and Rappaport (2001) ApJ 550, 172 (astro-ph/0009212).

Popov, S.B., Colpi, M., Treves, A., Turolla, R., Lipunov, V.M., and Prokhorov, M.E. (2000a) ApJ 530, 896.

Popov, S.B., Colpi, M., Prokhorov, M.E., Treves, A., and Turolla, R. (2000b) ApJ 544, L53.

Popov, S.B., and Prokhorov, M.E. (2000) A&A 357, 164.

Popov, S.B., M.E., and Khoperskov, A.V. (2001) (in preparation).

Rutledge, R.E. (2001) ApJ (in press), astro-ph/0101550.

Treves, A., Turolla, R., Zane, S., and Colpi, M. (2000) PASP 112, 297.

Verbunt, F. (2001) A&A 368, 137 (astro-ph/0012261).

Walter, F. (2001) ApJ 549, 433 (astro-ph/0009031).


format: text

text:
# Viscoelasticity near the gel-point: a molecular dynamics study

## Abstract

We report on extensive molecular dynamics simulations on systems of soft spheres of functionality $`f`$, i.e. particles that are capable of bonding irreversibly with a maximum of $`f`$ other particles. These bonds are randomly distributed throughout the system and imposed with probability $`p`$. At a critical concentration of bonds, $`p_c0.2488`$ for $`f=6`$, a gel is formed and the shear viscosity $`\eta `$ diverges according to $`\eta (p_cp)^s`$. We find $`s0.7`$ in agreement with some experiments and with a recent theoretical prediction based on Rouse dynamics of phantom chains. The diffusion constant decreases as the gel point is approached but does not display a well-defined power law.

The behavior of transport coefficients and elastic moduli near the gelation transition has been discussed in the literature for many years. To date no consensus on either the theoretical or experimental side has emerged as far as the critical behavior of these quantities is concerned. The phenomenology is as follows. As monomers or polymers are randomly crosslinked to each other in a melt, the shear viscosity $`\eta `$ increases with crosslink concentration $`p`$ and diverges at a critical concentration $`p_c`$ at which an amorphous rigid network is formed. Experiment and theory both yield $`\eta (p_cp)^s`$ but there is no general agreement regarding the value of the exponent $`s`$. Indeed there is good reason to expect at least two different universality classes: As de Gennes showed, vulcanization (crosslinking of very long chains) must be distinguished from gelation (crosslinking of short chains or monomers) as far as critical behavior is concerned. This conclusion is supported by recent experiments which show quite clearly that chain length is a relevant parameter.

Before describing our model and calculations, we discuss briefly the experimental and theoretical situation for gels putatively in the percolation (short chains) universality class. One group of experiments has produced exponent values for the shear viscosity in the range $`0.6s0.9`$. Another group  has reported values for $`s`$ in the range 1.1 — 1.3 and interpreted  these in terms of a Rouse model without hydrodynamics. On the theory side, de Gennes argued that the viscosity is analogous to the conductivity of a random mixture of normal conductors and superconductors, with an exponent $`s0.67`$. The aforementioned Rouse model prediction is $`s=2\nu \beta 1.35`$, where $`\nu 0.88`$ and $`\beta 0.41`$ are the correlation length and order parameter exponents of percolation theory in three dimensions. Finally, a recent theoretical analysis of a different model with Rouse-dynamics has predicted $`s=\varphi \beta 0.7`$ where $`\varphi 1.11`$ is the crossover exponent of a random resistor network. Given this wide disparity in both theoretical and experimental results, computer simulations may help to clarify the situation.

The shear modulus $`\mu `$ of a rigid network near the gel point is typically entropic in nature and it vanishes with a power law of its own as the gel point is approached from the rigid phase $`\mu (pp_c)^t`$. Recent numerical work on systems in the percolation universality class has provided evidence that $`t=f`$ in both two and three dimensions, where $`f`$ is the exponent that describes the critical behavior of a randomly disordered network of conductors and insulators near the percolation point. This result is consistent with another argument of de Gennes. In the dynamical scaling theory of the gelation transition, the two exponents $`t`$ and $`s`$ are not independent, but rather obey the sum rule $`s+t=z`$ where $`z`$ describes the divergence of the longest relaxation time in the incipient gel: $`t^{}=t_0(p_cp)^z`$. This connection allows an important consistency check between the results reported here and those of .

The model that we simulate is capable of describing the entire range from simple liquid to entropic solid. All particles interact through a soft sphere potential $`V_{sa}(r_{ij})=ϵ(\sigma /r_{ij})^{36}`$  with $`\sigma =1`$ and, for our simulations, $`k_BT/ϵ=1`$. If there are no other interactions, this system forms a simple three-dimensional liquid at least at low density. All of our simulations are done at a volume fraction $`\mathrm{\Phi }=\pi N\sigma ^3/6V=0.4`$ which is well away from the liquid-solid coexistence density. The viscosity of the system is progressively increased by introducing random crosslinks between particles. Specifically, the system of particles is initially placed on a simple cubic lattice that fills the cubic computational box. Each particle may bond with probability $`p`$ with each of its six nearest neighbors. This bonding is permanent and enforced by the spherically symmetric potential $`V_{nn}(r_{ij})=\frac{1}{2}k(r_{ij}r_0)^2`$ with $`k=5ϵ/\sigma ^2`$ and $`r_0=(\pi /6\mathrm{\Phi })^{1/3}\sigma `$. For $`p<p_c`$, where $`p_c0.2488`$ is the bond percolation probability of the simple cubic lattice, the system consists of finite clusters of varying masses. For $`p>p_c`$ the system is an entropic solid with nonvanishing shear modulus. This system has been previously studied by us for $`V_{sa}=0`$ and $`p>p_c`$. Farago and Kantor and Cohen and Plischke have shown that self-avoidance is irrelevant as far as the critical behavior of the elastic constants is concerned.

In the simulations, the system of particles is first equilibrated for $`5\times 10^4`$ time steps with Brownian dynamics. At the end of this equilibration time, the damping and thermal noise are turned off and the subsequent evolution is conservative. The equations of motion are integrated with a standard velocity Verlet algorithm  with a time step of $`\delta t=0.005\sqrt{m\sigma ^2/ϵ}`$. The shear viscosity $`\eta (p)`$ and the self-diffusion constant $`D(p)`$ are then calculated from the appropriate Green-Kubo formula :

$`\eta `$ $`=`$ $`\underset{t_{max}\mathrm{}}{lim}{\displaystyle _0^{t_{max}}}𝑑tC_{\sigma \sigma }(t)`$ (1)
$`=`$ $`\underset{t_{max}\mathrm{}}{lim}{\displaystyle \frac{1}{3Vk_BT}}{\displaystyle _0^{t_{max}}}𝑑t{\displaystyle \underset{\alpha <\beta }{}}\sigma ^{\alpha \beta }(t)\sigma ^{\alpha \beta }(0)`$ (2)
where

$`\sigma ^{\alpha \beta }={\displaystyle \underset{i=1}{\overset{N}{}}}mv_{i\alpha }v_{i\beta }{\displaystyle \underset{i<j}{}}{\displaystyle \frac{r_{ij\alpha }r_{ij\beta }}{r_{ij}}}V^{}(r_{ij})`$

and where $`V^{}`$ is the derivative of the total potential energy. Similarly, the diffusion constant is given by the familiar expression:

$$D=\underset{t_{max}\mathrm{}}{lim}\frac{1}{3N}\underset{i=1}{\overset{N}{}}_0^{t_{max}}𝑑t𝐯_i(t)𝐯_i(0).$$
(3)
It is well known that the velocity-velocity correlation function decays extremely slowly, typically with a ‘long time tail’ $`t^{3/2}`$ power law, even in simple liquids . We find the same long-time behavior in our simulations as well, remarkably for all values of $`p`$. There is considerably more disagreement regarding the stress-stress correlation function. Extended mode-coupling theory  suggests that close to the melting point the stress correlator decays exponentially at long times. Powles and Heyes  have found that both an exponential decay and a Lorentzian provide a reasonably good fit to their data in the dense liquid regime. In our situation where clusters of various sizes form the system, the decay is dramatically affected by crosslinking and becomes very slow close to the gel-point. Therefore in the evaluation of (1,3) we have used time series from $`t_{max}=300\tau `$ to $`t_{max}=1200\tau `$ where $`\tau `$ is the average time between collisions for particles in the uncrosslinked liquid ($`p=0`$). Even with such long runs, for $`p`$ close to $`p_c`$ an estimate of the residual integral to $`t=\mathrm{}`$ had to be added. This is discussed further below.

We have simulated systems consisting of $`N=L^3`$ particles with $`L=5,8,12`$ and $`20`$, the first three over the concentration range $`0p0.24`$, the fourth only for $`p0.20`$. For these relatively small systems, the probability that one of the clusters percolates in at least one of the three directions is an issue. If there is percolation, then the viscosity is not a well-defined quantity and the sample is characterized instead by a shear modulus. Therefore, we have eliminated all percolating samples from the calculation. For $`L=5`$, a non negligible number of samples percolates in at least one direction at $`p=0.15`$; at $`L=20`$ percolation becomes significant only at $`p=0.23`$. Figure 1 depicts the stress correlator for $`N=12^3`$ particles for a set of crosslink concentrations. The top panel shows this function over the entire range in time, the lower panel for $`t>2\tau `$. The top panel shows that the changes for short times between the liquid ($`p=0`$) and the incipient gel ($`p=0.24`$) are extremely small. The effect of increasing crosslink density is illustrated in the lower panel where it is clearly seen that the decay of the correlation function becomes progressively slower as the gel point is approached and that even at the longest time, for $`p0.2`$ the correlator is nonnegligible. Therefore, integrating $`C_{\sigma \sigma }(t)`$ only to $`t_{max}`$ would result in an underestimate of the shear viscosity. In order to capture the remaining contribution, we have fit $`C_{\sigma \sigma }(t)`$ with a stretched exponential over various windows $`(t_s,t_{max})`$ for starting values $`t_s>2\tau `$. Such a fit is shown in Figure 1b for the uppermost curve ($`p=0.24`$) over the range $`2\tau tt_{max}`$. The fit is essentially indistinguishable from the simulation. However, there is no fundamental reason to believe that a stretched exponential is the functional form of the long-time behavior of $`C_{\sigma \sigma }(t)`$. This is the reason that we have chosen several different starting points for the fit: the spread in values of the remaining integral of the fitting function from $`t_{max}`$ to infinity provides an estimate of the error associated with this part of the calculation.

In Figure 2 we display our data for the dimensionless shear modulus $`\eta \sigma ^2/\sqrt{mk_BT}`$ in two versions. In part (a) the raw data is shown as function of $`(p_cp)`$ together with a guide to the eye of the form $`a(p_cp)^{0.7}`$. This function clearly captures the general behavior of the data in the intermediate range of $`p`$. For $`p`$ close to zero, one would not expect the system to anticipate the formation of a gel at $`p0.25`$ and for $`p`$ close to $`p_c`$ finite-size effects are clearly evident. Part (b) of this figure attempts to collapse the data by means of the finite-size scaling ansatz $`\eta (L,p)=L^{s/\nu }\mathrm{\Psi }[L(pp_c)^\nu ]`$ with $`s=0.7`$ and $`\nu =0.876`$. Internal consistency requires that the scaling function have the asymptotic form $`\mathrm{\Psi }(x)x^{s/\nu }`$ for large $`x`$ and a line corresponding to this form is also shown on the figure.

Further support for the conclusion $`s0.7`$ comes from the complex frequency-dependent viscosity $`\eta ^{}(\omega )\eta ^{}(\omega )+i\eta ^{\prime \prime }(\omega )=G^{}(\omega )/(i\omega )`$ where $`G^{}`$ is the complex elastic modulus. The scaling ansatz for these functions  is $`lim_{\omega 0}\eta ^{}(\omega )(p_cp)^s`$ for $`p<p_c`$; $`lim_{\omega 0}G^{}(\omega )(pp_c)^t`$ for $`p>p_c`$ and, for frequencies $`\omega >\omega ^{}`$, where $`\omega ^{}`$ is a characteristic crossover frequency that approaches zero as $`pp_c`$, $`G^{}(\omega )(i\omega )^u`$. The connection between the critical behavior of the modulus in the rigid phase and the viscosity in the fluid phase is then expressed through the scaling relation $`u=t/(s+t)`$. Moreover, in the high frequency region, one expects $`\eta ^{}(\omega )`$ and $`\eta ^{\prime \prime }(\omega )`$ to both vary as $`\omega ^{u1}`$ and the ratio of the real and imaginary parts to obey $`u=2/\pi \mathrm{tan}^1\{\eta ^{}/\eta ^{\prime \prime }\}`$. In our previous work on the rigid phase , we have concluded that $`t2`$ in three dimensions. Therefore, with $`s0.7`$ we have $`u0.74`$. The frequency-dependent viscosity is plotted in Fig. 3 for $`L=12`$ and $`p>0.20`$. There is clearly a region of power-law behavior that extends to lower frequencies as the critical point is approached. This behavior is seen more clearly in $`\eta ^{}`$ than in $`\eta ^{\prime \prime }`$. Nevertheless, both pieces of the shear viscosity decrease in a way consistent with $`\omega ^{0.27}`$ in very satisfactory agreement with the foregoing analysis. As well, the ratio of $`\eta ^{}`$ to $`\eta ^{\prime \prime }`$ in this regime produces a second estimate $`u0.76`$.

Finally, in Figure 4 we show the dimensionless self-diffusion constant $`\sqrt{m/\sigma ^2k_BT}D`$ obtained from the velocity-velocity correlator (3) for all particles in the system. In this case, we show only the raw data. It seems clear from the behavior of $`D`$ in the critical region that a finite-size scaling analysis is unlikely to improve the collapse of the data. For $`p<0.2`$ the data are not inconsistent with a power law behavior of the form $`(p_cp)^{0.7}`$ but the evidence for this is weak at best. Moreover, the fact that the data for $`pp_c`$ are essentially independent of $`L`$ suggests that $`D(L\mathrm{},pp_c)`$ is finite. Precisely at the gel-point, in the thermodynamic limit, the system consists of a percolating cluster with fractal dimension $`D_F2.5`$. The particles that are not on the infinite cluster are organized into finite clusters of various sizes. Approximately $`18\%`$ remain as monomers that presumably are able to diffuse quite easily through the percolating cluster since this cluster contains holes on all length scales. This would account for the absence of critical behavior in $`D`$.

We are aware of one previous simulation that attempted to address the critical behavior of the shear viscosity near the gel point. Recently del Gado et al. simulated a very different model, namely particles confined to the sites of a lattice and randomly crosslinked to form clusters of various sizes. This system was then evolved by a bond fluctuation method and the diffusion constants $`D_m`$ of clusters of mass $`m`$ measured. They postulated the relation $`D(R)R^{(1+s/\nu )}`$, where $`R`$ is the radius of gyration of a cluster, and from this determined $`s1.3`$. We are not aware of any rigorous derivation of this connection between diffusion and viscosity. However, it may be that their model simply contains different physics.

In conclusion, we have obtained the shear viscosity of a viscoelastic liquid below the gel point for a conceptually simple model that we have previously studied in the rigid phase. We have presented evidence that the shear viscosity diverges with an exponent $`s0.7`$ consistent with recent theoretical work and some of the experimental data . Future extensions of this work will include an investigation of viscoelasticity in two dimensions and a study of dynamics in the gel phase.

One of us (BJ) thanks the Physics Department at Simon Fraser University, where this work was carried out, for its hospitality. This work was supported by the NSERC of Canada.


text:
# A Model for Neutrino Warm Dark Matter and Neutrino Oscillations

## Abstract

The muon- and tau-neutrinos with the mass in the keV range, which are allowed in a low reheating temperature cosmology, can compose the warm dark matter of the universe. A model of four light neutrinos including the keV scale $`\nu _\mu `$ and $`\nu _\tau `$ is studied, which combines the seesaw mechanism and the Abelian flavor symmetry. The atmospheric neutrino anomaly is due to the $`\nu _\mu \nu _\tau `$ oscillation. The solar neutrino problem is answered by the oscillation into the light sterile neutrino, where the SMA, LMA, and LOW-QVO solutions can be accommodated in our scenario.

preprint:   ASITP-2001-002 KIAS-P01007 hep-ph/0102246

Cosmological studies show increasing evidences that the dark matter is in fact warm, neither cold nor hot, so as to explain the observed structure and behavior of the expanding universe . In particle physics, as the candidates of the warm dark matter (WDM), gravitinos or sterile neutrinos with masses in the keV range have been proposed. In this paper, we consider the possibility that the ordinary active neutrinos are the WDM.

Experiments have provided various constraints on the neutrino mass spectrum. From the direct experimental search, we have $`m_{\nu _1}<2.5`$ eV, $`m_{\nu _2}<170`$ keV and $`m_{\nu _3}<15.5`$ MeV . More information comes from the neutrino oscillation experiments. The Super-Kamiokande (Super-K) data for the atmospheric neutrino anomaly suggest that the $`\nu _\mu `$ is maximally mixed with $`\nu _x`$ ($`xe`$) with $`\mathrm{\Delta }m_{\mu x}^22.2\times 10^3`$ eV2 . And the $`x=\tau `$ case is strongly favored. The solar neutrino deficit  may imply the oscillation of $`\nu _e`$ into $`\nu _y`$. There are several allowed parameter regions. For example, if the $`\nu _y`$ is a sterile neutrino, the currently favored solution is the small mixing angle (SMA) with $`\mathrm{\Delta }m_{ey}^25\times 10^6`$ eV2 and $`\mathrm{tan}^2\theta _{ey}10^3`$. If one experiment among Super-K, Ga and Cl is removed from the data analysis, the large mixing angle (LMA) solution with $`\mathrm{\Delta }m^210^510^4`$ eV2 and the low-mass and quasi-vacuum oscillation (LOW-QVO) solution with $`\mathrm{\Delta }m^210^{10}10^7`$ eV2 are also allowed. The LSND experiment has reported positive appearance results of $`\overline{\nu }_\mu \overline{\nu }_e`$ oscillations, which implies $`\mathrm{\Delta }m_{e\mu }^21`$ eV2 and $`\mathrm{sin}^2(2\theta _{e\mu })10^2`$ . However, a large part of its parameter space is excluded by the null results of the KARMEN data .

One of the most stringent constraints on the neutrino mass comes from the cosmological consideration to avoid over-closing the universe. In the standard cosmology , the stable neutrinos should be no heavier than $`20`$ eV. Recently it has been shown that if the reheating temperature $`T_{RH}`$ is low, the densities of neutrinos can become much less than usually assumed . The cosmological constraint is very much relaxed. Indeed, for $`1`$ MeV $`<T_{RH}<3`$ MeV, the abundance of tau- and muon-neutrinos is

$$\mathrm{\Omega }_{\nu _\tau }h^2=\mathrm{\Omega }_{\nu _\mu }h^2=\left(\frac{m_\nu }{4\mathrm{keV}}\right)\left(\frac{T_{RH}}{1\mathrm{MeV}}\right)^3.$$
(1)
This leads to Ref.  taking $`\nu _\mu `$ and $`\nu _\tau `$ as the WDM.

The above WDM consideration makes the neutrino mass pattern quite unique. First the Super-K results constrain the mass of both $`\nu _\mu `$ and $`\nu _\tau `$ around keV scale . Their mass-squared difference of order $`10^3`$ eV2 implies that the masses of $`\nu _\mu `$ and $`\nu _\tau `$ are highly degenerate. In this case, the solar neutrino problem can only be understood by introducing a light sterile neutrino, $`\nu _s`$. The $`\nu _e`$ (and in the large mixing case the $`\nu _s`$) should be lighter than eV order, for the consistence with the laboratory experiments of $`m_{\nu _e}`$. The LSND result, however, cannot be compatible with the presence of keV scale muon-neutrinos. Instead the mixing between the $`\nu _\mu `$ and $`\nu _e`$ should be very small, as expected from the keV scale mass of $`\nu _\mu `$ and the sub-eV scale mass of $`\nu _e`$.

In this paper, a four light neutrino model is built to give the above-mentioned neutrino mass pattern. How to construct the mass spectrum of four light neutrinos with hierarchies is a theoretically challenging problem. We extend a method proposed in Ref. . It simply incorporates the seesaw mechanism  with flavor symmetry . Introducing three right-handed neutrinos and assuming singularity in both the Dirac and Majorana mass matrices, the neutrino mass matrix has approximately the following form:

$$=\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& m_{22}& m_{23}\\ 0& 0& 0& 0& m_{32}& m_{33}\\ 0& 0& 0& 0& 0& 0\\ 0& m_{22}& m_{32}& 0& M_{22}& M_{23}\\ 0& m_{23}& m_{33}& 0& M_{23}& M_{33}\end{array}\right).$$
(2)
In this mass spectrum, there are two heavy neutrinos of masses $`M`$, two light neutrinos of masses $`m^2/M`$, and two massless neutrinos; three mass scales, $`M`$, $`m^2/M`$, and $`0`$, characterizes this model. Thus the four light neutrinos are naturally divided into two pairs with a mass gap of $`m^2/M`$, and each pair consists of two degenerate mass eigenstates.

Abelian flavor symmetry  can be used to generate such a neutrino mass matrix. Furthermore this symmetry should also provide the maximal mixing for the atmosphere neutrino anomaly which is not guaranteed in the form of $``$ in Eq. (2). A softly breaking of the symmetry is necessary to generate small masses for two massless neutrinos and to lift the degeneracy in each pair. In the following discussion, supersymmetry is implied. The flavor symmetry is spontaneously broken by a vacuum expectation value (VEV) of an electroweak singlet field $`X`$. As long as the flavor charges balance under the symmetry, the following interactions are allowed:

$$L_\alpha HN_\beta \left(\frac{X_{\mathrm{VEV}}}{\mathrm{\Lambda }}\right)^{m_{\alpha \beta }},MN_\alpha N_\beta \left(\frac{X_{\mathrm{VEV}}}{\mathrm{\Lambda }}\right)^{n_{\alpha \beta }},$$
(3)
where $`L_\alpha `$ ($`\alpha =e`$, $`\mu `$, $`\tau `$), $`H`$ and $`N_\alpha `$ denote the lepton doublets, a Higgs field, and the right-handed neutrino fields, respectively. The $`\mathrm{\Lambda }`$ is the flavor symmetry breaking scale, and the $`m_{\alpha \beta }`$ and $`n_{\alpha \beta }`$ are non-negative integers, required for the holomorphy of the superpotential. The order parameter for this new symmetry is defined by

$$\lambda \frac{X_{\mathrm{VEV}}}{\mathrm{\Lambda }}1,$$
(4)
which is a typical order of Cabbibo angle $`0.1`$

The assignment of the Abelian flavor charges relevant to Eq. (2) and to maximal $`\nu _\mu \nu _\tau `$ mixing is assumed to be

$`L_e(z+1/2),L_\mu (a),L_\tau (a),`$ (5)
$`E_e^c(z+9/2),E_\mu ^c(a+3),E_\tau ^c(a+2),`$ (6)
$`N_e(x1/2),N_\mu (a),N_\tau (a),X(1),`$ (7)
with the positive integers $`a`$ and $`x(>a>1)`$. The sign of the integer $`z`$ is to be set later. The $`E_\alpha ^c`$’s in Eq. (5) are the anti-particle fields of the SU(2) singlet charged leptons. In order to obtain the physical mixing angles of neutrinos, the charged lepton mass matrix should be simultaneously taken into account. The gauge and Higgs bosons possess vanishing flavor charges. Note that the flavor charges for the first generation are half-integers while those for the second and third generations are integers expressed by a single parameter $`a`$. Compared to analogous analysis for three light neutrino scenario, the choice of the flavor charges here is more limited. One tricky point is that one of the right-handed neutrino masses is made to be vanishingly small ($`m^2/M`$).

The flavor charge assignment in Eq. (5) produces the Dirac and Majorana mass matrices of neutrinos as

$$_D=m\left(\begin{array}{ccc}Y_{11}\lambda ^{x+z}& 0& 0\\ 0& 1& Y_{23}\lambda ^{2a}\\ 0& 0& 1\end{array}\right),_M=M\left(\begin{array}{ccc}\zeta _1\lambda ^{2x1}& 0& 0\\ 0& 0& 1\\ 0& 1& \zeta _4\lambda ^{2a}\end{array}\right),$$
(8)
and the mass matrix of charged leptons as

$$_l=m\left(\begin{array}{ccc}\eta _{11}\lambda ^5& 0& 0\\ 0& \eta _{22}\lambda ^3& \eta _{23}\lambda ^{2a+2}\\ 0& 0& \eta _{33}\lambda ^2\end{array}\right),$$
(9)
where $`Y`$’s, $`\zeta `$’s and $`\eta `$’s are order one coefficients. The $`_l`$ and $`_D`$ are almost diagonal, while the $`_M`$ is mainly off-diagonal. Therefore to leading order, the mass matrix of four light neutrinos in the $`(\nu _s,\nu _e,\nu _\mu ,\nu _\tau )`$ basis is obtained as,

$$_\nu ^{(0)}\frac{m^2}{M}\left(\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right).$$
(10)
The neutrino mass spectrum due to $`_\nu ^{(0)}`$ is

$$m_{\nu _s}=m_{\nu _e}=0,m_{\nu _\mu }=m_{\nu _\tau }=\frac{m^2}{M},\mathrm{sin}\theta _{\mu \tau }=\frac{1}{\sqrt{2}}.$$
(11)
It is to be noted that the maximal mixing between the $`\nu _\mu `$ and $`\nu _\tau `$ results from the flavor symmetry. The requirement of keV scale $`\nu _\mu `$ and $`\nu _\tau `$, i.e.,

$$\frac{m^2}{M}1\mathrm{keV},$$
(12)
is achieved via the seesaw mechanism with $`m300`$ GeV and $`M10^{11}`$ GeV. The $`ϵ`$, defined by the ratio $`m/M`$, is then of order $`10^{8.5}`$.

First let us examine the mass spectrum in the charged lepton sector. The eigenvalues of the $`_l`$ are of order $`\lambda ^5m`$, $`\lambda ^3m`$, and $`\lambda ^2m`$, which yields appropriate mass scales of the charged leptons. The $`_l`$ is diagonalized by

$$R_l^L_lR_l^R=\mathrm{Diag}(m_e,m_\mu ,m_\tau ),$$
(13)
where the $`R_l^L`$ diagonalizes the hermitian mass-squared matrix $`_l_l^{}`$.

In the neutrino sector the mass matrix of four light neutrinos can be obtained by the method described in Ref. . We finally have the following symmetric mass matrix of four light neutrinos in the $`(\nu _s,\nu _e,\nu _\mu ,\nu _\tau )`$ basis with $`\lambda ^x/ϵ\lambda ^\beta `$:

$$_\nu =\frac{m^2}{M}\left(\begin{array}{cccc}\lambda ^{2\beta 1}& \lambda ^{\beta +z}& 0& 0\\ \lambda ^{\beta +z}& 0& 0& 0\\ 0& 0& \lambda ^{2a}& 1\\ 0& 0& 1& 0\end{array}\right),$$
(14)
where the charged lepton mixing effects are incorporated so that the charged lepton fields have been rotated into mass eigenstates. In Eq. (14) the combinations of order one parameters such as $`Y_{ij}`$, $`\eta _i`$, and $`\zeta _i`$ are omitted, for simplicity.

In our scenario, there is no mixing between the two pairs due to our assignment of half-integer charges for the first generation, and integer charges for the second and third generations. The absence of the $`\nu _e\nu _\mu `$ mixing is compatible with the KARMEN data, predicts null results in future laboratory searching for the $`\nu _e\nu _\mu `$ or $`\nu _e\nu _\tau `$ mixing, and has no influence on astrophysical processes . Then the $`_\nu `$ is automatically factorized out into $`M_{\nu _s\nu _e}`$ and $`M_{\nu _\mu \nu _\tau }`$.

The requirement of Super-K data determines the value of $`a`$ as

$$\mathrm{\Delta }^2m_{23}\lambda ^{2a}\left(\frac{m^2}{M}\right)^23\times 10^3\mathrm{eV}^2,\mathrm{for}a=4.$$
(15)
The effective mass matrix of $`\nu _e`$ and $`\nu _s`$ in the unit of eV can be written by

$$M_{\nu _s\nu _e}\mathrm{eV}\left(\begin{array}{cc}\lambda ^{2\beta 4}& \lambda ^{\beta +z3}\\ \lambda ^{\beta +z3}& 0\end{array}\right).$$
(16)
The case of $`\beta 1<z`$ $`(\beta 1>z)`$ corresponds to the small (large) mixing angle solution for the solar neutrino problem. Some possible solutions are listed in the following table. Note that in the SMA case, the $`m_{\nu _e}`$ is the smaller mass since the $`M_{\nu _s\nu _e}`$ in Eq. (16) is in the $`(\nu _s,\nu _e)`$ basis.

| Solution | $`z`$ | $`\beta (x)`$ | $`m_{\nu _e}`$ (eV) | $`\mathrm{\Delta }m_{}^2`$ (eV2) | $`\mathrm{sin}^22\theta _{}`$ |
| --- | --- | --- | --- | --- | --- |
| SMA | 4 | 3.5 (12) | $`10^6`$ | $`10^6`$ | $`10^3`$ |
| LMA | 1 | 3.5 (12) | $`3\times 10^2`$ | $`3\times 10^5`$ | $`1`$ |
| LOW-QVO | 1 | 4.5 (13) | $`3\times 10^3`$ | $`3\times 10^8`$ | $`1`$ |
|  | 2 | 4.5 (13) | $`3\times 10^4`$ | $`3\times 10^9`$ | $`1`$ |
|  | 3 | 4.5 (13) | $`3\times 10^5`$ | $`3\times 10^{10}`$ | $`1`$ |
|  | 0 | 5.5 (14) | $`3\times 10^3`$ | $`3\times 10^{10}`$ | $`1`$ |
|  | $`1`$ | 5.5 (14) | $`3\times 10^2`$ | $`3\times 10^9`$ | $`1`$ |
|  | $`2`$ | 5.5 (14) | 0.3 | $`3\times 10^8`$ | $`1`$ |

It is to be checked in the SNO  and KamLAND experiments. For one of the LOW-QVO solution, the mass of the electron neutrino is just an order of magnitude lower than the current experiment limit. It is possible to probe this value in future experiments directly. In addition, it may have observable effect in the cosmic microwave background anisotropies. The $`\nu _e`$ and $`\nu _s`$ compose a Dirac neutrino which does not result in any observable neutrinoless double $`\beta `$ decay of relevant next generation experiments.

In summary, it would be simple if the muon- and tau-neutrinos are just the WDM. We have presented a neutrino model which generates keV scale $`\nu _\mu `$ and $`\nu _\tau `$ and meanwhile provides the neutrino oscillation solutions for the solar and atmospheric neutrino experiments. It combines the seesaw mechanism and the Frogatt-Nielsen mechanism. The key point is that the neutrino mass matrix has a singular form. A light sterile neutrino is obtained due to the flavor symmetry.

Finally some remarks should be mentioned.

* Compared to the model in Ref. , the light neutrino spectrum is similar to the ordinary $`2+2`$ light neutrino scheme . But the two neutrino pairs are more widely separated (by keV). For the solar neutrino problem, the SMA, LMA or LOW-QVO oscillation solutions into light sterile neutrinos can be accommodated in our scenario. This is achieved by taking the U(1) flavor charges to be half-integers for the first generation and to be integers for the second and third generations. Of course, the LSND results cannot be relevant in this model. But our results are compatible with the KARMEN data.
* The LMA and LOW-QVO solutions have been obtained in fact due to the singular seesaw mechanism. The singular seesaw mechanism was applied to the atmospheric neutrino anomaly in Ref. . However, the current Super-K data do not favor the $`\nu _\mu \nu _s`$ oscillation. What we have done in Ref.  and in this paper for LMA and LOW-QVO scenarios are to obtain a singular seesaw mechanism for the solar neutrinos.
* It seems a drawback to introduce both the seesaw and the Frogatt-Nielsen mechanism to make neutrinos light. However the seesaw mechanism itself cannot predict the neutrino flavor structure, especially the maximal $`\nu _\mu \nu _\tau `$ mixing. Another underlying motivation is to make the four light neutrino scenario easier to be understood in the framework of grand unification theories, which will be studied further.
* Other interesting aspects of keV neutrinos were discussed before. Such neutrinos may be responsible for the large velocity of pulsars . Because the neutrino densities in the early universe are much smaller than that assumed in the standard cosmology, the cosmological constraints for neutrinos from the big bang nucleosynthesis is different from that in the standard one . This needs more detailed study. In addition to the keV active neutrinos, certain amount of some other keV sterile neutrinos  may also contribute to the WDM.

###### Acknowledgements.

C.L. was supported in part by the National Natural Science Foundation of China with grant no. 10047005.