text:
# Probing Planck-scale physics with Extragalactic Sources?

## Abstract

At Planck-scale, spacetime is “foamy” due to quantum fluctuations predicted by quantum gravity. Here we consider the possibility of using spacetime foam-induced phase incoherence of light from distant galaxies and gamma-ray bursters to probe Planck-scale physics. In particular, we examine the cumulative effects of spacetime fluctuations over a huge distance. Our analysis shows that they are far below what is required in this approach to shed light on the foaminess of spacetime.

Planck-scale physics, quantum foam-induced phase incoherence of light, gamma-ray bursts, distant galaxies

It is generally believed that quantum gravity, the synthesis of quantum mechanics and general relativity, predicts that spacetime becomes “foamy” or “fuzzy” at the Planck scale given by the Planck time $`t_P=(\mathrm{}G/c^5)^{1/2}10^{44}s`$, Planck length $`l_P=ct_P10^{33}cm`$, and Planck energy $`E_P=\mathrm{}/t_P10^{28}eV`$. The fuzziness of spacetime leads to uncertainties in distance ($`l`$) measurements whose absolute value is given by $`\delta ll_P(l/l_P)^{1\alpha }`$ (similar uncertainties for time measurements) and uncertainties in energy ($`E`$) measurements given by $`\delta EE(E/E_P)^\alpha `$ (again $`\delta E`$ is an absolute value). Similar uncertainties for momentum measurements also obtain.(Ng & van Dam, 1994) The parameter, $`\alpha 1`$, specifies different quantum gravity models. The standard choice(Misner et al, 1973) of $`\alpha `$ is $`\alpha =1`$; the choice of $`\alpha =2/3`$ appears(Ng & van Dam, 2000; Ng, 2001) to be consistent with the holographic principle(’tHooft, 1993; Susskind, 1995) and black hole physics; $`\alpha =1/2`$ corresponds to the random-walk model found in the literature(Amelino-Camelia, 1994).

The ultra-high energy cosmic ray threshold anomalies (see, e.g., Ref. (Amelino-Camelia & Piran, 2001; Ng et al., 2001; Amelino-Camelia et al., 2002; Aloisio et al., 2002) and references therein) may have given us some tentalizing hints of Planckian fluctuations. But so far we lack direct experimental evidence.

Recently Lieu and Hillman(Lieu & Hillman, 2002, 2003) and then Ragazzoni, Turatto, and Gaessler (Ragazzoni et al., 2003) proposed a technique that has hitherto been overlooked to directly test Planck scale fluctuations. They argued that these fluctuations can cumulatively lead to a complete loss of phase coherence for radiation that has propagated a sufficiently large distance. As a result of this inferred phase scrambling attributed to Planck scale uncertainties (spacetime foam) these authors concluded that distant compact radiation sources should not produce the normal interference patterns (e.g., Airy rings) which are often observed. In this Letter, we examine their very interesting idea and extend their proposal to include the use of gamma-ray burst interferometry.

Since the cumulative effect of quantum foam on the phase coherence of light from extragalactic sources figures crucially in the method, we begin by establishing how large the cumulative effects are, especially for the three quantum gravity models mentioned above.

Consider a distance $`L`$ (the reader may anticipate this to denote the distance between the extragalactic source and the telescope). Let us divide it into $`L/\lambda `$ equal parts each of which has length $`\lambda `$ (the reader may anticipate this to denote the wavelength of the observed light from the distant source). In principle, any length suitably larger than $`l_P`$ can be used as $`\lambda `$. But, as shown below, the wavelength of the observed light is the natural choice. We already know that the absolute value of the uncertainty in distance $`L`$ is given by $`\delta Ll_P(L/l_P)^{1\alpha }`$.

But let us calculate it again by starting with $`\delta \lambda `$ for each part and then adding up the contributions to $`\delta L`$ from the $`L/\lambda `$ parts in $`L`$. In doing so we will find out how large the cumulative effects from the various $`\delta \lambda `$’s are. To gain insight into the process, let us first consider the $`\alpha =1/2`$ random-walk model of quantum gravity. In this model, the typical quantum fluctuation $`\delta \lambda `$ for each segment of length $`\lambda `$ is of order $`l_P(\lambda /l_P)^{1/2}`$ and takes on $`\pm `$ sign with equal probability. (In the terminology used in Amelino-Camelia et al. 2002, this is referred to as a “non-systematic” effect of quantum gravity.) To simplify this part of the argument, let us assume that $`\delta \lambda `$ takes on only two values, viz. $`\pm l_P(\lambda /l_P)^{1/2}`$ (instead of, say, a Gaussian distribution about zero, which is more likely). If the fluctuations from the different segments were all of the same sign (i.e., coherent), then together they would contribute $`\pm l_P(\lambda /l_P)^{1/2}\times (L/\lambda )`$ to $`\delta L`$. But both these two cases, yielding a linear $`L`$-dependence for $`\delta L`$, are ruled out for incoherent quantum gravity fluctuations. Each has a probablity of $`1/2^{L/\lambda }1`$ for $`(L/\lambda )1`$. (We note that $`(L/\lambda )10^{30}`$ for the example involving the active galaxy PKS1413+135 considered below.) Clearly we have here a one-dimensional random walk involving $`L/\lambda `$ steps of equal size ($`\delta \lambda `$), each step moving right or left (corresponding to $`+`$ or $``$ sign) with equal probability. The result is well-known: the cumulative fluctuation is given by $`\delta L\delta \lambda \times (L/\lambda )^{1/2}`$ which is $`l_P(L/ł_P)^{1/2}`$ as expected for consistency.

Alternatively, we can derive this cumulative factor $`(L/\lambda )^{1/2}`$ for the random walk model of quantum gravity by simply using the expressions for $`\delta \lambda `$ and $`\delta L`$ themselves. The cumulative factor $`𝒞_\alpha `$ defined by

$$\delta L=𝒞_\alpha \delta \lambda ,$$
(1)
is given, for the random walk ($`\alpha =1/2`$) case, by

$$𝒞_{\alpha =1/2}=\frac{\delta L}{\delta \lambda }=\left(\frac{L}{\lambda }\right)^{\frac{1}{2}},$$
(2)
since $`\delta Ll_P(L/l_P)^{1/2}`$ and $`\delta \lambda l_P(\lambda /l_P)^{1/2}`$. The result is as expected.

But while we have some intuitive understanding of the cumulative factor $`𝒞_{\alpha =1/2}`$ for the random walk model, viz., why it goes as the square-root of the number $`L/\lambda `$ of “steps”, we have much less intuition for the other cases. Nevertheless, since we know $`\delta L`$ and $`\delta \lambda `$, we can use Eq. (1) to find the cumulative factors, with the results

$$𝒞_\alpha =\left(\frac{L}{\lambda }\right)^{1\alpha },$$
(3)
in particular,

$$𝒞_{\alpha =2/3}=(L/\lambda )^{1/3},𝒞_{\alpha =1}=(L/\lambda )^0=1,$$
(4)
for the holographic, $`\alpha =2/3`$, case and the “standard”, $`\alpha =1`$, case respectively. Note that $`𝒞_{\alpha =1}=1`$ is independent of $`L`$. Strange as it may seem, the result is not unreasonable if we recall, for this “standard” model, $`\delta ll_P`$, independent of $`l`$. The crucial point to remember is that, for the three quantum gravity models under consideration, none of the cumulative factors is linear in $`(L/\lambda )`$. (In fact, according to Eq. (3), the cumulative effects are linear in $`L/\lambda `$ only for the physically unacceptable case of $`\alpha =0`$ for which $`\delta ll`$.) To obtain the correct cumulative factor (given by Eq. (3)) from what we may inadvertently think it is, viz., $`(L/\lambda )`$ (independent of $`\alpha `$), we have to put in the correction factor $`(\lambda /L)^\alpha `$.

With the correct cumulative factors for the various quantum gravity models at hand, we can now examine the prospect of probing Planck-scale physics by observing very distant sources. Consider the phase behavior of radiation with wavelength $`\lambda `$ received from a celestial source located at a distance $`L`$ away. Fundamentally, the wavelength defines the minimum length scale over which physical quantities such as phase and group velocities (and hence dispersion relations) can be defined. Thus, the uncertainty in $`\lambda `$ introduced by spacetime foam is the starting point for our analysis. A wave will travel a distance equal to its own wavelength $`\lambda `$ in a time $`t=\lambda /v_g`$ where $`v_g`$ is the group velocity of propagation, and the phase of the wave consequently changes by an amount

$$\varphi =2\pi \frac{v_pt}{\lambda }=2\pi \frac{v_p}{v_g},$$
(5)
(i.e., if $`v_p=v_g,\varphi =2\pi `$) where $`v_p`$ is the phase velocity of the light wave. Quantum gravity fluctuations, however, introduce random uncertainties into this phase which is simply

$$\delta \varphi =2\pi \delta \left(\frac{v_p}{v_g}\right).$$
(6)
We now argue that, due to quantum fluctuations of energy-momentum(Ng & van Dam, 1994), $`\delta EE(E/E_P)^\alpha `$ and $`\delta pp(pc/E_P)^\alpha `$, the standard radiation dispersion relation $`E^2c^2p^2=0`$ should be changed to $`E^2c^2p^2E^2(E/E_P)^\alpha `$. Recalling that $`v_p=E/p`$ and $`v_g=dE/dp`$, we obtain

$$\delta \left(\frac{v_p}{v_g}\right)\pm \left(\frac{E}{E_P}\right)^\alpha =\pm \left(\frac{l_P}{\lambda }\right)^\alpha ,$$
(7)
where we have used $`E/E_P=l_P/\lambda `$. We emphasize that this may be either an incremental advance or a retardation in the phase.

In travelling over the macroscopically large distance, $`L`$, from source to observer an electromagnetic wave is continually subjected to random, incoherent spacetime fluctuations. Therefore, by our previous argument, the cumulative statistical phase dispersion is $`\mathrm{\Delta }\varphi =𝒞_\alpha \delta \varphi `$, that is

$$\mathrm{\Delta }\varphi =2\pi a\left(\frac{l_P}{\lambda }\right)^\alpha \left(\frac{L}{\lambda }\right)^{1\alpha }=2\pi a\frac{l_P^\alpha L^{1\alpha }}{\lambda },$$
(8)
where $`a1`$. This, we believe, is our fundamental disagreement with the Lieu $`\&`$ Hillman paper, where they assume that the microscale fluctuations induced by quantum gravity into the phase of electromagnetic waves are coherently magnified by the factor $`L/\lambda `$ (see their equation (11)) rather than $`(L/\lambda )^{1\alpha }`$.

In stellar interferometry, following Lieu and Hillman’s(Lieu & Hillman, 2002, 2003) reasoning, for light waves from an astronomical source incident upon a two element interferometer to subsequently form interference fringes, it is necessary that $`\mathrm{\Delta }\varphi 2\pi `$. But the analysis of the principles of interferometry of distant incoherent astronomical “point” sources can be tricky. The local spatial coherence across an interferometer’s aperature for photons from a distant point source (i.e., plane waves) is a reflection of the fact that all photons have the same resultant phase differences across the interferometer. However, this local coherence can be lost if there is an intervening medium such as a turbulent plasma or spacetime foam capable of introducing small changes into the “effective” phases of the photon stream falling on the interferometer. Such spacetime foam-induced phase differences are themselves incoherent and therefore must be treated with the correct cumulative factors $`𝒞_\alpha `$ appropriate for the quantum gravity model under consideration.

Fluctuations due to quantum gravity are very minuscule, so they can be detected only if there is a huge cumulative effect from “summing” up the individual fluctuations. But since the cumulative factor for the “standard” model of quantum gravity (for which $`\alpha =1`$) is 1, i.e., there is no cumulative effect, obviously the proposed approach (of applying spacetime fluctuations on the phase coherence of light from extragalactic sources to probe the graininess of spacetime) cannot be used to rule out (or confirm) the $`\alpha =1`$ model. This is the first result of this Letter.

To rule out models with $`\alpha <1`$, the strategy is to to look for interference fringes for which the phase coherence of light from the distant sources should have been lost (i.e., $`\mathrm{\Delta }\varphi 2\pi `$) for that value of $`\alpha `$ according to theoretical calculations. Consider the example cited by Lieu and Hillman(Lieu & Hillman, 2003), involving the clearly visible Airy rings in an observation of the active galaxy PKS1413+135 ($`L`$ = 1.216 Gpc) by the HST at $`\lambda =1.6\mu m`$ wavelength(Perlman et al., 2002). For this example, Eq. (8) yields $`\mathrm{\Delta }\varphi 10\times 2\pi a`$ for the random walk $`\alpha =1/2`$ model and $`\mathrm{\Delta }\varphi 10^9\times 2\pi a`$ for the holography $`\alpha =2/3`$ model. Since we expect $`a1`$, the observation of Airy rings in this case would seem to marginally rule out the random walk model. (But of course proponents of the random walk model can equally claim that their favorite model is still marginally acceptable.) On the other hand, the holography model is obviously not ruled out. This finding contradicts the conclusion reached recently by Lieu and Hillman(Lieu & Hillman, 2003) who argued that the HST detection of Airy rings from PKS1413+135 has ruled out a majority of modern models of quantum gravity, including the “standard” $`\alpha =1`$ model. (Earlier, Lieu and Hillman(Lieu & Hillman, 2002) had claimed to have ruled out the $`\alpha =2/3`$ model by noticing that interference effects were clearly seen in the Infra-red Optical Telescope Array(van Belle et al., 2002) at $`\lambda =2.2\mu m`$ light from the star S Ser which is $`1`$ kpc away.) The resolution of this disagreement lies in the fact that Lieu and Hillman neglected to take into account the correction factor in estimating the cumulative effects of spacetime foam. This neglect resulted in their overestimate of the cumulative effects by a factor $`(L/\lambda )^\alpha `$: for the case of PKS1413+135, $`10^{15}`$ for $`\alpha =1/2`$, $`10^{20}`$ for $`\alpha =2/3`$, and $`10^{30}`$ for $`\alpha =1`$ respectively. Ragazzoni, Turatto, and Gaessler(Ragazzoni et al., 2003) also assumed that the cumulative factor is $`(L/\lambda )`$ rather than the correct factor $`(L/\lambda )^{1\alpha }`$. We do not agree with their claim either that the $`\alpha =2/3`$ model and the $`\alpha =1`$ model are ruled out. We note that Coule(Coule, 2003) has independently pointed out that “Planck scale is still safe from stellar images” using another argument.

Now consider a delta function type pulse of radiation from a source at a distance $`L`$. This pulse will spread in time because of quantum gravity fluctuations and the overall time dispersion in the pulse can be simply related to the aforementioned phase dispersion, i.e.,

$$\mathrm{\Delta }T\mathrm{\Delta }\varphi \left(\frac{\lambda }{2\pi v_p}\right)=\left(\frac{l_P}{L}\right)^\alpha \frac{L}{v_p}.$$
(9)
The width of the pulse increases with distance as $`L^{1\alpha }`$ but is independent of the wavelength (i.e., it is not dispersive in frequency space). For example, consider gamma ray bursts at a cosmological distance $`L10^{28}`$ cm. Then the overall time dispersion in the pulse is only given by an unobservably small $`\mathrm{\Delta }T10^{17}\times (10^{61})^\alpha `$ sec. Thus GB’s also do not offer a promising venue for testing for quantum foam, even at high Z’s.

In conclusion, we have examined the possibility of using spacetime foam-induced phase incoherence of light from distant galaxies and gamma- ray bursters to probe Planck-scale physics. These effects are real and are magnified over the large distances traversed by radiation from distant extragalactic sources. However, the effects of spacetime foam are incoherent and do not grow linearly with distance, instead increasing as $`L^{1\alpha }`$. As a result, the cumulative effects of spacetime fluctuations on the phase coherence of light are too small to be observable. Therefore, we do not conclude that modern theories of quantum gravity have been observationally ruled out.

This work was supported in part by the US Department of Energy and by the Bahnson Fund at UNC. YJN thanks D.H. Coule, E.S. Perlman, and G.T. van Belle for useful email correspondence.


format: text

text:
# Primordial Nucleosynthesis as a Test of the Friedmann Equation in the Early Universe

## I Introduction and Conclusions

The framework of our present understanding of the evolution of the universe is based on the hot big bang model. Assuming the properties of homogeneity and isotropy of the universe, the space-time metric has to be the Robertson-Walker metric

$$ds^2=dt^2a^2(t)\left\{\frac{dr^2}{1kr^2}+r^2d\theta ^2+r^2\mathrm{sin}^2\theta d\varphi ^2\right\}$$
(1)
Here, $`a(t)`$ is the scale factor, and $`k=+1,0`$ or $`1`$ corresponds to a closed, flat or open universe. The other fundamental ingredient for the dynamics of the universe are the Einstein equations. These equations, along with the metric (1), and in a universe containing a perfect fluid of total energy density $`\rho `$ and pressure $`p`$, imply the Friedmann equation (FE),

$$H^2=\frac{8\pi G}{3}\rho \frac{k}{a^2}$$
(2)
with $`H\dot{a}/a`$ the Hubble expansion rate. One also deduces the energy-momentum conservation equation

$$\dot{\rho }=3(\rho +p)H$$
(3)
To have a predictive model we need, in addition to Eqs. (2) and (3), to use other inputs, like the density of baryons, the density of dark matter and energy, etc. However, Eqs. (2) and (3) are very basic equations directly determining the evolution of the universe. In this sense, it is worth to test them as much as we can, since they could be altered by the presence of new physics Binetruy:1999ut ; Cline:1999ts ; Rafelski:2002gq ; Giudice:2000ex ; Freese:2002sq . In principle, we do not expect a radical change from the form (2), but the possibility of some alteration is open.

In this article we will concentrate on the FE (2), specifically on whether the evolution of the universe at its earliest times is given exactly by the FE. We will constrain potential modifications to the FE in the period of big bang nucleosynthesis (BBN), since this is the earliest epoch where we have observational information. We will show that we can test this equation with the present available data.

The period of BBN starts at $`T1`$ MeV with the decoupling of the weak interaction processes that allow $`np`$ inter-conversion and ends with the nuclear fusion production of the light elements, at about $`T50`$ keV.

It is an easy exercise to particularize the FE (2) in the radiation era, where BBN happens. For $`Tm_\mu `$, the relevant energy density appearing in (2) is given by

$$\rho =\rho _o\rho _\gamma +\rho _{\nu _e}+\rho _{\nu _\mu }+\rho _{\nu _\tau }+\rho _e$$
(4)
In (4) we have the contribution from photons,

$$\rho _\gamma =g_\gamma \frac{d^3k}{(2\pi )^3}\frac{E}{\mathrm{exp}(E/T)1}=2\frac{\pi ^2}{30}T^4$$
(5)
with the degrees of freedom $`g_\gamma =2`$ and $`E=|\stackrel{}{k}|`$. We also have the contribution of the three neutrino species (neutrinos plus antineutrinos) and of electrons plus positrons

$$\rho _i=g_i\frac{d^3k}{(2\pi )^3}\frac{E}{\mathrm{exp}[(E\mu _i)/T]+1}+g_i\frac{d^3k}{(2\pi )^3}\frac{E}{\mathrm{exp}[(E+\mu _i)/T]+1}$$
(6)
with the degrees of freedom $`g_\nu =1`$ and $`g_e=2`$, $`E=\sqrt{k^2+m_i^2}`$, and we allow for a particle-antiparticle asymmetry by introducing a chemical potential $`\mu _i`$.

We write the FE in the radiation dominated epoch as

$$H^2=\frac{8\pi G}{3}\rho _o$$
(7)
with $`\rho _o`$ given in (4) and where we have neglected the $`k`$-curvature term in (2).

The article is organized as follows. First, in section II, we add a new contribution to $`\rho _o`$ and constrain it when the neutrino/antineutrino asymmetry vanishes, $`\mu _\nu =0`$. We show that the new contribution is tightly constrained. In particular we follow the suggestion of the model presented in Rafelski:2002gq , and introduce modifications in $`G`$. More specifically, we let the second and third generation neutrinos to enter in Eq. (7) with a coupling $`G_2`$ and $`G_3`$, respectively, different from $`G`$. As a consequence of our general study, we can place strong constraints on this model. In section III we allow for $`\mu _\nu 0`$ and show that one may have larger departures from Eq. (7). However, even in this general case, we will be able to rule out the model of Rafelski:2002gq .

## II Testing the Friedmann equation with nondegenerate BBN

Modifications to the FE arise in a variety of contexts. For example, in brane-world models with extra dimensions one obviously finds departures from standard general relativity, and hence the FE is altered. One may have radical changes in this equation when compared to the standard equation (2) Binetruy:1999ut  but then it is difficult to reconcile the theoretical predictions with observations. Thus, the hard work is in the struggle to find scenarios where one gets an evolution equation that at first approximation is close to the FE, with modifications that can be tuned to be small Cline:1999ts . How small they must be ? We will give an answer using certain type of alterations of the FE (7), and constraining them in the radiation period using BBN. In this section we set $`\mu _{\nu _i}=0`$, but in the next section we will relax this condition.

Apart from brane-world models, we have also found inspiration in two other recent papers, where the general FE (2) is assumed to be valid, but where there are contributions to $`\rho `$ not having exactly the form (5), valid for radiation. First, in Giudice:2000ex , the authors show that the reheating temperature $`T_r`$ can be as low as $`T_r0.7`$ MeV (when BBN starts). For $`T`$ above $`T_r`$ we would have the decay of coherent oscillations of the inflaton, or the decay of a modulus or similar field. If $`T_r`$ indeed happens to be as low as speculated in Giudice:2000ex , it is conceivable that in the period of BBN, for $`T<T_r`$, apart from having radiation dominating the energy density $`\rho `$, we might still have some traces of the period $`T>T_r`$ amounting to some small contribution to $`\rho `$ that should be added to the r.h.s. of Eq. (7). Second, in another context, the so-called Cardassian expansion model has been shown to lead to a model in which the universe is flat, matter dominated and accelerating Freese:2002sq . This model simply adds a new contribution to the energy density on the r.h.s. of the FE (2), which would dominate at a late epoch in the evolution of the universe. In our article we consider similar modifications but, as we said, we constrain them in the BBN period.

With all this in mind, we add a new term to Eq. (7)

$$H^2=\frac{8\pi }{3}G\left[\rho _o+\lambda \left(\frac{T}{0.1\mathrm{MeV}}\right)^\gamma \right]$$
(8)
with $`\gamma `$ and $`\lambda `$ two arbitrary parameters, that we shall bound. Even if the type of modification that we introduce in (8) is arbitrary, some values of the parameter $`\gamma `$ have a particular physical significance. For example, the contribution of a term with $`\gamma =0`$ would correspond to a cosmological constant. For $`\gamma =4`$, the new term would be of the form of extra relativistic particles, and for $`\gamma =3`$ it would correspond to extra massive particles.

The new term in the FE modifies the expansion rate, which plays a crucial role in the dynamics during BBN. The universe, just before this period, consists of a hot dense plasma made of $`n`$, $`p`$, $`e^\pm `$, $`\nu `$ and $`\gamma `$. Neutrons and protons are in equilibrium through weak interactions. This equilibrium is lost when the universe is expanding so fast that particles cannot interact with each other. It happens when the expansion rate $`H`$ equals the weak interaction rate $`\mathrm{\Gamma }_W`$, a process that is known as the freeze out of the weak interaction. In general, any modification of the expansion rate changes the time of the freeze out.

After freeze out, the temperature of the cooling universe will eventually be low enough to permit the formation of the lightest nuclides. In order to calculate the abundances of these nuclides is important to know how many neutrons and protons are there when the nuclear reactions begin. The value $`n/p`$ has its equilibrium value, which decreases with time $`(m_n>m_p)`$, until the freeze out of weak interactions. Since then, neutrons are not destroyed or created by two body reactions, but can still be destroyed by neutron decay. Anyway, to know how many neutrons are there, one has to know at which temperature the weak reactions freeze out, because depending on the freeze out temperature, there will be more or less neutrons available to form the lightest elements.

Before we do the numerical analysis, we can understand some qualitative features of the effect of the new term in (8). We can see that for positive values of $`\lambda `$, its effect is to increase the expansion rate of the universe. Then, as we have said, the weak interaction will freeze out before the time of the freeze out in the standard case, and there will be more neutrons present in the universe at the moment the nuclides begin to form. So, the presence of a term with $`\lambda >0`$ will increase the light element yields .

The effect of a $`\lambda <0`$ term would be the contrary, but in that case there is the problem that for certain values of $`\lambda `$ and $`\gamma `$, $`H^2`$ becomes negative, which is physically unacceptable. This is indeed what happens. If we look at (8), we can see that for large values of $`\gamma `$ ($`\gamma >4`$) this occurs at high values of the temperature, i.e. when BBN begins ($`T110`$ MeV). On the other hand, for low or negative values of $`\gamma `$, $`H^2`$ becomes negative at low temperature, i.e. when BBN finishes ($`T1\mathrm{keV}`$). Anyway, for all values of $`\gamma `$, there are problems when $`\lambda <0`$. We could be able to study particular cases of $`\gamma `$ where $`\lambda <0`$ but, since we are interested in a systematic study, we will concentrate only on positive values of $`\lambda `$.

We have modified the standard numerical code Kawano , performing the modification (8) and we have obtained the new predicted values for the several light element yields. An input that is needed for our calculation is the baryon-to-photon number ratio $`\eta _Bn_B/n_\gamma `$. We use values consistent with the analysis of the cosmic microwave background anisotropies, where the parameter $`\eta _B`$ plays a role CMB ; Sarkar:2002er .

In Figs. 1 and 1, we can see the results for the deuterium and 4He abundances as a function of $`\lambda `$, for several values of $`\gamma `$. The horizontal line corresponds to the observational upper limit Sarkar:2002er . We are presenting the results for deuterium and 4He because these two elements constrain more $`\lambda `$ and $`\gamma `$ than the other elements. We can see that the abundances grow with $`\lambda `$, in agreement with what we have reasoned before.

In Fig. 2, the observationally allowed region for $`\lambda `$ and $`\gamma `$ is presented. Solid, dotted, and dashed lines correspond to the constraints coming from 4He, deuterium, and 7Li, respectively. The intersection of the constraints, giving the allowed region is shadowed. From this figure, we can see that the constraint coming from 4He is important for larger values of $`\gamma `$ and the one coming from deuterium is more relevant for smaller values of $`\gamma `$. We can understand this feature in the following way. When the nuclear reactions that produce the light elements begin, the first element formed from neutrons and protons is deuterium. Through 3He and tritium, all the deuterium converts into 4He, which is the most stable nuclide. It is when the reactions that convert deuterium into 4He freeze out that primordial deuterium begins to form. Therefore, primordial deuterium is formed after all the primordial 4He is formed. For large values of $`\gamma `$, the modification $`\lambda T^\gamma `$ we have introduced is important at the beginning of BBN (high temperature), i.e., during the formation of 4He. For small values of $`\gamma `$, the modification is more important at the end of BBN (low temperature), i.e., during the formation of deuterium (when all the 4He is already formed). This feature can also be seen in Fig. 1. In Fig. 1, the maximum allowed value of $`\lambda `$, given by the observational constraint from deuterium, is smaller for low values of $`\gamma `$ and grows as $`\gamma `$ becomes larger. The case of 4He, Fig. 1, is similar for the low values of $`\gamma `$, but for $`\gamma >3`$ the effect of the term $`\lambda T^\gamma `$ is more important, and the maximum allowed value of $`\lambda `$ is more constrained.

From Fig. 2 we see that there are strong bounds on $`\lambda `$, so we conclude that departures from the FE are very much constrained. To weight the importance of the new term $`\lambda \left(T/0.1\mathrm{MeV}\right)^\gamma `$ in (8) we can compare it numerically with the standard contribution $`\rho _o`$ in (7). For the maximum allowed values of $`\lambda `$, and making the comparison at the standard BBN temperature $`T=0.1\mathrm{MeV}`$, we find that the new term does not exceed $`10^1`$ times the standard term.

Carroll and Kaplinghat Carroll:2001bv  have done a similar study, motivated by the same concern than we have discussed in this section. They use the following generalization of the FE

$$H=\left(\frac{T}{1\mathrm{MeV}}\right)^\alpha H_1$$
(9)
and bound the parameters $`\alpha `$ and $`H_1`$ using BBN. The form (9) is suitable to test radical departures from the FE as suggested by some work in brane-world models Binetruy:1999ut . However, brane-world models developed in a series of papers Cline:1999ts  that come closer to the standard FE are, in our opinion, more promising. Thus, we think that, in this sense, our parameterization is more suitable. More support to our point of view and the form we use for the modified FE is the fact that, as it follows from our study, possibles departures from the FE have to be small due to the observational data on BBN. Even if both parameterizations are not mathematically completely equivalent, one may establish connections between the parameterization in Eq. (9) and ours, in Eq. (8), when using suitable expansions in some limits. In this sense, our work complements the study in Carroll:2001bv .

Another related work has recently been done by Zahn and Zaldarriaga Zahn:2002rr . They show how future cosmic microwave background radiation experiments will probe the FE (2). If we take their modification and test it in the radiation era using BBN (it would correspond to the particular case $`\gamma =4`$ in our Eq. (8)) we reach bounds that are on the same order of magnitude that the values found in Zahn:2002rr .

### A Particular Case : Non-Universal Gravity

In a recent paper Rafelski:2002gq , Rafelski discusses a model where one would expect the fine structure constant $`\alpha `$ to vary in time. His work has been prompted by the recent measurements of $`\alpha `$ at cosmological time scales that may show this time variation effect Webb:2000mn . He proposes that gravity couples much strongly (a factor $`10^3`$) to the second and third generation than it does to the first. In Rafelski:2002gq  it is claimed that the model is not excluded, since there is a lack of laboratory experimental constraints on how gravity couples to particles of the second and third generation (in the laboratory one measures the value of the gravitational strength through observations with matter of the first generation, namely electrons, protons and neutrons.)

However, the claim in Rafelski:2002gq  is partially not true. While certainly the laboratory constraints on how muons, taus, and second and third generation quarks couple to gravity are practically absent, surprisingly enough there are strong constraints on the neutrino sector. Gasperini Gasperini:zf  demonstrated that if gravity is not universally coupled to leptonic flavors, the gravitational field may contribute to neutrino oscillations. Laboratory experiments on neutrino oscillations place very strong bounds on the strength of the violation of universal gravity Gasperini:zf ; Minakata:1994kt .

BBN also offers a test of this exotic possibility, since muonic and tauonic neutrinos participate actively in the expansion of the early universe. The FE in a model of non-universal gravity reads

$$H^2=\frac{8\pi }{3}\left[G(\rho _\gamma +\rho _e+\rho _{\nu _e})+G_2\rho _{\nu _\mu }+G_3\rho _{\nu _\tau }\right]$$
(10)
where $`G_2`$ and $`G_3`$ are the gravitational constants for the second and third generation and $`G`$ is the usual Newton gravitational constant. As before, we have neglected the curvature term, which is not important in the radiation dominated epoch.

We will now show that with some approximations we can use the general results that we have presented before in this section. First, we consider non-degenerate neutrinos (in the next section we shall drop this assumption). In this case, and if we ignore the small temperature difference between $`\nu _e`$, $`\nu _\mu `$, and $`\nu _\tau `$ we have $`\rho _{\nu _e}=\rho _{\nu _\mu }=\rho _{\nu _\tau }\rho _\nu `$, with

$$\rho _\nu =\frac{7}{8}\frac{\pi ^2}{30}T^4$$
(11)
Then we can write (10) as

$$H^2=\frac{8\pi G}{3}\left[\rho _o+\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G}\rho _\nu \right]$$
(12)
where $`\rho _o`$ is given in (4) and we have defined $`\mathrm{\Delta }G_2`$ and $`\mathrm{\Delta }G_3`$ as $`\mathrm{\Delta }G_iG_iG`$. Since here we are interested in getting an approximate bound on $`\mathrm{\Delta }G_2`$ and $`\mathrm{\Delta }G_3`$, we consider that the temperature of the neutrinos and the photons are of the same order. Thus we can see that Eq. (12) is a particular case of the general modification (8) with $`\gamma =4`$ and

$$\lambda =\frac{7}{8}\frac{\pi ^2}{30}\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G}(0.1\mathrm{MeV})^4$$
(13)
Now, from Fig. 2, in the case $`\gamma =4`$, we can estimate the limit $`\lambda 10^5\mathrm{MeV}^4`$, which constrains the deviation from universal gravity to values of the order

$$\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G}1$$
(14)
But before comparing this result with what is needed in Rafelski:2002gq , we will see in the next section that if we allow for a neutrino degeneracy, all these constraints become weaker.

## III Testing the Friedmann equation with degenerate BBN

Up to this point we have taken $`\mu _{\nu _i}=0`$, so we have not considered a neutrino/antineutrino asymmetry. In this section, we would like to study how the results change in the case of degenerate neutrinos, i.e., $`\mu _{\nu _i}0`$. To reduce the number of parameters, we put $`\mu _{\nu _e}=\mu _{\nu _\mu }=\mu _{\nu _\tau }\mu _\nu `$, namely that the asymmetry between neutrinos and antineutrinos is the same for the three families. Actually, it has been recently pointed out Dolgov:2002ab  that this condition is a consequence of oscillations among the different families of neutrinos. Because of the dependence of $`\mu _\nu `$ on the temperature, it is convenient to use the ratio $`\xi _\nu \mu _\nu /T_\nu `$ which is a constant quantity (again, we ignore a small temperature difference between $`\nu _e`$ and $`\nu _\mu ,\nu _\tau `$).

The introduction of a neutrino degeneracy strongly affects the production of $`{}_{}{}^{4}\mathrm{He}`$. This can be easily understood because for this element, following a standard qualitative analysis Kolb:vq , one can get an approximate analytic expression for $`Y_p`$, the primordial $`{}_{}{}^{4}\mathrm{He}`$ mass fraction. It reads

$$Y_p=\frac{2\left(n/p\right)_{Nuc}}{1+\left(n/p\right)_{Nuc}}$$
(15)
where $`\left(n/p\right)_{Nuc}`$ is the ratio between neutrons and protons at the moment the $`{}_{}{}^{4}\mathrm{He}`$ is formed. In the case of degenerate neutrinos this ratio is given by 111Actually, it is the degeneracy of $`\nu _e`$ that matters, but we have put $`\xi _{\nu _e}=\xi _{\nu _\mu }=\xi _{\nu _\tau }\xi _\nu `$.

$$\left(n/p\right)_{Nuc}\mathrm{exp}\left(\frac{m_nm_p}{T_{Nuc}}\xi _\nu \right)$$
(16)
where we can see that degeneracy affects $`Y_p`$ exponentially. From Eqs. (15) and (16) we can see that positive values of $`\xi _\nu `$ diminish the predicted abundance $`Y_p`$, an effect that goes in the opposite direction compared to the consequences of the modification (8), when $`\lambda >0`$. This can be seen in Fig. 3 where a value of $`\xi _\nu =0.06`$ increases the maximum allowed value of $`\lambda `$ by about one order of magnitude (such value of $`\xi _\nu `$ is consistent with recent analysis, see Dolgov:2002ab ). As expected from our discussion above, the allowed region coming from the $`{}_{}{}^{4}\mathrm{He}`$ restriction is significatively enlarged, while the restrictions coming from the other elements remain nearly constant. We conclude that for $`\xi _\nu 0`$, larger values of $`\lambda `$ are allowed, so one has much more freedom to modificate Eq. (7) than in the case $`\xi _\nu =0`$.

### Constraints on Non-Universal Gravity from Degenerate BBN

In the same way as we did in the precedent section, we could treat the model of non-universal gravity as a particular case of (8), and from Fig. 3 we could limit the parameter $`\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G}`$ for the value $`\xi _\nu =0.06`$. But here, taking profit that in this case we have fixed one parameter ($`\gamma =4`$), we can make a more general analysis of how $`\xi _\nu `$ affects the restrictions on $`\lambda `$.

We notice that Eq. (12) is still valid in the degenerate case (recall we set $`\mu _{\nu _e}=\mu _{\nu _\mu }=\mu _{\nu _\tau }\mu _\nu `$) but now, to get $`\rho _\nu `$, we have to use Eq. (6) with the temperature $`T`$ equal to the neutrino temperature, as it should be. We make a numerical analysis with the free parameters $`\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G}`$ and $`\xi _\nu `$. In Fig. 4 we show the allowed region for the parameter space $`(\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G},\xi _\nu )`$ for the value $`\eta _B=4\times 10^{10}`$. We can see that constraints coming from the observational limits of the different elements bound these parameters into a closed region. Notice that if one only considered the restrictions coming from one element, the allowed region would be no longer closed. The reason for this is that the effects of $`\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G}`$ and $`\xi _\nu `$ on each element abundance go in opposite directions and may compensate each other. As we have noted before in the case of $`{}_{}{}^{4}\mathrm{He}`$, positive values of $`\xi _\nu `$ diminish the predicted abundance $`Y_p`$. On the other hand, positives values of $`\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G}`$ increase the expansion rate of the universe, $`H`$, which leads to an increase of $`Y_p`$, compensating the change induced by $`\xi _\nu `$.

In Fig. 5 we show the allowed region for different values of $`\eta _B`$. We note that for any $`\eta _B`$, there are values of $`\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G}`$ and $`\xi _\nu `$ that give theoretical predictions compatible with the primordial element observations. Rather than give accurate bounds on $`(\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G},\xi _\nu )`$, here we are interested in pointing out that the deviation from universal gravity is roughly limited to values $`|\frac{\mathrm{\Delta }G_2+\mathrm{\Delta }G_3}{G}|20`$. At the light of this result, we conclude that the values of $`G_2`$ and $`G_3`$ that are postulated in Rafelski:2002gq , i.e., about a factor of 1,000 bigger than the value of the Newton gravitational constant (for the first family), are excluded.

Our bounds on non-universal gravity are much weaker than the bounds coming from oscillation experiments, which are at the level $`10^{14}`$ Minakata:1994kt . However, the bounds presented in this work may have some interest due to the following reasons. First of all, our limits and the limits from oscillations are obtained at very different times; BBN period ($`z10^9`$) versus today ($`z=0`$). Also, there might be instances where oscillation experiments do not bound non-universal gravity. Two such cases are discussed by Gasperini Gasperini:zf . First, if there is no mixing between neutrino weak eigenstates and the total energy eigenstates there is no gravitational induced oscillations, and thus no bounds. Second, there might be an anomalous gravitational coupling of the neutrino with respect to the charged fields but with the three $`\nu `$ coupling equally, namely, $`\nu _e`$,$`\nu _\mu `$, and $`\nu _\tau `$ all couple with strength $`G^{}G`$. Then again there are no gravitationally induced oscillations. However, our limits still apply in these two cases.

###### Acknowledgements.

We acknowledge support by the CICYT Research Project FPA2002-00648, by the EU network on Supersymmetry and the Early Universe (HPRN-CT-2000-00152), and by the Departament d’Universitats, Recerca i Societat de la Informació, Project 2001SGR00188.


text:
# 1 Fortran program for Cray-T3E parallel computing, averaging over number processors each simulating num random walkers. Without the complications of parallel computation, loops 1 and 9 as well as all commands involving shmem, barrier, node, number, avs can be omitted.

Anomalous diffusion at percolation threshold in high dimensions on $`10^{18}`$ sites

Dirk Osterkamp, Dietrich Stauffer and Amnon Aharony1

Institute for Theoretical Physics, Cologne University, D-50923 Köln, Euroland

1 visiting from School of Physics and Astronomy, Raymond and Beverly Sackler Faculty of Exact Sciences, Tel Aviv University, Ramat Aviv, Tel Aviv 69978, Israel

e-mail: osterkamp@gmx.de, stauffer@thp.uni-koeln.de, aharony@post.tau.ac.il

Using an inverse of the standard linear congruential random number generator, large randomly occupied lattices can be visited by a random walker without having to determine the occupation status of every lattice site in advance. In seven dimensions, at the percolation threshold with $`L^7`$ sites and $`L420`$, we confirm the expected time-dependence of the end-to-end distance (including the corrections to the asymptotic behavior).

Keywords: Monte Carlo, ant in labyrinth, random number generators, anomalous diffusion

For transport in disordered media, like the ant-in-the-labyrith model of random percolation , one usually first constructs the disordered medium, and then starts the transport medium. In a hypercubic lattice of linear size $`L`$ in $`d`$ dimensions, one first has to decide about the status of $`L^d`$ sites. In two dimensions this is efficient since for long enough time every site is visited. In higher dimensions, however, only a small fraction of the lattice is visited and it is more efficient to determine the status of a site only when it is visited first, keeping this status fixed for later visits. We achieve this aim here for the case that every site is randomly allowed for the random walker with probability $`p`$ and forbidden with probability $`1p`$. We apply this method to seven-dimensional percolation  at its threshold $`p_c=0.088951`$ where an infinite network of allowed sites is just possible . We found good agreement of the new method $`(L420)`$ with the traditional one $`(L=23)`$, while for smaller $`L`$ finite-size effects are seen.

Lattice sites $`i=1,2,\mathrm{}L^d`$ are numbered in helical order such that the neighbours of site $`i`$ are $`i\pm 1,i\pm L,i\pm L^2,\mathrm{},i\pm L^{d1}`$. In the normal technique, we fill one site after the other using consecutive random numbers $`0<x_i<1`$; if $`x_i<p`$ the site is allowed. With a simple linear congruential random number generator like

$$j_{i+1}=kj_i\mathrm{mod}m;m=2^n$$

for integers $`j_i`$ instead of real numbers $`x_i`$, on computers with at least $`n`$ bits per word, this method means that a step of the random walker to the right ($`ii+1`$) corresponds to a multiplication of $`j_i`$ by $`k`$, a step to the left to a multiplication  by the inverse $`k`$* such that $`kk`$* = 1 modulo $`m`$. This inverse exists if and only if $`k`$ and $`m`$ are relatively prime; then $`k`$* can easily be calculated by the extended Euclid algorithm . Tables of pairs $`k,k`$* were given by L’Ecuyer  or can be easily calculated by a Fortran program available from DO; for example $`k`$ = 16807 works with $`k`$* = 1278498327 for $`n=32`$ and $`k`$* = 4409460005719528983 for $`n=63`$. We used on a 64-bit Cray-T3E

$`k=3512401965023503517,k`$* = 3753721746144068021,

table 5 of . If a step upwards, $`ii+L`$, corresponds to $`L`$ multiplications with $`k`$ (modulo $`m`$), then a downwards step corresponds to $`L`$ multiplications with $`k`$* (modulo $`m`$), i.e. $`L`$ pseudo-divisions by $`k`$. The other directions involve higher powers of $`L`$. In this way an arbitrary walk can be followed in a reproducible way by the proper number of multiplications with $`k`$ and $`k`$\*, always using the modulo restriction. The number of sites is then only limited by the period of the random number generator; we used a multiplier  with good spectral results and maximal period $`2^{61}`$. This method can be extended to the more general generator $`j_{i+1}=kj_i+c\mathrm{mod}m`$.

We use helical boundary conditions in $`d1`$ directions and unlimited extent in the remaining direction, i.e. site $`i`$ is thought to have as its right neighbour the site $`i+1`$ even if $`i`$ is an integer multiple of $`L`$ and lies at the right boundary. And if due to a large jump the new site has $`i0`$ or $`i>L^d`$ we do not have to put $`i`$ back into the traditional interval $`1iL^d`$ through $`ii\pm L^d`$ since the new method does not actually store an array with index $`i`$ from 1 to $`L^d`$. In this sense our lattice size is infinite in one direction. (It is easiest to imagine a planar lattice, with sites numbered in a typewriter fashion.)

In our 32-bit Fortran program, instead of the modulo function we used the automatic omission of leading bits if a product of two integers gives an integer with more than 32 bits. Since Fortran does not have the unsigned long integer type of C, the resulting products between $`2^{31}`$ and $`+2^{31}`$ are often negative. We tested by comparison with the traditional method that the new technique works nevertheless. For 64 bits, we masked off the leading bit, thus working modulo $`2^{63}`$ with only positive integers. This program including all shmem commands for communication between different Cray processors still fits onto 68 lines, Fig. 1. It can also be used for less than seven dimension by reducing the parameter idim; we made test runs with $`1,100,000,000^2,\mathrm{\hspace{0.33em}1},050,000^3,\mathrm{\hspace{0.33em}32800}^4,\mathrm{\hspace{0.33em}4100}^5,\mathrm{\hspace{0.33em}1001}^6`$ sites, as well as $`180^8`$.

```
c     diffusion in d < 9 dimensions on random lattice; < 10^18 to avoid repeats
      implicit none
      real p,factor,av(30),avs(30),time
      integer L,idim,maxtime,num,i,iwalker,maxt,iter,itime
      integer idir,isum,number,node
      integer*8 ishift(0:15),mult(0:15),ip,imod,iseed,isee2,ibm,ibmj,new
      integer shmem_n_pes, shmem_my_pe, info, iadd, barrier
      common /t3e/ av, avs

      data L/420/,idim/7/,p/0.088951/,maxtime/17/,num/1000/
      data av/30*0.0/,iseed/50923/,isee2/4711/
      data imod/’7FFFFFFFFFFFFFFF’x/
      node=shmem_my_pe()
      number=shmem_n_pes()
      time=irtc()*3.33e-9
      if(node.eq.0) print *,p,L,idim,iseed,maxtime,num,number

      ibm =2*iseed-1
      ibmj=2*isee2-1
      do 1 iter=0,node
        ibm =ibm *65539
 1      ibmj=ibmj*65539
      mult(0)=3512401965023503517
      mult(1)=3753721746144068021
      do 2 i=2,idim*2-1
 2      mult(i)=iand(mult(i-2)**L,imod)
      ip=p*imod
      factor=1.0/(number*num)

      do 8 iwalker=1,num
        do 3 i=1,idim
          ishift(2*i-2)=0
 3        ishift(2*i-1)=0
 4      ibmj=iand(ibmj*16807,imod)
        if(ibmj.gt.ip) goto 4
        maxt=2
        do 8 iter=1,maxtime
            if(iter.gt.2) maxt=maxt*2
c           evaluation of distances only if time is power of two
            do 6 itime=1,maxt
 5            ibm=ibm*16807
              idir=ishft(ibm,-60)
              if(idir.ge.2*idim) goto 5
              new=iand(ibmj*mult(idir),imod)
              if(new.gt.ip) goto 6
              ibmj=new
              ishift(idir)=ishift(idir)+1
 6          continue
            isum=0
            do 7 i=1,idim
 7            isum=isum+(ishift(2*i-2)-ishift(2*i-1))**2
            av(iter)=av(iter)+isum
 8    continue

      info = barrier()
      if(node.gt.0) stop
      if (number.gt.1) then
        do 9 iadd = 1,number-1
          call shmem_get(avs,av,maxtime,iadd)
          do 9 i=1,maxtime
 9          av(i)=av(i)+avs(i)
      endif
      do 10 iter=1,maxtime
 10     print *, iter,av(iter)*factor,2**iter
      time = irtc() * 3.33e-9 -time
      print *,’CPU= ’,4*time,number
      stop
      end

```

Now we repeat the standard scaling theory  for the average end-to-end distance $`r(t,p)`$ and enlarge it by predicting correction terms.

Let $`R_s,<r_s(t)>,<r(t)>`$ be the average radius of gyration of a cluster containing $`s`$ sites, the average end-to-end distance travelled on such a cluster in time $`t`$, and this distance averaged over all starting points on all clusters. Analogous notations hold for the squared distances; these averages $`<r^2(t)>`$ over the squares do in general not scale like the squares of the averages $`<r(t)>`$. The number of clusters containing $`s`$ sites each is $`n_s=s^\tau f[(p_cp)s^\sigma ]`$ with the standard critical exponents $`\sigma =1/(\beta \delta ),\tau =2+1/\delta `$. In general dimensions $`d`$, we have $`R_s=s^{\sigma \nu }h[(p_cp)s^\sigma ]`$ with the correlation length exponent $`\nu `$, and assume

$$<r_s^2(t)>=R_s^2g(t/s^x,R_s/s^\sigma ).$$

Right at $`p=p_c`$ this assumption simplifies to

$$<r_s^2(t)>=R_s^2G(t/s^x),$$

(varying as $`R_s^2`$ for $`ts^x`$, when every site on the cluster has been visited many times, and as $`t^{2\sigma \nu /x}=t^{2/d_w}`$, similar to the anomalous behavior on the infinite cluster, for $`1ts^x`$). For $`p<p_c`$ and for sufficiently long times one has only the former behavior: the random walker has visited the whole finite cluster, and thus $`r_s(t=\mathrm{})R_s`$. Here, $`f,h,g,G,g_1,g_2,\mathrm{}`$ are suitable scaling functions.

Thus below $`p_c`$ for long times we have

$$<r(t)>=\underset{s}{}R_ssn_sg_1[(p_cp)s^\sigma ]$$

$$<r^2(t)>=\underset{s}{}R_s^2sn_sg_2[(p_cp)s^\sigma ]$$

in $`d`$ dimensions.

For $`d>6`$ the critical exponents should be those of the Bethe lattice, with $`\tau =5/2,\sigma =\nu =1/\delta =1/2,n_ss^{5/2}\mathrm{exp}[\mathrm{const}(p_cp)^2s]`$. Also, $`d_w=6,x=3/2`$. Then

$$<r(t)>\mathrm{Const}+\mathrm{}(p_cp)^{1/2}$$

$$<r^2(t)>\mathrm{log}(p_cp)+\mathrm{}$$

for $`d>6,p<p_c`$. Right at $`p=p_c`$ the powers of $`p_cp`$ are replaced by the proper powers of $`t1/(p_cp)^3`$:

$$<r(t)>\mathrm{Const}+O(1/t^{1/6})$$
$`(1)`$
$$<r^2(t)>\mathrm{log}t$$

More precisely, in $`<r^2(t)>_ss^1g(t/s^x,\mathrm{const})`$ we can approximate the sum by $`_ss^1`$ with an upper limit for $`s`$ of order $`t^{1/x}`$, giving ln(const $`t`$) or Const + ln $`t`$.

For the third-leading term in $`<r^2(t)>`$ at $`p=p_c`$, we note that the leading correction to scaling for $`d>6`$ comes from the leading irrelevant parameter $`w`$, which represents the probability of having vertices with three bonds . This implies corrections to any leading power law behavior of the form $`(1+\mathrm{const}\times w^2X^{6d})`$, where $`X`$ represents an appropriate length scale, which could be $`r(t),R_s`$ or $`t^{1/d_w}`$. Specifically, at $`p_c`$ we expect corrections to $`<r^2(t)>`$ of relative order $`t^{(6d)/6}`$, which become $`1/t^{1/6}`$ at $`d=7`$, just as for $`<r(t)>`$ above. The two corrections differ for $`d>7`$. Thus finally at $`d=7`$ we expect:

$$<r^2(t)>\mathrm{log}t+\mathrm{Const}+O(1/t^{1/6})$$
$`(2)`$
Obviously, Eq. (1) is easier to test and reasonably confirmed by Fig. 2, using $`p_c=0.088951`$ from . For Eq. (2) we see in Fig. 3 a logarithmic long-time behaviour for $`L=23`$ (traditional method) and $`L=420`$ (new method), while for short times and/or small lattices deviations exist. Figure 4, similar to Fig. 2, shows that $`<r^2(t)>=1.25\mathrm{ln}(t)7.8+8/t^{1/6}`$ is consistent with our data for large systems over 8 decades in time. ($`<r^2(t=1)=p`$ exactly, but the curve shown in Fig.4 extrapolates to a value at $`t=1`$ differing from this value by 0.2,) The fit is not shown in Fig. 3 since it would barely be distinguishable there from the curve for $`L=420`$. Less extensive simulations in eight dimensions, $`L=180`$, nicely confirm $`<r^2>\mathrm{log}(t)`$ \+ const, but are not accurate enough to distinguish between correction exponents 1/3 and 1/6. In contrast, for seven dimensions 1/6 fits over a wider range than 1/3.

In summary, the new method gives results consistent with the old one but allows for enormously larger lattice sizes; and these results are consistent with scaling theory.

We thank Humboldt Foundation and GIF for supporting this collaboration, and R.M. Ziff, P. L’Ecuyer and P.M.C. de Oliveira for encouraging discussion.


text: # 1 Characteristics of H II regions in NGC 3256.


text:
# Single spin- and chiral-glass transition in vector spin glasses in three-dimensions

## Abstract

Results of Monte Carlo simulations of XY and Heisenberg spin glass models in three dimensions are presented. A finite size scaling analysis of the correlation length of the spins and chiralities of both models shows that there is a single, finite-temperature transition at which both spins and chiralities order.

There is now general consensus that a spin glass transition occurs in three dimensional Ising spin glasses with short range interactions at finite temperature $`T_{SG}`$. The most convincing work is that of Ballesteros et al. Ballesteros et al. (2000) who performed a finite-size scaling analysis of the correlation length, $`\xi _L`$, in samples of different sizes $`L`$. Data for the dimensionless ratio $`\xi _L/L`$ is found to intersect cleanly at $`T=T_{SG}`$, as expected at a second order phase transition.

The situation is much less clear, however, for vector spin glasses. Early work on XY Morris et al. (1986); Jain and Young (1986) and Heisenberg Morris et al. (1986); Olive et al. (1986) models indicated a zero temperature transition, or possibly a transition at a very low but non-zero temperature. However, following earlier work of Villain Villain (1975), which emphasized the role of “chiralities” (Ising-like variables which describe the handedness of the non-collinear spin structures), Kawamura and Tanemura Kawamura and Tanemura (1987) proposed a chirality transition at $`T=T_{CG}(>0)`$, even though the spin glass transition temperature, $`T_{SG}`$, is assumed to be zero. This scenario requires that spins and chiralities decouple at long length scales. Kawamura and collaborators have given numerical evidence for this scenario both for XY Kawamura and Li (2001) and Heisenberg Kawamura (1998); Hukushima and Kawamura (2000) models.

However, the absence of a spin glass transition in vector spin glass models has been challenged. For the XY case, Maucourt and Grempel Maucourt and Grempel (1998) and subsequently Akino and Kosterlitz Akino and Kosterlitz (2002) found evidence for a possible finite $`T_{SG}`$ from zero temperature domain wall calculations. Furthermore, by studying the dynamics of the XY spin glass in the phase representation, Granato Granato (2000) found that the “current-voltage” characteristics exhibited scaling behavior which he interpreted as a transition in the spins as well as the chiralities. For the Heisenberg model, Matsubara et al. Matsubara et al. (2001); Endoh et al. (2001), and Nakamura and Endoh Nakamura and Endoh (2002) have argued that the spins and chiralities order at the same low but finite temperature.

Since the most successful approach to demonstrate a finite temperature transition in the Ising case has been the scaling of the correlation length Ballesteros et al. (2000) it seems useful to perform a similar analysis for vector spin glass models. Furthermore, one can calculate the correlation length for both the spins and chirality, and so perform the same analysis for both types of ordering. Here, we present results of these calculations for the XY and Heisenberg models. For both models, we find a single transition for both spins and chiralities at low but finite temperature.

We take the standard Edwards-Anderson spin glass model

$$=\underset{i,j}{}J_{ij}𝐒_i𝐒_j,$$
(1)
where the $`𝐒_i`$ are $`n`$-component vectors of unit length at the sites of a simple cubic lattice, and the $`J_{ij}`$ are nearest neighbor interactions with zero mean and standard deviation unity. We consider both the XY model ($`n`$=2), and the Heisenberg model ($`n=3`$). Periodic boundary conditions are applied on lattices with $`N=L^3`$ spins.

The spin glass order parameter generalized to wave vector $`𝐤`$, $`q^{\mu \nu }(𝐤)`$, is defined to be

$$q^{\mu \nu }(𝐤)=\frac{1}{N}\underset{i}{}S_i^{\mu (1)}S_i^{\nu (2)}e^{i𝐤𝐑_i},$$
(2)
where $`\mu `$ and $`\nu `$ are spin components, and “$`(1)`$” and “$`(2)`$” denote two identical copies of the system with the same interactions. From this we determine the wave vector dependent spin glass susceptibility $`\chi _{SG}(𝐤)`$ by

$$\chi _{SG}(𝐤)=N\underset{\mu ,\nu }{}[\left|q^{\mu \nu }(𝐤)\right|^2]_{\mathrm{av}},$$
(3)
where $`\mathrm{}`$ denotes a thermal average and $`[\mathrm{}]_{\mathrm{av}}`$ denotes an average over disorder. The spin glass correlation length is then determined Palassini and Caracciolo (1999); Ballesteros et al. (2000), from

$$\xi _L=\frac{1}{2\mathrm{sin}(k_{\mathrm{min}}/2)}\left(\frac{\chi _{SG}(0)}{\chi _{SG}(𝐤_{\mathrm{min}})}1\right)^{1/2},$$
(4)
where $`𝐤_{\mathrm{min}}=(2\pi /L)(1,0,0)`$.

For the XY model the chirality of a square is Kawamura and Li (2001)

$$\kappa _i^\mu =\frac{1}{2\sqrt{2}}\underset{l,m}{\overset{}{}}\mathrm{sgn}(J_{lm})\mathrm{sin}(\theta _l\theta _m),$$
(5)
where $`\theta _l`$ is the angle characterizing the direction of spin $`𝐒_l`$, and the sum is over the four bonds around the elementary plaquette perpendicular to the $`\mu `$ axis and whose “bottom left” corner is site $`i`$. The chiral glass susceptibility is then given by

$$\chi _{CG}^\mu (𝐤)=N[\left|q_c^\mu (𝐤)\right|^2]_{\mathrm{av}},$$
(6)
where the chiral overlap $`q_c^\mu (𝐤)`$ is given by

$$q_c^\mu (𝐤)=\frac{1}{N}\underset{i}{}\kappa _i^{\mu (1)}\kappa _i^{\mu (2)}e^{i𝐤𝐑_i}.$$
(7)
We define the chiral correlation lengths $`\xi _{c,L}^\mu `$ by

$$\xi _{c,L}^\mu =\frac{1}{2\mathrm{sin}(k_{\mathrm{min}}/2)}\left(\frac{\chi _{CG}(0)}{\chi _{CG}^\mu (𝐤_{\mathrm{min}})}1\right)^{1/2},$$
(8)
in which $`\chi _{CG}(𝐤=0)`$ is independent of $`\mu `$. Note that $`\xi _{c,L}^\mu `$ will, in general, be different for $`\widehat{\mu }`$ along $`𝐤_{\mathrm{min}}`$ (the $`\widehat{x}`$ direction) and perpendicular to $`𝐤`$. We denote these two lengths by $`\xi _{c,L}^{}`$ and $`\xi _{c,L}^{}`$ respectively.

For the Heisenberg spin glass, Kawamura Kawamura (1998) defines the local chirality in terms of three spins on a line as follows:

$$\kappa _i^\mu =𝐒_{i+\widehat{\mu }}𝐒_i\times 𝐒_{i\widehat{\mu }}.$$
(9)
The chiral glass susceptibilities and correlation lengths are then given in terms of the $`\kappa _i^\mu `$ by Eqs. (6)-(8), as for the XY model.

We use parallel tempering Hukushima and Nemoto (1996); Marinari (1998) Monte Carlo to go down to the low temperatures that are needed, and study sizes from $`L=4`$ to 12. To test for equilibration Katzgraber et al. (2001) we require that the following relation Katzgraber and Young (2002),

$$[q_lq_s]_{\mathrm{av}}=\frac{2}{z}T[U]_{\mathrm{av}},$$
(10)
valid for a Gaussian bond distribution, is satisfied. Here $`U`$ is the energy per spin, $`q_l=(1/N_b)_{i,j}𝐒_i𝐒_j^2`$ is the “link overlap”, $`q_s=(1/N_b)_{i,j}(𝐒_i𝐒_j)^2`$ where $`N_b=(z/2)N`$ is the number of nearest neighbor bonds, and $`z(=6\text{here})`$ is the lattice coordination number. We averaged over 1000 samples, except for the following cases: XY, $`L=12`$, 601 samples, Heisenberg: $`L=8`$, 436 samples, and $`L=12`$, 331 samples. The number of sweeps that each set of spins performed varied from 6000 for the small sizes to 300000 for $`L=12`$.

Since $`\xi _L/L`$ is dimensionless it has the finite size scaling form

$$\frac{\xi _L}{L}=\stackrel{~}{X}\left(L^{1/\nu }(TT_{SG})\right),$$
(11)
where $`\nu `$ is the correlation length exponent. Note that there is no power of $`L`$ multiplying the scaling function $`\stackrel{~}{X}`$, as there would be for a quantity with dimensions. There are analogous expressions for the chiral correlation lengths. From Eq. (11) it follows that the data for $`\xi _L/L`$ different sizes come together at $`T=T_{SG}`$. In addition, they are also expected to splay out again on the low-$`T`$ side if there is spin glass order below $`T_{SG}`$.

Next we discuss the results, starting with the XY spin glass. Data for the various correlation lengths are shown in Fig. 1. One sees that the chiral correlation lengths are smaller than the spin glass correlation length at the higher temperatures, but increase faster on lowering $`T`$, such that, for a given $`L`$, all the lengths become comparable at the lowest temperatures simulated. Furthermore the two chiral correlation lengths (parallel and perpendicular to $`𝐤`$) become indistinguishable at lower $`T`$ and larger sizes, as one would expect.

The data for $`\xi _L/L`$, shown in Fig. 2, intersects at a well defined temperature $`0.34`$ and splay out at lower temperatures, which, according to Eq. (11), implies a transition at this temperature. We find

$$T_{SG}=0.34\pm \mathrm{\hspace{0.17em}0.02}\text{XY spin glass}.$$
(12)
Note the $`L=12`$ data intersects at somewhat lower $`T`$, implying that corrections to finite size scaling may still be significant for this range of sizes. Fig. 2 provides compelling evidence, in our view, that there is finite spin glass transition temperature in a three-dimensional XY spin glass, in contrast to the claim in most of the literature. The inset to Fig. 2 shows the data collapses well according to Eq. (11) with $`\nu =1.2\pm 0.2`$. In this paper, error bars do not include systematic effects which are hard to estimate. Given the lower intersection point of the $`L=12`$ data, it is possible that $`T_c`$ could be lower than that estimated here, in which case the value of $`\nu `$ would be increased perhaps to the Ising valueBallesteros et al. (2000) $`2.15\pm 0.15`$.

Data for the chiral correlation lengths for the XY model is shown in Fig. 3. Though the intersections are not quite as clean, they occur at about the same temperature as found above for the spin glass correlation length. For both spin and chiral correlations, they occur at a slightly lower $`T`$ for the $`L=12`$ data. Altogether, the evidence is good that there is a finite chiral glass transition at (or very close to) $`T_{SG}`$. Collapsing the data we find the chiral correlation length exponent is $`\nu _c=1.3\pm 0.3`$, which is compatible with our estimate for the spin correlation exponent $`\nu `$. Note that if the spins order then the chiralities must also order, assuming a non-collinear state, and so $`T_{CG}T_{SG}`$.

We are not aware of any estimates of transition temperatures for the XY spin glass with Gaussian couplings, though for the $`\pm J`$ model Kawamura and Li Kawamura and Li (2001) find $`T_{CG}=0.39\pm 0.03`$, somewhat higher than ours. Since, for the Ising spin glass, $`T_c`$ is somewhat higher for the $`\pm J`$ model than the Gaussian model, our estimate of $`T_c`$ is probably compatible with Kawamura and Li’s. However, we emphasize that, in contrast to them, we find simultaneous ordering of the spins and chiralities.

Next, we go on to our results for the Heisenberg spin glass. As for the XY model, the spin glass correlation length is larger at higher temperatures but the chiral correlation length grows faster and is comparable to the spin correlation length at the lowest temperatures. Figure 4 shows data for $`\xi _L/L`$, which intersect at a common temperature indicating a finite spin glass transition temperature $`T_{SG}`$ which we estimate to be

$$T_{SG}=0.16\pm \mathrm{\hspace{0.17em}0.02}\text{Heisenberg spin glass}.$$
(13)
The inset to Fig. 4 shows the data collapses quite well according to Eq. (11) with $`\nu =1.1\pm 0.2`$.

Figure. 5 shows that the data for the chiral correlation lengths indicate a transition at about the same value, with the intersections being cleaner for the parallel than for the perpendicular correlation length. Our estimate for the chiral correlation length exponent is $`\nu _c=1.3\pm 0.3`$. As for the XY spin glass, $`T_c`$ may be somewhat lower than that found here, which would lead to a larger value of $`\nu `$.

Our value for the transition temperature agrees well with values of $`T_{CG}`$ given by Kawamura Kawamura (1998), $`0.157\pm 0.01`$, and Hukushima and Kawamura Hukushima and Kawamura (2000), $`0.160\pm 0.005`$, though, unlike those authors, we claim that the spins, as well as the chiralities, order at this temperature. For the $`\pm J`$ model, Endoh et al. Endoh et al. (2001) find a spin glass transition at $`T_{SG}=0.19\pm 0.02`$ while Nakamura and Endoh Nakamura and Endoh (2002) find both chiral and spin glass ordering for $`T0.21`$.

To conclude, by analyzing data for the spin and chiral correlation lengths, we have argued that there is a single phase transition, at which both spins and chiralities order, in the XY and Heisenberg spin glasses in three dimensions. In our view, the evidence for a spin glass transition is at least as strong as that for a chiral glass transition. Spin–chirality decoupling does not seem to occur. The present work used quite modest work station facilities, so it would be feasible to extend these results to larger sizes by a major computational effort.

Why has this simple picture of a single finite temperature transition in vector spin glasses in three dimensions not been generally accepted before? One reason is that $`T_{SG}`$ is very low compared with the mean field value $`T_{SG}^{MF}`$ ($`1.22`$ for XY and $`0.82`$ for Heisenberg). Until the advent of parallel tempering Hukushima and Nemoto (1996) it was difficult to reach the actual $`T_{SG}`$ in simulations. Furthermore, the commonly used Binder ratio does not seem to be very useful Shirakura and Matsubara (2002) for vector spins, and $`T=0`$ domain wall calculations are plagued by uncertainties over the optimal choice of boundary conditions Akino and Kosterlitz (2002). We argue that, as for the Ising spin glass Ballesteros et al. (2000), finite size scaling of the correlation length is the optimal technique to use.

###### Acknowledgements.

We would like to thank Helmut Katzgraber for helpful discussions, and both he and H. Kawamura for comments on a preliminary version of the manuscript. We acknowledge support from the National Science Foundation under grant DMR 0086287.


text:
# Josephson Effect in Pb/I/NbSe2 Scanning Tunneling Microscope Junctions

## I Introduction

Scanning tunneling microscopy (STM) has been proven to be an invaluable tool in the study of superconductors by serving as a probe of the density of states with very high energy resolution and sub-nanometer spatial resolutionHess90 ; Yazdani97 ; Pan00 ; Hoffman02 . STM has also been successfully used to create tunnel junctions on materials that do not readily form high quality planar junctions (e.g. in MgB2)Rubio01 , or in cases where a particular small region of the sample is of interestTruscott99 .

In recent years it has been realizedPan98 ; Naaman01a ; Smakov01  that the use of superconducting (SC) STM tips may provide a complementary technique to conventional STM studies, since in SC STM both quasiparticle and Josephson tunneling can be measured locally with a scannable superconductor/insulator/superconductor (S/I/S) junction. Here we report on experiments we performed using SC STM tips. We measure the tunneling characteristics of Pb and NbSe2 samples at T$``$2 K for various values of the junctions’ normal state resistances. For junction resistances below 100 k$`\mathrm{\Omega }`$ we observe fluctuation dominated Josephson effects, and estimate the Josephson IcRN product, a physically important quantity, for tunneling into NbSe2.

## II Experimental Technique

Superconducting tips are prepared according to Ref. Naaman01a,  by deposition of Pb(5000 Å)/Ag(30 Å) proximity bilayers on conventional PtIr tips. Fig. 1(a) shows a typical dI/dV curve for tunneling at 2.1 K between a SC tip and a Pb(5000 Å)/Ag(30 Å) sample deposited on a graphite substrate at the same time as the tips were made. Sharp peaks are seen at voltages $`eV=\pm 2\mathrm{\Delta }_{Pb}`$, and the Pb phonon structure is clearly observed, confirming the high quality of the junction and that single-step tunneling is the dominant conduction mechanism. At these thicknesses (30 Å) the Ag contributes very little to the junction characteristics. The method described above for SC tip fabrication is robust and reproducible, yielding over 80% working tips, with virtually indistinguishable SC properties.

The measurement of the Josephson effect in STM junctions is challenging due to their small size and the associated high junction resistance. The energy scale for the phase coupling of a Josephson junction is the Josephson binding energy, $`E_J=\mathrm{\Delta }R_Q/2R_N`$, where $`R_Q=h/4e^26.45`$ k$`\mathrm{\Omega }`$, $`\mathrm{\Delta }`$ is the SC gap, and $`R_N`$ is the normal state resistance of the junction. In junctions for which $`k_BTE_J`$, phase locking is destroyed by thermal fluctuations, and no dc supercurrent can be observed. Indeed, for Pb junctions with resistances $``$100 M$`\mathrm{\Omega }`$, the Josephson binding energy is $`E_J/k_B0.5`$ mK, and we see no trace of the dc Josephson effect at our base temperature of 2.1 K (Fig. 1(a)).

Several authors have shownIvanchenko69 ; Ambegaokar69  that when $`E_J`$ is smaller than, but comparable to $`k_BT`$, the phase dynamics in a small Josephson junction is dominated by thermal fluctuations. In this regime pair tunneling would be observed, but the pair current will be dissipative, with the voltage drop proportional to the rate of change of the relative phase across the junction. Ivanchenko and Zil’bermanIvanchenko69  found that the current-voltage ($`IV`$) relation for small $`E_J/k_BT`$ due to this incoherent pair tunneling has the form:

$$I(V)=\frac{I_c^2Z_{env}}{2}\frac{V}{V^2+V_p^2};$$
(1)
where $`I_c=2eE_J/\mathrm{}`$, $`Z_{env}`$ is the effective impedance of the junction’s environment, and $`V_p=(2e/\mathrm{})Z_{env}k_BT_n`$ with $`T_n`$ defined as the effective noise temperature.

Fig. 1(b) shows representative $`IV`$ curves taken on a Pb/Ag film at 2.1 K for various junction resistancesnote1  below 120 k$`\mathrm{\Omega }`$. A current peak in the $`IV`$ characteristics emerges from the quasiparticle background as the junction resistance is lowered. The $`IV`$ curves can be fitted with good agreement to Eq. 1 for $`|V|<`$1.0 mV, with $`AI_c^2Z_{env}/2`$ and $`V_p`$ as the only fitting parameters. The dependence of the quantity $`\sqrt{(4e/\mathrm{})A/V_p}`$ on $`G_N=1/R_N`$ should be linear with a slope equal to $`I_cR_N/\sqrt{k_BT_n}`$. Fig. 2 shows such a plot for several Pb/Ag samples and a NbSe2 sample. Knowing the noise temperature $`T_n`$, one can infer the Josephson product $`I_cR_N`$.

Fluctuation dominated dc and ac Josephson effects were observed and carefully analyzed in Ref. Naaman01b, , confirming that the features we observe in the tunneling characteristics of these STM junctions are indeed signatures of pair tunneling. We have shownNaaman01b  that the noise parameters $`T_n`$ and $`Z_{env}`$ depend only on the experimental setup and not on the sample studied or the tip used. Once these parameters are determined by measuring a sample with a known $`I_cR_N`$ (e.g. Pb), this quantity can be inferred for any sample measured under the same experimental conditionsnote4 .

## III Results and Discussion – NbSe2

We apply the method described above to the measurement of the Josephson product for tunneling between a SC STM tip and a 2H–NbSe2 crystal. 2H–NbSe2 is a layered material that undergoes a charge density wave (CDW) transition at 33 K and a superconducting transition at 7.2 K. Niobium diselenide has been studied extensively in the pastHess90 ; Truscott99 ; Yokoya01 ; Dordevic01 ; HudsonPhD , but to our knowledge, no reliable Josephson tunneling data exist for this system. While NbSe2 appears to be a conventional BCS superconductor, it shows anisotropy of the SC order parameterHudsonPhD , transport and optical propertiesDordevic01 , as well as coexisting CDW and SC states. These make NbSe2 an interesting system, especially since the high-Tc cuprates are believed to share some of these propertiesHoffman02 .

$`IV`$ curves were measured at 2.1 K on a 2H–NbSe2 single crystalOglesby94  cleaved at room temperature in vacuum ($`5\times 10^7`$ torr). For junction resistances larger than 40 k$`\mathrm{\Omega }`$, the scaled tunneling current $`I(V)\times R_N`$ is independent of $`R_N`$. We use these curves to represent the background quasiparticle current, with no Josephson contribution, which we then subtract from the $`IV`$ curves measured at lower $`R_N`$ to find the contribution of pair tunneling to the total current. This is done because junctions with NbSe2 as one of the electrodes have considerable quasiparticle weight below the gap even at this temperature of 2.1 K.

Fig. 2 shows a plot constructed from the fits to Eq. 1 for several Pb/Ag films and the NbSe2 data shown in the inset. From linear fits to the data we obtain the slopes $`0.75\pm 0.03`$ mV1/2 and $`0.46\pm 0.03`$ mV1/2 for Pb/Ag and NbSe2 respectively. The ratio of $`I_cR_N`$(Pb/I/NbSe2) and $`I_cR_N`$(Pb/I/Pb) is thus $`0.61\pm 0.05`$. Since $`I_cR_N`$(Pb/I/Pb)=1.671 mV is knownnote2  from the Ambegaokar-Baratoff (AB) formulaAmbegaokar63 , we get $`I_cR_N`$(Pb/I/NbSe2)=$`1.02\pm 0.08`$ mV.

The SC gap in NbSe2 can be described quite well by an anisotropic s-wave gap functionHudsonPhD , with $`\mathrm{\Delta }(0)`$ ranging from 0.7 to 1.4 meV over the Fermi surface (FS). Estimates using the AB formula for Pb/I/NbSe2 junctions yield $`I_c^{AB}R_N`$1.25 mV, significantly higher than the one we observe, even when using the smallest gap parameter in NbSe2. This suggests that the gap anisotropy alone cannot explain our result and we must look for additional possible explanations.

A further reduction in $`I_cR_N`$ may result from the existence of a CDW and its effect on the SC state. Gabovich et al. Gabovich90  calculated the Josephson product for CDW and SDW superconductors and found that it decreases from the AB limit by an amount that depends on the ratio between the CDW and the SC energy gaps, and the degree to which the FS is affected by the formation of a CDW gap. In light of our results it appears that explanations along these lines should be sought, especially since the effects of a CDW may be important in the interpretation of Josephson data in high-Tc cuprates.

In summary, we measured the Josephson effect in junctions formed between SC STM tips and Pb and NbSe2 samples. In NbSe2 we find the Josephson product $`I_cR_N=1.02\pm 0.08`$ mV. This result is smaller than expected from the AB formula, but may be described within the theory of Gabovich et al. for Josephson tunneling in CDW superconductors. These results suggest that the interplay between the SC state and the CDW may have an important role in Josephson tunneling.

###### Acknowledgements.

We would like to thank L. Bokacheva, W. Teizer, and the UCSD Physics Electronics Shop. This work was supported by DOE Grant number DE-FG03-00ER45853.


text:
# Untitled Document

About new quasilinear model of heat conduction with finite velocity of heat front movement.

A.N. Skripka

KCK Soft, Kiev, Ukraine

skripka@ukrpost.net

A new quasilinear mathematical model of heat conduction with finite velocity of heat front movement is offered.

1. The analysis of existing models of heat conduction.

Now as basic mathematical model in the theory of heat conduction the linear parabolic differential equation is used. In case of one spatial dimensional:

$$\frac{T}{t}=D_T\frac{^2T}{^2x}.$$
$`(1)`$
Also it can be presented as system of two equations:

$$c\rho \frac{T}{t}=div\stackrel{}{J},$$
$`(2)`$
$$\stackrel{}{J}=\stackrel{}{i}\lambda \frac{T}{x}.$$
$`(3)`$
The following designations are entered: $`T(t,x)`$ \- field of temperature, $`\stackrel{}{J}(t,x)`$ \- vector of heat flow, $`\lambda `$ \- coefficient of heat conduction, $`D_T=\lambda /(c\rho )`$ \- coefficient of temperature conduction, $`\rho `$ \- density, $`c`$ \- thermal capacity, $`div=/x`$, $`\stackrel{}{i}`$ \- unit vector.

Equation (1) is obtained taking into account suppositions about instantaneous relaxation of the heat flow, that contradicts physical conceptions about the nature of heat. This leads to infinite velocity of front heat movement in the given model.

Sometimes model the nonlinear heat conduction equation:

$$\frac{T}{t}=\frac{}{x}(D_T(T)\frac{T}{x}),$$
$`(4)`$
where the coefficient of heat conduction is function of temperature, or linear hyperbolic equation Cattaneo is used. But both these models have basic lacks, which do not allow to widely use them at the decision of practical tasks.

In this work new heat conduction equation is obtained, which is a generalization of the linear equation (1). For this purpose equation (1) is transformed to the nonlinear equation in such a manner, that quantitatively its solutions differ from solutions of equation (1) a little, but the new equation will describe the process of heat conductivity qualitatively precisely.

In the new equation there is a parameter, connected with velocity of heat conduction front.

2. The new model of heat conduction and investigation of its local groups of symmetry.

Let’s change the equation (3) and have it presented as

$$\stackrel{}{J}=\stackrel{}{i}\lambda (\frac{T}{x}aarctg(\frac{1}{a}\frac{T}{x})),$$
$`(5)`$
where $`a`$ \- parameter, which has dimension of a gradient of temperature. From equation (5) it is clear, that flow $`\stackrel{}{J}`$ on the absolute value never exceeds the flow described in equation (3) more than $`0.5a\pi `$ and relative value of the flows differences approaches zero, when temperature gradient increases tending to infinity or this gradient being close to zero. Equation (5) and equation (2) will give the new equation of heat conduction:

$$\frac{T}{t}=D_T\frac{(T/x)^2}{(T/x)^2+a^2}\frac{^2T}{x^2}.$$
$`(6)`$
Equation (6) is invariant concerning the local groups of symmetry:

$$X_1=_t,X_2=_x,X_3=_T,X_4=2t_t+x_x+T_T.$$
$`(7)`$
Using invariants $`T/\sqrt{t}=const`$ and $`x/\sqrt{t}=const`$ generated by the operator of dilatation $`X_4`$ it is possible to obtain replacement of variables $`z=x/\sqrt{t},`$ $`T(x,t)=\sqrt{t}f(z).`$ Substituting them into equation (6), we obtain the ordinary differential equation for function $`f(z):`$

$$2D_T\frac{^2f}{z^2}(\frac{f}{z})^2=(fz\frac{f}{z})((\frac{f}{z})^2+a^2).$$
$`(8)`$
3. Numeric solution of equation (8).

Consider boundary value problem for equation (6) with initial condition $`T(t=0,x)=0,`$ boundary conditions $`T(t,x=0)=B\sqrt{t},T(t,x=0)/x=C,`$ $`x0`$ and assume, that $`D_T=1.`$

At reduction of equation (6) to equation (8) boundary conditions will be transformed to a kind of $`f(0)=B,f(0)/z=C.`$ In figure 1 results of the numerical solution of the equation (8) at $`B=1`$, $`C=1`$, $`a^2=0.001,z=0\mathrm{}5`$ are shown. Velocity of movement of heat front being $`V_0=V_0(t)`$ it is possible to obtain the following:

$$f(z=z_0)=0z_0=x_0/\sqrt{t}V_0=dx_0/dt=z_0/(2\sqrt{t}).$$

Function $`f=f(z)`$ accepts zero value at $`z1.55`$ and, consequently, velocity of movement of heat front will be $`V_01.55/(2\sqrt{t}).`$

This solution qualitatively precisely describes the process of heat conduction in real physical media.

4. Example of non-local symmetry of the equation (6).

Equation (6) by means of non-local transformation $`H(t,x)=T/x`$ can be transformed to a kind of:

$$\frac{H}{t}=\frac{}{x}\left(\frac{D_T(H/x)^2}{(H/x)^2+a^2}\frac{H}{x}\right).$$

This is a nonlinear equation of heat conduction, in which the coefficient of heat conduction is function from gradient of temperature.

5. Quasilinearity of equation (6).

Numerical researches of equation (8) have shown,that its solutions have physical sense at $`0<a0.1`$. It is possible to assume, that this parameter and in the initial equation (6) also will accept only such small values. Then equation (6) can be presented as a system of the equations:

$$\frac{T}{t}D_T\frac{^2T}{x^2},\text{ for }\left(\frac{^2T}{x^2}\right)^2>>a^2,$$
$`(9)`$
$$\frac{T}{t}0,\text{ for }\frac{^2T}{x^2}0.$$
$`(10)`$
Equations (9) and (10) are close to a linear kind and, therefore, equation (6) can be named as the quasilinear equation of heat conduction.

1.Tikhonov A.N., Samarsky A.A. Equations of mathematical physics. – Moscow: Nauka, 1966. – 660p. (in Russian)
2.Lykov A.V. Theory of heat conduction. – Moscow: Vyshaya shkola, 1967. – 559p. (in Russian)