text:
# Excitation lines and the breakdown of Stokes–Einstein relations in supercooled liquids

## I Introduction

Normal liquids exhibit homogeneous behavior in their dynamical properties over length scales larger than the correlation length of density fluctuations. For example, the Stokes–Einstein relation that relates the self–diffusion constant $`D`$, viscosity $`\eta `$, and temperature $`T`$,

$$D\frac{T}{\eta },$$
(1)
is usually accurate Hansen and McDonald (1994); Balucani and Zoppi (1994). This relation is essentially a mean field result for the effects of a viscous environment on a tagged particle. In recent experimental studies, it has been reported that the Stokes–Einstein relation breaks down as the glass transition is approached in supercooled liquid systems Ediger (2000); Fujara et al. (1992); Chang and Sillescu (1997); Cicerone and Ediger (1996); Blackburn et al. (1996); Swallen et al. (2003). Translational diffusion shows an enhancement by orders of magnitude from what would be expected from Eq. (1) Hodgdon and Stillinger (1993); Stillinger and Hodgdon (1994); Tarjus and Kivelson (1995); Liu and Oppenheim (1996); Xia and Wolynes (2001). Here, we show that this breakdown is due to fluctuation dominance in the dynamics of low temperature glass formers. These pertinent fluctuations are dynamic heterogeneities Perera and Harrowell (1996, 1999, 1998); Kob et al. (1997); Donati et al. (1998, 1999a, 1999b); Glotzer (2000). Thus, the Stokes–Einstein breakdown is one further example of the intrinsic role of dynamic heterogeneity in structural glass formers Palmer et al. (1984); Garrahan and Chandler (2002, 2003).

In the treatment we apply, dynamic heterogeneity is a manifestation of excitation lines in space–time Garrahan and Chandler (2002). This picture leads to the prediction of dynamic scaling in supercooled liquids, $`\tau (l)l^z`$. Here, $`\tau (l)`$ is the structural relaxation time for processes occurring at length scale $`l`$, and $`z`$ is a dynamic exponent for which specific results have been established Garrahan and Chandler (2002, 2003); Whitelam et al. (2003). This picture and its predicted scaling results differ markedly from those derived with the view that glass formation is a static or thermodynamic phenomenon Debenedetti and Stillinger (2001); Sastry et al. (1998); Büchner and Heuer (1999, 2000); Schroder et al. (2000); Keyes and Chowdhary (2002); Kivelson et al. (1995); Xia and Wolynes (2000). It also differs from mode coupling theory which predicts singular behavior at non–zero temperature Götze and Sjögren (1992); Kob (2003).

This paper is organized as follows. In Sec. II, we introduce our model for a supercooled liquid with a probe molecule immersed in the liquid. Simulation results are given in Secs. III and IV. Section IV also provides analytical analysis of the diffusion coefficient and the Stokes–Einstein violation, and explains the origin of the decoupling of transport coefficients based on the excitation line picture of trajectory space. Comparison of our theory with recent experimental results is carried out in Sec. V. We conclude in Sec. VI with a Discussion.

## II Models

We imagine coarse graining a real molecular liquid over a microscopic time scale (e.g., larger than the molecular vibrational time scale), and also over a microscopic length scale (e.g., larger than the equilibrium correlation length). In its simplest form, we assume this coarse graining leads to a kinetically constrained model Garrahan and Chandler (2002, 2003); Fredrickson and Andersen (1984); Jäckle and Eisinger (1991); Ritort and Sollich (2003) with the dimensionless Hamiltonian,

$$H=\underset{i=1}{\overset{N}{}}n_i,(n_i=0,1).$$
(2)
Here, $`n_i=1`$ coincides with lattice site $`i`$ being a spatially unjammed region, while $`n_i=0`$ coincides with it being a jammed region. We call $`n_i`$ the “mobility field”. The number of sites, $`N`$, specifies the size of the system. From Eq. (2), thermodynamics is trivial, and the equilibrium concentration of defects or excitations is

$$c=n_i=\frac{1}{1+\mathrm{exp}(1/\stackrel{~}{T})},$$
(3)
where $`\stackrel{~}{T}`$ is a reduced temperature. We make explicit connection of $`\stackrel{~}{T}`$ with absolute temperature later when comparing our theory with experimental results.

The dynamics of these models obey detailed balance and local dynamical rules. Namely,

$$n_i=0\begin{array}{c}\stackrel{k_i^{(+)}}{}\\ \underset{k_i^{()}}{\overset{}{}}\end{array}n_i=1,$$
(4)
where the rate constants for site $`i`$, $`k_i^{(+)}`$ and $`k_i^{()}`$, depend on the configurations of nearest neighbors. For example, in dimension $`d=1`$,

$`k_i^{(+)}`$ $`=`$ $`e^{1/\stackrel{~}{T}}f(n_{i1},n_{i+1}),`$ (5)
$`k_i^{()}`$ $`=`$ $`f(n_{i1},n_{i+1}),`$ (6)
where $`f(n_{i1},n_{i+1})`$ reflects the type of dynamical facilitation. In the Fredrickson–Andersen (FA) model Fredrickson and Andersen (1984), a state change is allowed when it is next to at least one defect. The facilitation function in this case is given by,

$$f_{\mathrm{FA}}(n_{i1},n_{i+1})=n_{i1}+n_{i+1}n_{i1}n_{i+1}.$$
(7)
In the East model Jäckle and Eisinger (1991), dynamical facilitation has directional persistence. The facilitation function in this case is

$$f_{\mathrm{East}}(n_{i1},n_{i+1})=n_{i1}.$$
(8)
In order to study translational diffusion in supercooled liquids, we extend the concept of dynamic facilitation to include a probe molecule. The dynamics of a probe will depend on the local state of the background liquid. When and where there is no mobility, the diffusive motion of the probe will be hindered. When and where there is mobility, the probe molecule will undergo diffusion easily. As such, in a coarse–grained picture, the probe molecule is allowed to jump from lattice site $`i`$ to a nearest neighbor site when site $`i`$ coincides with a mobile region, $`n_i=1`$. In order to satisfy detailed balance, we further assume that the probe molecule can move only to a mobile region, i.e.,

$$x(t+\delta t)=x(t)\pm \delta xn_xn_{x\pm \delta x},$$
(9)
where $`x(t)`$ denotes the position of the probe at time $`t`$. Units of time and length scales are set equal to a Monte Carlo sweep and a lattice spacing, respectively.

## III Computer Simulations

Using the rules described in Sec. II, we have performed Monte Carlo simulations of diffusion of a probe molecule in the FA and East models for various temperatures. For the purpose of numerical efficiency, we have used the continuous time Monte Carlo algorithm Bortz et al. (1975); Newman and Barkema (1999). In the all systems, $`N`$ was chosen as $`N=100/c`$, and the simulations were performed for total times $`𝒯100\tau `$, with $`\tau `$ being the relaxation time of the model. Averages were performed over $`10^3`$ to $`10^5`$ independent trajectories.

In Fig. 1, we show typical trajectories of probe molecules in the FA and East models. In the high temperature case, trajectory space is dense with mobile regions and there are no significant patterns in space–time. As such, the dynamics is mean–field like. It is for this reason that the relaxation time in this case is inversely proportional to the equilibrium probability of excitation, $`c`$ (see, for example, Ref. Berthier and Garrahan (2003a)). The probe molecule executes diffusive motion, without being trapped in immobile regions for any significant period of time.

The low temperature dynamics is different. Mobility is sparse, defects tend to be spatially isolated at a given time, and trajectory space exhibits space–time patterns. See Fig. 1(b) and (c). Because of the facilitation constraint, an immobile region needs a nearest mobile region to become mobile at a later time. The excitations therefore form continuous lines and bubble–like structures in trajectory space. While inside a bubble, the probe molecule will be immobilized. See, for example, the segment of the trajectory of a probe molecule for $`0<t<500`$ in Fig. 1(b). Due to exchanges between mobile and immobile regions, an immobile region can become mobile after a period of time. At that stage the probe molecule can perform a random walk until it is again in an immobile region. The motion of a probe molecule will manifest diffusive behavior over a time long enough for many dynamical exchanges to occur. In the East model at low temperatures such as pictured in Fig. 1(c), the bubbles in space–time form hierarchical structures Garrahan and Chandler (2002).

Figure 2 plots mean square displacements of probe molecules for the FA and East models for three different cases pictured in Fig. 1. In the high temperature case, the mean square displacement reaches its diffusive linear regime after a very short transient time. In the low temperature case, the probe molecule in the East model case reaches the diffusive regime after a longer time and over a larger length scale than that in the FA model with the same reduced temperature.

## IV Stokes–Einstein Violation

### IV.1 Diffusion Coefficient

Figure 3 plots the diffusion coefficient of a probe molecule for the FA and the East models. The diffusion coefficient is determined from the mean square displacement,

$$D=\underset{t\mathrm{}}{lim}\frac{[\mathrm{\Delta }x(t)]^2}{t},$$
(10)
where $`\mathrm{\Delta }x(t)=x(t)x(0)`$. Error estimates for our simulations are no larger than the size of the symbols.

In the FA model, the diffusion coefficient exhibits Arrhenius behavior for $`\stackrel{~}{T}<1`$. This behavior reflects the fact that relaxation dynamics in the FA model is similar to that of a strong liquid. In this regime, over more than 4 orders of magnitude in $`D`$, the slope of $`\mathrm{log}D`$ vs $`\stackrel{~}{T}^1`$ is close to 2. This result is consistent with the expected low temperature scaling,

$$D_{\mathrm{FA}}c^2\mathrm{exp}(2/\stackrel{~}{T}),$$
(11)
as discussed in the next subsection. In the East model case, also pictured in Fig. 3, the diffusion coefficient decreases more quickly than Arrhenius. This super–Arrhenius behavior is due to the hierarchical nature of dynamics in the East model Sollich and Evans (1999).

Comparing the diffusion coefficients with the relaxation times of the background liquids demonstrates Stokes–Einstein violation in both models. The relaxation times, $`\tau `$, of the FA and the East models at different temperatures have been determined in prior work Buhot and Garrahan (2001); Berthier and Garrahan (2003b). When the Stokes–Einstein relation is satisfied, $`D\tau \mathrm{const}`$. This behavior occurs in the FA and East models when $`\stackrel{~}{T}>2`$, but Fig. 4 shows that $`D\tau `$ is enhanced from that behavior by 2 or 3 orders of magnitude when $`\stackrel{~}{T}<1`$. Bear in mind, these deviations from Stokes–Einstein are $`d=1`$ results. The appropriate generalization of the FA model to $`d=3`$ does not exhibit such large deviations. On the other hand, we expect that generalizations of the East model, which is hierarchical and therefore fragile, will have weak dimensional dependence and continue to exhibit large deviations for $`d=3`$. We turn to the arguments that explain these claims now.

### IV.2 Scaling Analysis

For high temperatures, the local mobility field will tend to be close to its mean value, $`c`$. As such, both the relaxation mechanism of the material and the diffusional motion of the probe molecule make use of the same local mobility fields. For this reason, the diffusion coefficient and the relaxation time scale are strongly coupled in this regime, leading to the Stokes–Einstein relation.

At low temperatures, however, the dynamics of the system is not so simply related to the mean mobility field. Here, the fluctuations of bubble structures dominate. The relaxation time of the background liquid will approximately scale as the longest temporal extension of bubbles. The persistence time of an individual lattice site, $`t_{\mathrm{pers}}`$, is the time for which that site makes its first change in state. Its typical size will be intimately tied to the structural relaxation time of the liquid. For the FA model in $`d=1`$,

$$\tau t_{\mathrm{pers}}c^3.$$
(12)
See, for example, Refs. Ritort and Sollich (2003); Garrahan and Chandler (2002).

This result is consistent with a simple argument concerning diffusive motions of excitation lines in the low temperature FA model Garrahan and Chandler (2002). In particular, the structural relaxation times in the FA model is given by the time in which a typical bubble structure looses its identity through wandering motions of excitation lines. The excitation line has a local diffusivity of $`𝒟c`$. (We use caligraphic $`𝒟`$ to distinguish this diffusion constant for excitations from that for particles, $`D`$.) In order to form a bubble, an excitation line needs to wander distance of the order of the typical length between defects, $`l_{\mathrm{eq}}c^1`$. Therefore, the mean relaxation time is given by $`\tau l_{\mathrm{eq}}^2/𝒟c^3`$.

When the probe molecule is at the boundary of a bubble, it may not need to wait until the bubble closes in order to undergo diffusion; rather, it can remain within mobile cells and diffuse around the boundary of the bubbles. In this way, translational diffusion will be more facilitated than structural relaxation, leading to an enhanced diffusion in the fluctuation dominated low temperature region. Specifically, consider the dynamical exchange times, i.e., the times between flipping events for a given lattice site. See Fig. 5. $`t_0`$ is such a time duration for an $`n_i=0`$ state and $`t_1`$ is such a time duration for an $`n_i=1`$ state. The probe molecule can move only while in a mobile region. Further, the mean square displacement of the probe will be proportional to the number of diffusive steps that a probe molecule will take during the trajectory, $`𝒩`$,

$$[\mathrm{\Delta }x(t)]^2𝒩\frac{𝒯}{t_0+t_1}.$$
(13)
Here, $`𝒯`$ is the length of a long trajectory in the FA model. The average duration of the defect state, $`t_1`$, is inversely proportional to the probability of a lattice site being mobile, $`c`$, times the flip rate, $`k_i^{()}`$. Since $`k_i^{()}𝒪(1)`$, we have,

$$t_1c^1.$$
(14)
From detailed balance, therefore,

$$t_0c^2.$$
(15)
Since $`t_1t_0`$ in the low temperature region, Eqs. (13)–(15) give

$$D_{\mathrm{FA}}\frac{(\mathrm{\Delta }x)^2}{𝒯}\frac{1}{t_0}c^2.$$
(16)
This result explains Eq. (11). Together with Eq. (12), it leads to

$$D_{\mathrm{FA}}\tau ^\xi ,$$
(17)
with $`\xi =2/3`$ in the $`d=1`$ FA model case. This scaling is to be contrasted with the Stokes–Einstein result, $`\xi =1`$.

Numerical simulation Garrahan and Chandler (2003) and renormalization group analysis Whitelam et al. (2003) of higher dimension generalizations of the FA model indicate that for $`d=3`$, $`\tau c^{2.1}`$. However, the scaling $`Dc^2`$ remains true for all dimensions as it is based solely on detailed balance. Thus, for $`d=3`$, $`\xi 0.95`$. In other words, there is only a weak breakdown in the Stokes–Einstein relation for strong liquids in $`d=3`$.

In the East model case, both the diffusion coefficient and the relaxation time show super-Arrhenius behavior. The hierarchical, fractal structure of pattern development in trajectory space for the East model does not allow a simple scaling analysis of the diffusion coefficient, and it is not obvious whether temperature independent scaling exists. One can define a temperature dependent scaling exponents, $`\alpha (\stackrel{~}{T})`$ and $`z(\stackrel{~}{T})`$,

$`D_{\mathrm{East}}`$ $``$ $`c^{\alpha (\stackrel{~}{T})},`$ (18)
$`\tau `$ $``$ $`l^{z(\stackrel{~}{T})},`$ (19)
so that

$$D_{\mathrm{East}}\tau ^{\alpha (\stackrel{~}{T})/z(\stackrel{~}{T})}.$$
(20)
Interestingly, our numerical results indicate that $`\xi =\alpha /z0.73`$ is independent of temperature as shown in Fig. 6. This exponent, $`\xi 0.73`$ for the $`d=1`$ East model, is very close to what many experiments and simulations have found for three–dimensional glass forming liquids. For example, a recent experiment finds that $`\xi 0.77`$ in the self–diffusion of tris-naphthylbenzene(TNB) Swallen et al. (2003). It was found that $`\xi 0.75`$ in a molecular dynamics simulation of Lennard–Jones binary mixture Yamamoto and Onuki (1998) and a recent detailed scaling analysis of numerical results shows $`\xi 0.65`$ Berthier (2003).

Presumably, such good agreement of scaling relation between the $`d=1`$ East model and higher dimension systems arises due to directional persistence of facilitation in the fragile liquid Garrahan and Chandler (2002, 2003). This persistence in higher dimensions causes motion to be effectively one–dimensional Garrahan and Chandler (2003). Therefore, dimensionality is not very significant for fragile glass formers. As such, for fragile systems, we expect that the scaling relation of the Stokes–Einstein violation will be reasonably well described by the $`d=1`$ East model. Based on this expectation, we further purse the comparison between theory and experiment.

## V Comparison with Experiment

Swallen et al. Swallen et al. (2003) measured the self–translational diffusion coefficient of TNB near the glass transition temperature. They observed an increase of $`D\eta /T`$ from its high temperature limit by a factor of 400 near the glass transition temperature. In order to compare our results with these experiments, we need to determine the excitation concentration, $`c`$, as function of temperature. Since TNB behaves as a fragile liquid, we determine the excitation concentration as a function of temperature by fitting the viscosity data of TNB Plazek and Magill (1968) with the generalization of the East model formulas to higher dimensions Garrahan and Chandler (2003). Namely,

$$\mathrm{ln}\tau \frac{1}{d\mathrm{ln}2}\left[\mathrm{ln}(g/c)\right]^2,$$
(21)
where $`g`$ is the number of equally likely persistence directions on a cubic lattice, and

$$\mathrm{ln}(c)=\mathrm{ln}(c_R)J\left(\frac{1}{T}\frac{1}{T_R}\right).$$
(22)
The parameter $`J`$ is the energy scale associated with creating a mobile region from an immobile region, and $`T_R`$ is an appropriate reference temperature. Details on the fitting can be found in Ref. Garrahan and Chandler (2003). Taking $`g=8`$ (the cubic lattice value) and $`T_R`$ as the temperature at which $`\mathrm{log}\tau `$ is half the value of $`\mathrm{log}\tau (T_g)`$, we determine that $`J/T_g21.7`$, and $`\mathrm{log}_{10}(c_R/g)1.28`$. The reduced temperature, $`\stackrel{~}{T}`$ of the East model is related to absolute temperature by

$$\frac{1}{\stackrel{~}{T}}=J\left(\frac{1}{T}\frac{1}{T_R}\right)+\mathrm{ln}(g/c_R).$$
(23)
Once we have determined the excitation concentration as a function of the temperature, we can compare experimental data with our computed results for the Stokes–Einstein violation in the East model case. Based on the argument that the scaling relation of the Stokes–Einstein violation ($`D\tau ^\xi `$) remains robust in higher dimensions and from the dimensional dependence of Eq. (21), we expect

$$\mathrm{ln}(D\tau )|_{d=3}\frac{1}{3}\mathrm{ln}(D\tau )|_{d=1}.$$
(24)
In Fig. 7 we use this relationship to compare the extent of the Stokes–Einstein violation of the experimental system with our East model results. The agreement between the two appears excellent.

## VI Discussion

There have been previous theoretical studies on the violation of the Stokes–Einstein relation in supercooled liquid systems. For example, Kivelson and Tarjus have argued that the Stokes–Einstein violation can be understood from their “frustration–limited domain” model for supercooled liquids Tarjus and Kivelson (1995); Kivelson et al. (1995). Assuming a distribution of local relaxation times associated with domain structures, this model describes the translational diffusion and viscosity as corresponding to different averaging process of such a distribution. Their idea contrasts to ours in that the domain structure in their work is purely static, and the exchange between different domains are not considered.

Hodgdon and Stillinger have proposed a fluidized domain model Hodgdon and Stillinger (1993); Stillinger and Hodgdon (1994). In their work, it is assumed that the system consists of a sparse collection of fluid–like domains in a background of more viscous media, and fluid–like domains appear and disappear with a finite life–time and rate. Relaxation times are determined by the rate of appearance of the fluid–like domains, while translation diffusion also depends on the life–time of the domains. To the extent that these domains refer to space–time and not simply space, this picture is not inconsistent with ours. Xia and Wolynes have applied the so–called “random first order transition theory” Xia and Wolynes (2000) to the Hodgdon–Stillinger model Xia and Wolynes (2001). In this case, the picture is both mean field and static and decidedly contrary to our fluctuation dominated and dynamic view.

From the perspective that Stokes–Einstein violation is a manifestation of fluctuation dominated dynamics, one expects that similar decoupling behavior occurs between other kinds of transport properties near the glass transition. The extent to which such decoupling can appear depends upon microscopic details in the specific transport properties and materials under study. For example, molecular rotations of a probe will be coupled to the mobility field, but less so than translations. Indeed, single molecule experiments indicate that rotations persist in both mobile and immobile regions of a glass former Deschenes and Bout (2001, 2002); Bartko et al. (2002). Rotational motions can therefore average the effects of dynamic heterogeneity to a greater extent than translational motions. As such, decoupling of rotational relaxation from structural relaxation can be more difficult to detect than violations of the Stokes–Einstein relation. Precisely how such effects might be detected seems worthy of further theoretical analysis.

###### Acknowledgements.

We are grateful to M.D. Ediger, D.A. VandenBout, L. Berthier and S. Whitelam for discussions. This work was supported at Berkeley by the Miller Research Fellowship (YJ) and by the US Department of Energy Grant No. DE-FG03-87ER13793 (DC), at Oxford by EPSRC Grant No. GR/R83712/01 and the Glasstone Fund (JPG), and at Nottingham by EPSRC Grant No. GR/S54074/01 (JPG).


format: text

text:
# 1 The exclusion plot in the parameter space for millicharged particles. Light grey area is excluded by accelerator experiments and BBN. Dashed region is excluded by the relic abundance of millicharged particles in models without paraphoton. Part of this region above the dash-dotted line is excluded by the relic abundance in models with paraphoton (assuming 𝛼'=0.1). On the left of the dotted line millicharged particles cannot be thermally produced in the Early Universe. Dark grey and dashed dark grey areas are the previously allowed regions which are now excluded by Eq. () in models without and with paraphoton, respectively. On the right of the dark solid line millicharged particles are coupled to baryons (see Eq. () in the text).

Narrowing the window for millicharged particles by CMB anisotropy
S. L. Dubovskya111e-mail: sergd@ms2.inr.ac.ru, D. S. Gorbunova222e-mail: gorby@ms2.inr.ac.ru, G. I. Rubtsova,b333e-mail: grisha@ms2.inr.ac.ru

aInstitute for Nuclear Research of the Russian Academy of Sciences,
60th October Anniversary prospect 7a, Moscow 117312, Russia
bMoscow State University, Department of Physics,
Vorobjevy Gory, Moscow, 119899, Russia

## Abstract

We calculate the cosmic microwave background (CMB) anisotropy spectrum in models with millicharged particles of electric charge $`q10^610^1`$ in units of electron charge. We find that a large region of the parameter space for the millicharged particles exists where their effect on the CMB spectrum is similar to the effect of baryons. Using WMAP data on the CMB anisotropy and assuming Big Bang nucleosynthesis value for the baryon abundance we find that only a small fraction of cold dark matter, $`\mathrm{\Omega }_{mcp}h_0^2<0.007`$ (at 95% CL), may consists of millicharged particles with the parameters (charge and mass) from this region. This bound significantly narrows the allowed range of the parameters of millicharged particles. In models without paraphoton millicharged particles are now excluded as a dark matter candidate. We also speculate that recent observation of 511 keV $`\gamma `$-rays from the Galactic bulge may be an indication that a (small) fraction of CDM is comprised of the millicharged particles.

Search for particles carrying small but non-vanishing electric charge (millicharged particles) has long history. If observed, millicharged particles would either cause serious doubts on the concept of Grand Unification or imply the existence of a new massless gauge boson — paraphoton . Furthermore, the existence of millicharged particles would hint towards processes with apparent electric charge non-conservation, like electron or proton decay to “nothing” .

There are various constraints on the parameters (charge and mass) of millicharged particles, coming from collider and laboratory experiments and from cosmology and astrophysics (see, e.g., Refs.  for the latest results and Refs.  for reviews), see Fig. 1. Interestingly, reported bounds did not exclude a possibility  that a significant part (or even all) of the cold dark matter (CDM) is comprised of millicharged particles.

The constraints on the parameters of millicharged particles are somewhat different in theories with and without paraphoton. Without paraphoton, two domains in the parameter space of millicharged particles are allowed. The first one corresponds to heavy particles with tiny electric charge (left upper corner in Fig. 1) which would never be produced thermally in the early Universe. This region is far beyond the reach of future collider and laboratory experiments. In the current Letter we are concerned with another region. This is a narrow window of relatively light particles with masses $`M10^3÷10^2`$ GeV and charges $`q10^6÷10^1`$. Larger charges $`q0.1`$ are ruled out by limits on the cosmic ray fluxes of fractionally charged particles . It is worth noting that particles with $`q0.2`$ are also excluded by measurements of the width of Z-boson (cf. Ref. ) if one makes use of the latest data .

In model with paraphoton, and for not very small values of the paraphoton coupling constant $`\alpha ^{}`$, millicharged particles annihilate mainly into pairs of paraphotons. As a result, their annihilation in the early Universe is more efficient and the cosmological bound coming from the relic abundance depends on the value of $`\alpha ^{}`$ and generically is less restrictive (see Fig. 1) than in the model without paraphoton.

It was noted in Ref. , that there is a part of the parameter space for the millicharged particles where they do not decouple from the acoustic oscillations of the baryon-photon plasma at recombination, and it was suggested that the effect of these particles on the CMB anisotropy spectrum may be similar to the effect of baryons. The purpose of this Letter is, using the recent precise CMB data from WMAP , to set an upper limit on the fraction of millicharged particles in CDM and to narrow the allowed window for millicharged particles. Assuming the standard Big Bang nucleosynthesis (BBN) value for the baryon abundance, $`\mathrm{\Omega }_bh_0^2=0.0214\pm 0.0020`$ , we arrive at the following constraint on the millicharged particle abundance,

$$\mathrm{\Omega }_{mcp}h_0^2<0.007(95\%\mathrm{CL}),$$
(1)
if millicharged particles are coupled to baryons at recombination. The latter condition is satisfied on the right of the dark solid line in Fig. 1. We see that the upper limit (1) apllies in the whole allowed window for millicharged particles in models without paraphoton. It is also worth stressing that in models with paraphoton the domain of applicability of the upper bound (1) does not depend on the value of $`\alpha ^{}`$. Using the Lee-Weinberg formula  for the relic abundance one then translates the upper bound (1) into the lower limit on the annihilation cross section of millicharged particles. As a result, the upper bound (1) excludes most part of the allowed window for the millicharged particles in models without paraphoton, leaving a small allowed region with masses in the range $`m10^110`$ GeV and charges in the range $`q10^310^1`$. In models with paraphoton, the bound (1) translates into meaningful limit on the parameters of millicharged particles as shown in Fig. 1; the excluded region depends on the value of $`\alpha ^{}`$.

Let us proceed to the derivation of the upper bound (1). To calculate the CMB anisotropy spectrum we adapt the CMBFAST code  which solves numerically the set of kinetic equations  for the linear perturbations in the primordial plasma. To take into account the presence of millicharged particles we extend this set by adding the kinetic equations for the millicharged component, modify the equations for the baryon component to take into account the elastic scattering off millicharged particles and include millicharged component contribution to the energy-momentum tensor. The rest of the perturbation equations are the same as in Ref. . The Compton scattering off millicharged particles is negligible, since the corresponding cross section is suppressed by the fourth power of the charge $`q`$.

We work in synchronous gauge and consider primordial plasma in the expanding Universe with scale factor $`a(\tau )`$ (where $`\tau `$ is a conformal time) normalized to unity at present time. Let $`T_f`$, $`\rho _f`$, $`\stackrel{}{v}_f`$ be the temperature, density and velocity of the $`f`$-th component of the plasma. In particular, $`f=e,b,\gamma ,mcp`$ for electrons, baryons, photons and millicharged particles, respectively. In what follows bar denotes space averaging. The standard variables describing fluid perturbations are $`\delta _f(\stackrel{}{k},\tau )=[\rho _f(\stackrel{}{k},\tau )\overline{\rho }_f(\tau )]/\overline{\rho }_f(\tau )`$ and $`\theta _f(\stackrel{}{k},\tau )=ik_jv_f^j(\stackrel{}{k},\tau )`$ where $`k_i`$ is conformal momentum.

Before recombination, the interaction between non-relativistic electrons and protons is strong enough to ensure that electron and baryon components have equal velocities, $`\theta _e=\theta _b`$. This makes it possible to use tight coupling approximation and consider electrons and protons as single baryon fluid . Then the set of equations for baryons and millicharged particles reads (cf. Ref. )

$`\dot{\delta }_b`$ $`=\theta _b{\displaystyle \frac{1}{2}}\dot{h},`$
$`\dot{\theta }_b`$ $`={\displaystyle \frac{\dot{a}}{a}}\theta _b+c_s^2k^2\delta _b+{\displaystyle \frac{4\overline{\rho }_\gamma }{3\overline{\rho }_b}}an_e\sigma _T(\theta _\gamma \theta _b)`$
$`+a\mathrm{\Gamma }_{mcp}\mathrm{\Omega }_{mcp}(\theta _{mcp}\theta _b),`$ (2)
$`\dot{\delta }_{mcp}`$ $`=\theta _{mcp}{\displaystyle \frac{1}{2}}\dot{h},`$
$`\dot{\theta }_{mcp}`$ $`={\displaystyle \frac{\dot{a}}{a}}\theta _{mcp}+c_{s,mcp}^2k^2\delta _{mcp}`$
$`+a\mathrm{\Gamma }_{mcp}\mathrm{\Omega }_b(\theta _b\theta _{mcp}),`$ (3)
where $`h`$ is the longitudinal metric perturbation, dot stands for derivative with respect to the conformal time $`\tau `$; $`c_s`$, $`c_{s,mcp}`$ are the sound velocities in the baryon and millicharged components, $`n_e`$ is the number density of electrons and $`\mathrm{\Gamma }_{mcp}`$ is the velocity transfer rate for millicharged particles due to scattering off baryons and electrons. The latter is given by

$$\mathrm{\Gamma }_{mcp}=\underset{x=e,p}{}\frac{n_x}{\mathrm{\Omega }_b}\frac{}{v_{M,x}}\mathrm{\Delta }v_M𝑑\sigma _{M,x}|_{v_{M,x}=0},$$
(4)
where brackets stand for thermal averaging, $`v_{M,e(p)}`$ is relative velocity and $`\mathrm{\Delta }v_M`$ is velocity transfer in a single process of scattering; $`d\sigma _{M,e}`$ ($`d\sigma _{M,p}`$) is the Rutherford cross section for millicharged particles scattering off electron (proton).

The Rutherford cross section is singular at zero scattering angle, but due to Debye screening the integral in Eq. (4) is cut at the value of the scattering angle equal to the Debye angle, $`\theta _D=\sqrt{2\pi \alpha n_e/T^2m_e}`$. As a result one arrives at the following expression for the velocity transfer rate in the case of thermal equilibrium

$$\mathrm{\Gamma }_{mcp}=\frac{4\sqrt{2\pi }\alpha ^2q^2\rho _{crit}}{3Mm_pa^{3/2}(\tau )T_0^{3/2}}|\mathrm{ln}\theta _D|(\sqrt{\mu _{M,e}}+\sqrt{\mu _{M,p}}),$$
(5)
where $`\mu _{M,e(p)}`$ is the reduced mass of a millicharged particle and electron (proton), $`\alpha `$ is the fine structure constant and $`T_02.726`$ K is the present CMB temperature. It is straightforward to generalize Eq. (5) to nonequilibrium case when electrons, protons and millicharged particles have different temperatures. In that case the value of $`\mathrm{\Gamma }_{mcp}`$ is larger than the one given by Eq. (5), so Eq. (5) may be used as a lower estimate of the interaction rate, which is sufficient for our purposes.

We solve the system of the kinetic equations starting from the early moment of time $`\tau _i`$ and using inflationary initial conditions ,

$`\delta _{mcp}`$ $`=\delta _b={\displaystyle \frac{3}{4}}\delta _\gamma ={\displaystyle \frac{1}{2}}C(k\tau _i)^2,`$
$`\theta _{mcp}`$ $`=\theta _\gamma =\theta _b={\displaystyle \frac{1}{18}}Ck(k\tau _i)^3,`$
where constant $`C`$ determines the overall normalization.

The last terms in the r.h.s. of Eqs. (2), (3) tend to equalize the velocities of the baryon and millicharged components of the fluid. This results in the kinetic equilibrium, $`\theta _{mcp}=\theta _b`$, provided $`\mathrm{\Gamma }_{mcp}`$ is large enough. In this case, the perturbation equations are difficult to solve numerically, because the kinetic relaxation rate for the interaction of millicharged particles and baryons is much larger than the rates of other processes. To deal with this situation we make use of the zeroth order tight coupling approximation adopting the method of Peebles and Yu . Namely, we expand Eqs. (2), (3) to the zeroth order in $`\mathrm{\Gamma }_{mcp}^1`$, setting $`\theta _b=\theta _{mcp}=\theta `$. Then we exclude $`\mathrm{\Gamma }_{mcp}`$ from Eqs. (2), (3) and arrive at the following equation

$`\dot{\theta }=`$ $`{\displaystyle \frac{\dot{a}}{a}}\theta +{\displaystyle \frac{\mathrm{\Omega }_b}{\mathrm{\Omega }_b+\mathrm{\Omega }_{mcp}}}{\displaystyle \frac{4\overline{\rho }_\gamma }{3\overline{\rho }_b}}an_e\sigma _T(\theta _\gamma \theta )`$ (6)
$`+{\displaystyle \frac{\mathrm{\Omega }_bc_s^2\delta _b+\mathrm{\Omega }_{mcp}c_{s,mcp}^2\delta _{mcp}}{\mathrm{\Omega }_b+\mathrm{\Omega }_{mcp}}}k^2.`$
The CMB spectrum obtained in this approximation agrees with the solution of the original set of equations at the level of one percent provided that

$$\mathrm{\Gamma }_{mcp}(\tau _{rec})(\mathrm{\Omega }_b+\mathrm{\Omega }_{mcp})H(\tau _{rec})^1250,$$
(7)
where $`\tau _{rec}`$ is conformal time at recombination and $`H(\tau _{rec})`$ is the Hubble parameter. In what follows we discuss a region of parameters where tight coupling condition (7) is satisfied (in particular, this region covers the whole allowed window in models without paraphoton) and comment on the rest of the parameter space in due course.

To compare the results of our simulation with the CMB data we consider the flat $`\mathrm{\Lambda }`$CDM model with the number of massless neutrino species $`N_\nu =3`$. We perform a scan over the space of the cosmological models by varying parameters from minimal to maximal values as given in Table 1. All priors are at 95% CL.

We assume helium fraction at the moment of recombination $`Y_{He}=0.24`$. We have checked that reionization effects are irrelevant here, the reason being that reionization affects the spectrum only at the lowest values of the multipole moments.

For each cosmological model we calculated likelihood to the WMAP data . As a result we arrived at the limit given by Eq. (1), which means that no model exists in the considered parameter range with larger values of $`\mathrm{\Omega }_{mcp}h_0^2`$ and likelihood better than 5%. We also checked that additional CMB data from CBI  and ACBAR  do not improve the limit (1).

Using the Lee–Weinberg formula we translate the limit given by Eq. (1) into the bound on the parameter space. The corresponding excluded areas are shown in Fig. 1. In models with paraphoton we extend our analysis to the region where the tight coupling condition (7) is not satisfied. We checked that actually the upper bound (1) applies provided that the following less restrictive condition holds

$$\mathrm{\Gamma }_{mcp}(\tau _{rec})(\mathrm{\Omega }_b+\mathrm{\Omega }_{mcp})H(\tau _{rec})^12.5.$$
(8)
This inequality is true on the right of the black solid line in Fig. 1.

Coming back to the tight coupling regime, let us illustrate the upper bound (1) pictorially by generating a set of models with random parameters in tight coupling regime, assuming uniform distribution of models in the parameter space within ranges given in Table 1. For each model we calculated likelihood to the observed CMB spectrum and plot the resulting distribution of models in the ($`\mathrm{\Omega }_bh_0^2`$$`\mathrm{\Omega }_{mcp}h_0^2`$) plane in Fig. 2.

One observes that the CMB spectrum is approximately degenerate along the straight line

$$\mathrm{\Omega }_bh_0^2=0.0221.1\mathrm{\Omega }_{mcp}h_0^2.$$
(9)
This degeneracy is in agreement with the expectation of Ref. , that the effect of millicharged particles is similar to that of baryons. Thus in models with millicharged particles, the CMB data  determine actually the sum $`(\mathrm{\Omega }_b+\mathrm{\Omega }_{mcp})h_0^2=0.022\pm 0.001`$ (68% CL). Combining this value with the lower limit $`\mathrm{\Omega }_bh_0^2>0.019`$ from BBN one arrives at the upper bound very similar to Eq. (1). This serves as a qualitative explanation of our result.

Another illustration of the approximate degeneracy (9) is shown in Fig. 3, where two CMB anisotropy spectra calculated for different models on the degeneracy line are shown. One observes that the two spectra almost coincide in the region of the first and second acoustic peaks. However, the degeneracy is no longer present at higher multipoles. This is due to the fact that the electroneutrality of the plasma implies that the electron number density is proportional to the baryon density. Hence, replacing a certain amount of baryons by millicharged particles results in the enhancement of the Silk damping at small scales.

With future precise data for high values of $`l`$, one will be able to set a constraint on the value of $`\mathrm{\Omega }_{mcp}`$ using the CMB data only, without reference to BBN results. To check this we created a simulated dataset, which contains the same values of $`l`$ as in the WMAP data up to $`l=500`$, and then with the step $`\mathrm{\Delta }l=50`$ up to $`l=1600`$. The CMB anisotropy spectral coefficients $`C_l`$’s were taken from the best fit  to WMAP data. The error bars for these coefficients were assumed to be equal to cosmic variance. Repeating the above procedure for this dataset we obtained that an upper limit $`\mathrm{\Omega }_{mcp}h_0^2<0.003`$ can be placed. Further improvement of this limit turns out to be impossible due to the new approximate degeneracy,

$`\mathrm{\Omega }_bh_0^2=0.0220.65\mathrm{\Omega }_{mcp}h_0^2,`$
$`n_s=0.94+8.0\mathrm{\Omega }_{mcp},`$ (10)
arising at smaller values of $`\mathrm{\Omega }_{mcp}h_0^2`$.

To conclude, we note that when translated into the parameter space, the limit (1) is especially interesting for the models without paraphoton, where it excludes most part of the window with not very heavy particles and substantial electric charges. To completely close the window, sensitivity to millicharged particle abundance at the level of $`\mathrm{\Omega }_{mcp}h_0^2310^4`$ would be required, which cannot be achieved with future CMB experiments due to the degeneracy (10). Determination of the baryon abundance from the BBN is not accurate enough to improve the situation. Hopefully, the rest of the window will be explored by future accelerator and/or laboratory experiments.

Finally, recently it was suggested  that the flux of the 511 keV $`\gamma `$-rays from the Galactic bulge detected by the INTEGRAL satellite  may be explained by the annihilation of the $`1÷100`$ MeV dark matter particles into $`e^+e^{}`$ pairs, provided their annihilation cross section $`\sigma \beta `$ and abundance $`\mathrm{\Omega }`$ satisfy $`\left(\frac{\sigma \beta }{\mathrm{pb}}\right)\left(\frac{1\mathrm{M}\mathrm{e}\mathrm{V}}{M}\right)^2\left(\frac{\mathrm{\Omega }^2}{\mathrm{\Omega }_{CDM}^2}\right)10^{(3.5÷4.5)}`$. Intriguingly, this condition holds in the left corner of the parameter space for millicharged particles without paraphoton allowed by Eq. (1) (say, for $`q=310^3`$, $`M=100`$ MeV). One is tempted to speculate that the observation of the 511 keV line is an indication that a (small) fraction of CDM is comprised of the millicharged particles. This possibility will be checked by the future CMB data.
Acknowledgments. We would like to thank V.Rubakov and S.Sibiryakov for useful discussions. This work was supported in part by the RFBR02-02-17398 and GPRFNS-2184.2003.2 grants. The work of D.G. was also supported in part by the GPRF grant MK-2788.2003.02.


text:
# Small-core photonic crystal fibers with weakly disordered air-hole claddings

## Abstract

Motivated by recent experimental work by Folkenberg et al. we consider the effect of weak disorder in the air-hole lattice of small-core photonic crystal fibers. We find that the broken symmetry leads to higher-order modes which have generic intensity distributions resembling those found in standard fibers with elliptical cores. This explains why recently reported experimental higher-order mode profiles appear very different from those calculated numerically for ideal photonic crystal fibers with inversion and six-fold rotational symmetry. The splitting of the four higher-order modes into two groups fully correlates with the observation that these modes have different cut-offs.

: J. Opt. A: Pure Appl. Opt.

The advances in fabrication techniques and structural control of both index-guiding as well as photonic band-gap photonic crystal fibers (PCFs) is believed to improve the agreement between experiments and theory based on ideal structures (for a recent review of the field we refer to Ref.  and references therein). The recent investigation of micro-bending induced attenuation in PCFs  is an example of the successful merging of experiments and simulations as well as recent experiments  on multi-mode cut-off in small-core PCFs and the good agreement with theoretical predictions . On the other hand, in the latter experiments the observed higher-order mode profiles were very different from those of an ideal structure despite the fact that the studied fibers appeared to have highly regular air-hole lattices. However, the relative strength of disorder may still be significant and in this Letter we numerically study the effect of such variations on both fields and propagation constants.

We consider the type of PCF first studied by Knight et al.  for which the ideal structure has inversion and six-fold rotational symmetry, see Fig. 1. In the first row of Fig. 2 we show the electrical field intensity $`|𝑬_m^{(0)}|^2`$ (the superscript emphasizes the ideal or “unperturbed” structure) of the first six ($`m=1,2,3\mathrm{}6`$) eigenmodes at $`\lambda =780\mathrm{nm}`$ of an ideal PCF with $`\mathrm{\Lambda }=1.4\mu \mathrm{m}`$ and $`d/\mathrm{\Lambda }=0.67`$ corresponding to one of the fibers studied in Ref. . The results are based on a numerical fully-vectorial solution of Maxwell’s equations in a plane-wave basis with periodic boundary conditions  and basis vectors which support the inversion and six-fold rotational symmetry of the dielectric structure. For the refractive index $`n`$ we use $`n=1`$ in air and for silica we use a Sellmeier expression with $`n=1.45367`$. The fundamental modes ($`m=1`$ and 2) have close-to-Gaussian profiles .

Historically, their degeneracy has been widely debated (see Refs.  and references therein) but recently it was finally proved by Steel et al.  how their degeneracy follows from group-theory because $`𝑬_1^{(0)}`$ and $`𝑬_2^{(0)}`$ do not support the full symmetry. Recently, the degeneracies of higher-order modes were addressed by Guobin et al. , but only in the case of a PCF supporting 10 eigenmodes. The general picture is quite complicated and the degeneracies and order depend on the number of guided eigenmodes . In a picture of true eigenmodes this can be understood as follows; for a given number $`M`$ of guided eigenmodes the eigenvalues are minimized with the constraint that the eigenfields are orthonormal. When the number is changed to, say, $`M+1`$ the orthogonalization affects all $`M+1`$ eigenvalues and eigenfields and thus the picture of degeneracies may in general also change compared to the picture for $`M`$ eigenmodes.

In our case we numerically find that while $`𝑬_3^{(0)}`$ and $`𝑬_4^{(0)}`$ support the full symmetry $`𝑬_5^{(0)}`$ and $`𝑬_6^{(0)}`$ do not. Group theory thus predicts that the PCF have a set of modes which are either pair-wise two-fold degenerate ($`m=1,2`$ and $`m=5,6`$) or non-degenerate ($`m=3`$ and $`m=4`$). This contrasts the situation of a step-index fiber (SIF) in the weakly-guiding approximation where the first higher-order modes are four-fold degenerate. The 3rd and 4th modes are nearly degenerate in the sense that our numerical studies show that the splitting $`\beta _4^{(0)}\beta _3^{(0)}`$ is much smaller than the spacings $`\beta _5^{(0)}\beta _4^{(0)}`$ and $`\beta _3^{(0)}\beta _2^{(0)}`$ to the two neighboring groups of degenerate modes. The grouping of the four higher-order modes into two groups is experimentally reflected in different cut-off wavelengths of the two groups . Compared to the picture in the weakly-guiding approximation, the splitting of the higher-order modes in PCFs originates from the inversion and six-fold rotational symmetry and the magnitude of the splitting is driven by the high index contrast between air and silica. We note that for a standard fibre, a sufficiently high index contrast will eventually also split the higher-order modes.

From perturbation theory (for the problem of shifting boundaries see Ref. ) it is obvious that for guided modes the six inner air-holes are of most importance to the optical properties whereas disorder in the air-hole lattice further out in the cladding mainly affects the guided modes in terms of e.g. leakage loss . Thus, in the search of over-all qualitative effects of disorder it is sufficient to only let the six inner air-holes (see Fig. 1) have small random displacements $`\mathrm{\Delta }𝒓_i`$ ($`i=1,2,\mathrm{}6`$) with respect to the ideal positions of the air-holes as well as small random deviations $`\mathrm{\Delta }d_i`$ in their diameters. We consider the situation where the displacements and deviations in diameter are Gaussianly distributed (zero mean and widths $`\sigma _\mathrm{\Lambda }`$ and $`\sigma _d`$) and mutually uncorrelated, i.e.

$`<\mathrm{\Delta }𝒓_i>=<\mathrm{\Delta }d_i>=<\mathrm{\Delta }𝒓_i\mathrm{\Delta }d_j>=0`$
$`<\mathrm{\Delta }𝒓_i\mathrm{\Delta }𝒓_j>=\sigma _\mathrm{\Lambda }^2\delta _{ij},<\mathrm{\Delta }d_i\mathrm{\Delta }d_j>=\sigma _d^2\delta _{ij}.`$ (1)
In the following we study a large ensemble (more than 1000 members) of fibers numerically.

Since the structural imperfection is typically a small perturbation of the ideal structure we have $`<\beta _m>\beta _m^{(0)}`$. The grouping discussed above thus also holds approximately for the average propagation constants $`<\beta _m>`$ as we also have confirmed numerically (not shown). In Fig. 3 we show the distribution of the propagation constant for the first six eigenmodes at $`\lambda =780\mathrm{nm}`$ of a PCF with $`\mathrm{\Lambda }=1.4\mu \mathrm{m}`$ and $`d/\mathrm{\Lambda }=0.67`$. For the disorder we have used $`\sigma _\mathrm{\Lambda }=\sigma _d=0.02\times \mathrm{\Lambda }`$ corresponding to 2% structural variations of the six inner air holes. The fundamental mode only has a small fraction of the field near to the air-hole silica interfaces compared to higher-order modes, see row (a) in Fig. 2. From perturbation theory one thus expects the fundamental mode to be less sensitive to disorder compared to higher-order modes. Indeed, from Fig. 3 it is seen that the distribution $`P(\beta _m)`$ has a common width for $`m=1,2`$ which is significantly smaller than the common width for the higher-order modes ($`m=3,4,5,6`$).

For a non-zero $`\sigma `$, symmetries are broken and symmetry-related degeneracies are in principle lifted. For $`m=1,2`$ this splitting is often referred to as form-birefringence quantified by $`\mathrm{\Delta }n_{12}=(\beta _1\beta _2)/k`$ where $`k`$ is the free-space wave number. However, since $`<𝑬_1^{(0)}\left|\delta \epsilon \right|𝑬_1^{(0)}><𝑬_2^{(0)}\left|\delta \epsilon \right|𝑬_2^{(0)}>`$ for a scalar perturbation $`\delta \epsilon `$ there is no significant form-birefringence ($`\mathrm{\Delta }n_{12}0`$) to lowest order in perturbation theory. Indeed, for the studied disorder strength we numerically find an extremely narrow distribution. The width is vanishing (within our numerical accuracy) corresponding to a delta-function distribution $`P(\beta _1\beta _2)\delta (\beta _1\beta _2)`$.

The effect of disorder on the fields is illustrated in Fig. 2. Rows (b)-(e) show examples of the electrical field intensity for four random configurations of the inner ring of air-holes. In general, we find that the modes in the ideal structure, row (a), deform into profiles having two maxima with a node in between similarly to the deformation of the modes in elliptical core SIFs . This picture correlates perfectly with the modes observed experimentally in Ref. . Fig. 4 shows near-fields recorded at $`\lambda =780\mathrm{nm}`$ for a PCF with $`\mathrm{\Lambda }1.41\mu \mathrm{m}`$ and $`d/\lambda 0.66`$. Here, the fundamental mode, panel (a), corresponds to the linear combination $`𝑬𝑬_1+𝑬_2`$ and for the higher-order modes, panels (b) and (c) corresponds to $`𝑬𝑬_3+𝑬_4`$ and $`𝑬𝑬_5+𝑬_6`$, respectively.

In conclusion we have demonstrated how even weak disorder in the air-hole lattice leads to deformation of the higher-order modes in small-core photonic crystal fibers. In the presence of disorder the higher-order modes resemble those in elliptical core standard fibers. Our findings explain the recently reported higher-order mode profiles and the grouping of higher-order modes correlates with the observation that they have different cut-offs.

We acknowledge useful discussions with B. T. Kuhlmey. M. D. Nielsen and K. P. Hansen are financially supported by the Danish Academy of Technical Sciences and J. Lægsgaard by the Danish Technical Research Council.


text:
# The intrinsic intensity modulation of PSR B1937+21 at 1410 MHz

## 1 Introduction

This letter reports on observations of the radio pulsar PSR B1937+21 taken with a center frequency of 1410 MHz at the 305m Arecibo radio telescope in Puerto Rico. Previous observations of this source at 430 MHz revealed a behavior that was completely different from other known sources (Jenet et al., 2001). At 430 MHz this source exhibits almost no detectable pulse-to-pulse fluctuations. Occasional bursts of radio radiation, or “giant pulses”, are observed but they are restricted to small regions in pulse phase (Kinkhabwala & Thorsett, 2000; Cognard et al., 1996). The observations presented here show that this remarkable pulse-to-pulse stability of the non-giant pulse emission also occurs at 1410 MHz. Edwards & Stappers (2003) attempted to measure the modulation properties of PSR B1937+21 at 1300 MHz with the Westerbork synthesis radio telescope. Unfortunately, the signal-to-noise ratio in their observations was not high enough to measure the pulse-to-pulse properties.

Observations of bright pulsars have shown that the shapes and intensities of individual pulses are unique, although they average together to form a stable mean profile. The characteristic widths of individual pulses, typically referred to as sub-pulses, are usually smaller than the average profile width. This typical behavior has not been detected in the pulses from PSR B1937+21. Instead, this source appears to be highly stable in the sense that each pulse is identical to the next.

It will be shown that the stability of this source at both frequencies can be understood within the context of the sparking gap model of pulsar emission originally proposed by Ruderman & Sutherland (1975) and later reworked in more detail by Gil & Sendyk (2000, hereafter GS00). In this model, the stability is a result of a large number of emission regions, or “sparks”, occurring on the polar cap. The GS00 model predicts that if the emission is stable at a given frequency, then it will be stable at all higher frequencies. Since this stability will also depend on the velocity of the emitting plasma, an upper bound may be placed on its Lorentz factor.

The sparking-gap model of Ruderman & Sutherland (1975) and its revised form (GS00) have received much attention in recent years. They have been used to interpret the sub-pulse properties of slow pulsar conal emission (Edwards, Stappers, & van Leeuwen, 2003; Rankin & Ramachandran, 2003; van Leeuwen, Stappers, Ramachandran, &Rankin, 2003; Asgekar & Deshpande, 2001; Deshpande & Rankin, 2001) and of millisecond pulsars core emission (Edwards & Stappers, 2003). These models have also been used in pulsar population studies (Arzoumanian, Chernoff, & Cordes, 2002; Fan, Cheng, & Manchester, 2001). Given their widespread use, it is important to develop observational tests which can determine the validity of these models. Since the GS00 model predicts that the radio emission from PSR B1937+21 should be stable at all frequencies higher than 430 MHz, high frequency observations of this source will test the current form of the GS00 model.

The next section describes the 1410 MHz observations performed at the Arecibo radio observatory. It also describes the techniques used to measure the intrinsic fluctuation properties of this source. Section 3 discusses the sparking gap model in more detail with emphasis on the frequency dependence of the modulation index. This letter is summarized in section 4.

## 2 Observations and Analysis

The data were taken at the 305 m Arecibo radio telescope using the 1400 MHz Gregorian dome receiver. The observation lasted two hours and approximately 4.5 million pulses were obtained. The exact center frequency used was 1410 MHz. Both circular polarizations were two-bit complex sampled at a rate of 10 MHz (100 ns) and recorded to tape using the Caltech baseband recorder (Jenet et al., 1998). Further processing of the data were performed at the Caltech Center for Advanced Computation and Research (CACR) using a 256 processor Hewlett-Packard Exemplar. The two-bit complex samples were unpacked and assigned optimum values in order to minimize signal distortion (Jenet & Anderson, 1998). The dual polarization voltage data was adjusted using an empirically derived cross-talk matrix (Stineberg, 1982). The effect of the Earth’s motion around the Sun was removed by resampling the complex voltage data at a rate necessary to transform the data into the barycentric frame. This rate was calculated using the software package TEMPO111http://pulsar.princeton.edu/tempo. The effects of interstellar dispersion were removed by coherently dedispersing the data (Jenet et al., 1998; Hankins & Rickett, 1975) at a dispersion measure of 71.0249 pc $`\text{cm}^3`$ (Jenet et al., 2001).

As with the 430 MHz observations, the single pulses from PSR B1937+21 at 1410 MHz are not bright enough to use conventional single pulse analysis techniques. Hence, the low intensity analysis methods used by Edwards & Stappers (2003) and Jenet et al. (2001) are used here. Typically, a pulsar exhibits pulse-to-pulse shape variability and intensity fluctuations. Three statistical techniques are used to determine if these phenomena are occurring in this source. First, the phase resolved modulation index is calculated in order to determine if any pulse-to-pulse amplitude fluctuations are present. Second, the average intensity fluctuation spectrum (Backer, 1973) is calculated to determine the time scale of the modulation. Third, the average autocorrelation function (ACF) and the ACF of the average profile are calculated and compared to each other in order to determine if there are any pulse-to-pulse shape variations occurring. Pulse-to-pulse fluctuations may be intrinsic to the source or due to propagation through the ISM. The magnitude of the modulation index, the time scale of the modulation, and the existence of any pulse shape variations can be used to discriminate between intrinsic and ISM induced pulse-to-pulse variations.

The average pulse profile, $`<I_s(\varphi )>`$, normalized by its peak value is shown in figure 1. $`I_s(\varphi )`$ is the signal intensity at pulse phase $`\varphi `$ and the angle brackets, $`<>`$, represent averaging over the ensemble of 4.5 million single pulses. The average profile consists of two components separated by about 520 milli-periods (mP). Note that only a small temporal region of length 425 $`\mu `$s (273 mP) is shown around each component with a time resolution of .42 $`\mu s`$. The giant pulses occur between 35 and 40 mP near the primary component and between 560 and 565 mP in the secondary component. The autocorrelation function of the average profile, $`C_{<I_s>}(\mathrm{\Delta }\varphi )`$, is calculated for each component separately and compared to the average autocorrelation function, $`<C_{I_s}(\mathrm{\Delta }\varphi )>`$, of each component. These correlation functions are defined as:

$`C_{<I_s>}(\mathrm{\Delta }\varphi )`$ $`=`$ $`{\displaystyle \frac{1}{\varphi _1\varphi _0}}{\displaystyle _{\varphi _0}^{\varphi _1}}<I_s(\varphi )><I_s(\varphi +\mathrm{\Delta }\varphi )>𝑑\varphi `$ (1)
$`<C_{I_s}>(\mathrm{\Delta }\varphi )`$ $`=`$ $`{\displaystyle \frac{1}{\varphi _1\varphi _0}}<{\displaystyle _{\varphi _0}^{\varphi _1}}I_s(\varphi )I_s(\varphi +\mathrm{\Delta }\varphi )𝑑\varphi >,`$ (2)
where $`\varphi _0`$ and $`\varphi _1`$ are the starting and ending pulse phases of the given emission component, respectively, and $`\mathrm{\Delta }\varphi `$ is the phase lag. Note that cyclic boundary conditions are assumed when evaluating the above expressions for the case when $`\varphi +\mathrm{\Delta }\varphi `$ lies outside the interval $`[\varphi _0,\varphi _1]`$. If the individual pulses have identical shapes from pulse-to-pulse, then these two autocorrelation function will have identical shapes. Figure 2 shows both autocorrelation functions plotted together for each component. When calculating the average ACF, $`<C_{I_s}>`$, the high resolution ACF was rebined to a time resolution of 3.2 $`\mu s`$ (2.05 mP) in order to increase the signal-to-noise ratio. All the of ACFs were normalized by their value at $`\mathrm{\Delta }\varphi =2.05`$ mP. This corresponds to the second phase-lag bin of the rebined average ACF and was chosen to avoid the noise spike occurring in the first bin of the rebined average ACF. The solid line represents $`C_{<I_s>}(\mathrm{\Delta }\varphi )`$ while the filled circles represent $`<C_{I_s}(\mathrm{\Delta }\varphi )>`$. Each polarization was analyzed separately and the resulting ACFs were averaged together. The error bars where calculated using the number of points averaged together within each bin and the root-mean-square (RMS) of the data within each bin. The error bars represent the $`95\%`$ confidence range. Note that even though the giant pulses are included in the ACF analysis, they will have little effect since they are very narrow and they occur infrequently.

The phase resolved modulation index is defined as

$$m(\varphi )=\frac{\sqrt{<I_s(\varphi )^2><I_s(\varphi )>^2}}{<I_s(\varphi )>}.$$
(3)
$`m(\varphi )`$ is a measure of the pulse-to-pulse variation in intensity. Again, both polarization were analyzed separately and the resulting modulation indices were averaged together. In order to increase the signal-to-noise ratio, both $`<I_s(\varphi )>`$ and $`<I_s(\varphi )^2>`$ were rebined to a time resolution of 3.3 $`\mu s`$ (2.1 mP) before calculating $`m(\varphi )`$. The error bars where calculated using the number of points averaged within each bin and the RMS values of $`<I_s(\varphi )>`$ and $`<I_s(\varphi )^2>`$ within each bin. The error bars represent the $`95\%`$ confidence range. If the amplitude of each pulse is constant and the only remaining fluctuations are due to a Gaussian statistical process, then $`m(\varphi )`$ will be equal to 1222Note that this definition of $`m`$ is slightly different from the typical definition which averages down the noise before calculating $`m`$. In that case, $`m=0`$ if no pulse-to-pulse fluctuations were present.. Figure 3 plots $`m(\varphi )`$ for both components. Only the phase region where a significant measurement could be made is shown. The modulation index is consistent with being constant across each component as expected if all pulses have the same shape. The average value of the modulation index is 1.37 $`\pm 0.06`$ and 1.54 $`\pm 0.2`$ for the main and secondary components, respectively.

Since pulse-to-pulse modulation is detected, it is necessary to determine if it is intrinsic to the pulsar or a result of propagation through the ISM. The expected modulation index due to the ISM alone is given by (Jenet et al., 2001):

$$m=\sqrt{1+\frac{2}{N}},$$
(4)
where $`N`$ is the number of scintills in the band. $`N`$ may be estimated by (Cordes et al., 1990):

$$N=1+\eta \frac{B}{\delta \nu },$$
(5)
where $`B`$ is the observing bandwidth, $`\delta \nu `$ is the decorrelation bandwidth and $`\eta `$ is a packing fraction which lies between 0.1 and 0.2. Using a decorrelation bandwidth of 1.1 MHz at 1410 (Cordes, 1986), m is expected to be within 1.38 and 1.48. The two dashed lines in Figure 3 represent these two values. The time scale of the fluctuations were determined by calculating the average intensity fluctuation spectrum for both pulse components. Each spectrum is consistent with zero power down to a frequency of $`0.02`$ $`s^1`$. At this point, the spectrum starts to rise. This type of spectrum is indicative of the short time scale fluctuations induced by the ISM. The characteristic time scale as determined by that frequency where the power falls to half its maximum value is given by 125 s, consistent with that measured by Cordes (1986). Since the magnitude and time scale of the observed fluctuations are consistent with that expected for the ISM, it is concluded that these fluctuations are primarily due to ISM propagation and not due to intrinsic pulse-to-pulse fluctuations.

## 3 Intensity modulation and the sparking gap model

These results can be understood in the context of the GS00 model. In this model, a region of low charged particle density forms just above the surface of the polar cap. This “vacuum gap” is postulated to contain a number of localized plasma discharges or “sparks”. The plasma generated in each spark moves upwards following the local magnetic field. At some point, processes within the spark associated plasma generate the observed radiation with higher frequencies being emitted at lower altitudes (Melikidze, Gil & Pataraya, 2000; Gil, Lyubarski & Melikidze, 2003). The minimum distance between two adjacent sparks on the pulsar’s surface is approximately given by the height of the vacuum gap, $`h`$ (see GS00 for details). Assuming that the polar cap is populated as densely as possible by these sparks, the total number of sparks on the cap is estimated by $`a^2`$ where $`a`$ is given by

$$a=\frac{r_p}{h},$$
(6)
$`r_p=10^4R_6^{1.5}P^{0.5}`$ cm is the radius of the polar cap, $`R_6`$ is the neutron star radius in units of $`10^6`$ cm, and $`P`$ is the pulsar period in seconds. $`a`$ is called the complexity parameter since it was shown by GS00 that $`a`$ is correlated with the profile morphology described by Rankin (1983). Conal single profiles have the lowest values of $`a`$, followed by conal doubles and multiples and then by triple and core singles. Since $`a`$ is proportional to the number of sparks that would emit into the observers line of site, the pulse-to-pulse modulation is expected to be anticorrelated with $`a`$. Possible evidence for this correlation is discussed in (Jenet & Gil, 2003). PSR B1937+21 has the highest complexity parameter of all pulsars whose single pulse properties are known. Hence, its stability is a result of the large number of sparks occurring on its polar cap.

A large number of sparks on the surface is not quite enough to guarantee that the individual pulses will be stable. The emission from each spark-associated plasma column must overlap each other in order to average out the fluctuations. Since the spark plasma is moving out along the magnetic field at relativistic speeds, the angular width of the radiation pattern, $`\mathrm{\Theta }_r`$, will be given by

$$\mathrm{\Theta }_r=\frac{1}{\gamma },$$
(7)
with $`\gamma =1/\sqrt{1(v/c)^2}`$, $`v`$ is the velocity of the spark associated plasma in the radio emission region, and $`c`$ is the speed of light. As the spark plasma moves up along the magnetic field lines, its radiation will be beamed in a direction tangent to the local magnetic field line. Hence, the radiation from two sparks that originate near each other on the polar cap will be emitted into different directions. For a dipolar field, it turns out that the angle, $`\mathrm{\Delta }\theta `$, between the radiation beams from two adjacent sparks will be given by

$$\mathrm{\Delta }\theta =\frac{0.015}{R_6}\frac{\mathrm{\Delta }d}{r_p}r_6^{\frac{1}{2}}P^{\frac{1}{2}},$$
(8)
where $`\mathrm{\Delta }d`$ is the distance between adjacent sparks on the polar cap (see GS00 for details), $`r_6`$ is the emission altitude normalized by the neutron star radius $`R=R_610^6`$ cm. In order for the emission from adjacent sparks to overlap, the following condition must be met

$$\mathrm{\Delta }\theta <\mathrm{\Theta }_r.$$
(9)
When the above equation is satisfied, the emission from adjacent sparks will overlap and the corresponding radio emission will show very little pulse-to-pulse intensity fluctuations provided that $`a`$ is large. From Equations 6, 7, and 8 one can rewrite equation 9 as

$$\gamma <\frac{aR_6}{0.015}\left(\frac{P}{r_6}\right)^{\frac{1}{2}},$$
(10)
where $`\mathrm{\Delta }d`$ was taken to be equal to $`h`$. If the emission from a given altitude satisfies the above condition, then emission from lower altitudes will also satisfy this condition. Hence, if emission from a given altitude is observed to be stable, then emission from all lower altitudes will also be stable. Observations of other pulsars point to the existence of a “radius-to-frequency” mapping which states that higher frequencies are emitted at lower altitudes (e.g. Kijak & Gil, 2003). Assuming that this holds for PSR B1937+21, the pulses emitted at 1410 MHz are expected to be at least as stable as those emitted at 430 MHz.

Since the observations are in agreement with the sparking gap model described above, one can go one step further and use Equation 10 to place an upper bound,$`\gamma _{max}`$, on the Lorentz factor of the emitting plasma. The following form of the radius-to-frequency map will be used (Kijak & Gil, 2003):

$$r_6(\nu _{GHz})=40\nu _{GHz}^{0.26}\dot{P}_{15}^{0.07}P^{0.3},$$
(11)
where $`\nu _{GHz}`$ is the observing frequency in GHz and $`\dot{P}_{15}`$ is the pulsar period derivative in $`10^{15}`$ s/s. For PSR B1937+21, $`r_6(0.43)=3.8`$ and $`a=23`$. Hence, Equation 10 yields $`\gamma _{\text{max}}=31R_6`$. This estimate assumes the standard vacuum gap model described by RS75. For the case of the near threshold vacuum gap (NTVG) model (Gil & Melikidze, 2002), $`a`$ can be up to about 3 times higher then that estimated from equation 6 since $`a`$ will depend on the actual value of the surface magnetic field line radius of curvature. In this case, $`\gamma _{\text{max}}=91R_6`$.

## 4 Summary and Discussion

Radio observations of the millisecond pulsar PSR B1937+21 were taken at a center frequency of 1410 MHz. Statistical techniques were used to determine if this source exhibits any of the standard single pulse phenomenology that has been seen in other radio pulsars. Remarkably, none of the usual phenomena were detected. The data are consistent with each pulse having the same shape and intensity. Previously, Jenet et al. (2001) showed that this source exhibits the same single pulse stability at 430 MHz. This stability is currently unique to this source. Hence, PSR B1937+21 maybe the prototype for a new class of radio pulsars. On the other hand, the sparking gap model provides a unified framework that can explain this source’s behaviour as an extreme case of the normal pulsar emission process. In the context of the GS00 model, the stability is a result of a large number of sparks occurring on the polar cap together with a low Lorentz factor of the emitting plasma. Since the 430 MHz observations show no pulse-to-pulse fluctuations, the GS00 model predicts that all higher frequencies will be just as stable provided higher frequencies are emitted at lower altitudes. This result is independent of the exact physical model used to describe the vacuum gap. Future observations of this source at several frequencies above 430 MHz will further test the GS00 model. If fluctuations are seen, then the current form of the sparking gap model is incorrect. Observations of this source at frequencies below 430 MHz will enable a better measurement of the Lorentz factor. If pulse-to-pulse fluctuations are seen at lower frequencies, a lower bound may be placed on $`\gamma `$ for a given model of the vacuum gap.

Part of this research was performed at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. GJ acknowledges the support of the Polish State Committee for scientific research under grant 2 P03D 008 19.