text:
# Correlation between characteristic energies in non–𝑠-wave pairing superconductors

## 1 Introduction

In earlier work , we have discussed the possible correlation between the critical temperature $`T_c`$ in anisotropic superconductors and a ‘natural’ energy scale $`\epsilon _c\mathrm{}^2/m^{}l_c^2`$ involving the effective mass $`m^{}`$ and some characteristic length $`l_c`$. Here, by anisotropic superconductors we mean a superconductor characterized by non-spherically symmetric pairing, giving rise to an anisotropic $`𝐤`$-dependence of the gap energy in momentum space. Well-known instances of such materials are the heavy Fermion compounds, most of them being characterized by a $`p`$-wave order parameter , and the high-$`T_c`$, cuprates, whose order parameter displays $`d`$-wave symmetry . The most natural physical choice for the characteristic length $`l_c`$ entering the definition of $`\epsilon _c`$ above was the coherence length $`\xi `$ . In the case of anisotropic superconducting materials in the presence of magnetic fluctuations, we later correlated the coherence length $`\xi `$ to the spin-fluctuation temperature $`T_{\mathrm{sf}}`$ as $`k_\mathrm{B}T_{\mathrm{sf}}\mathrm{}^2/m^{}\xi ^2`$ .

The possible enhancement of the critical temperature $`T_c`$ in an anisotropic superconductor was already studied in the early work of Markowitz and Kadanoff  within a weak-coupling, BCS-like approach, and more recently revived in Refs.  as possibly relevant for the high-$`T_c`$ cuprates. A generalization of Markowitz and Kadanoff’s results to the strong and intermediate coupling regime however showed that in the strong coupling limit anisotropy is effectively averaged out, and the gap tends to become isotropic .

In the case of the high-$`T_c`$ cuprates, the strong coupling limit corresponds to the underdoped region of their phase diagram, which is usually interpreted in terms of a crossover between Bose-Einstein condensation of strongly coupled preformed pairs above $`T_c`$, and weak coupling, BCS-like superconductivity in the overdoped regime . On the other hand, the coherence length $`\xi `$ can be continuously connected to the characteristic size of the preformed bosonic pairs in the normal state of (underdoped) high-$`T_c`$ cuprates. Therefore, the coherence length may serve to parametrize the crossover from weak to strong coupling, with $`k_\mathrm{F}\xi `$ decreasing in going from weak coupling, characterized by large superconducting pair fluctuations, to strong coupling, $`k_\mathrm{F}`$ being the Fermi momentum  (see also Refs. , and refs. therein). In the underdoped regime, a relatively short coherence length is consistent with the idea of preformed pairs localized in real space . In such a limit, it is interesting to investigate how the internal structure of these preformed pairs is related to the overall symmetry of the many-body order parameter .

Here, we derive an explicit expression for the coherence length in terms of the pair binding energy for the Cooper problem in anisotropic superconductors. Such an expression \[see Eq. (13) below\] explicitly contains the quantum number $`\mathrm{}`$ of the pair relative angular momentum, which is usually employed to parametrize the anisotropic character of the order parameter, $`\mathrm{}=0,1,2`$ corresponding to $`s`$-, $`p`$-, and $`d`$-wave symmetry, respectively. While this expression correctly reduces to the standard one for isotropic $`s`$-wave superconductors, in the case $`\mathrm{}>0`$ it agrees qualitatively with the phenomenological dependence of $`k_\mathrm{B}T_c`$ on the characteristic energy $`\epsilon _c`$, proposed in Ref.  for the heavy Fermion compounds as well as for the high-$`T_c`$ cuprates. Such a dependence of the characteristic energy for superconductivity on $`\epsilon _c`$ is then qualitatively confirmed by numerically solving the self-consistent gap equations for the maximum gap at $`T=0`$ in the case of an anisotropic superconductor in the crossover between the weak- and strong-coupling limits, as a function of the dimensionless crossover parameter $`k_\mathrm{F}\xi `$, now employing a more general definition of the coherence length $`\xi `$ for anisotropic superconductors.

The paper is organized as follows. In Sec. 2, after briefly reviewing the Bethe-Goldstone equation for the Cooper problem in the case of isotropic superconductors, we generalize and solve it for non–$`s`$-wave superconductors, characterized by a superconducting instability in the $`(\mathrm{},m)`$ channel of relative angular momentum quantum numbers. Our results are discussed in relation to the earlier phenomenological findings in Ref. . In Sec. 3 we analyze the gap equations arising from the corresponding many-body problem at the mean-field level. Although we now have to resort to numerical integration (at least in the non–$`s`$-wave case), our results qualitatively confirm the approximate relation between characteristic energies derived in Sec. 2, also in the crossover from weak to strong coupling. In Sec. 4 we eventually summarize and propose some directions for future work.

## 2 Bethe-Goldstone equation for non–$`s`$-wave superconductors

The Bethe-Goldstone equation for the Cooper problem  in momentum space reads:

$$(\epsilon 2\xi _𝐤)\psi _𝐤=\underset{𝐤^{}}{}V_{\mathrm{𝐤𝐤}^{}}\psi _𝐤^{},$$
(1)
where $`\psi _𝐤`$ is the Fourier transform of the pair wave-function with wave-vector $`𝐤`$ (here, we assume a spin-singlet state, with zero total momentum), $`\xi _𝐤=\mathrm{}^2k^2/(2m^{})\mu `$ is the energy of a single electron with respect to the chemical potential $`\mu `$, $`m^{}`$ is the effective mass, $`\epsilon `$ is the binding energy of the electron pair, and $`V_{\mathrm{𝐤𝐤}^{}}`$ is the Fourier transform of the electron-electron interaction (see also Ref.  for a pedagogical review). In the weak-coupling limit, we may safely assume $`\mu =\epsilon _\mathrm{F}`$ at $`T=0`$, with $`\epsilon _\mathrm{F}`$ the Fermi energy. However, we anticipate that this identification will be superseded in Sec. 3, where the crossover from weak to strong coupling will be addressed more in detail. In the case of anisotropic superconductors, we adopt the spirit of the Anderson-Brinkman-Morel model for $`p`$-wave superfluidity in 3He , and expand the electron-electron interaction in spherical harmonics around the Fermi surface as

$$V_{\mathrm{𝐤𝐤}^{}}=\frac{1}{\mathrm{\Omega }}\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\underset{m=\mathrm{}}{\overset{\mathrm{}}{}}V_\mathrm{}mY_\mathrm{}m(\widehat{𝐤})Y_\mathrm{}m(\widehat{𝐤}^{}),$$
(2)
for $`|\xi _𝐤|`$, $`|\xi _𝐤^{}|<\mathrm{\Lambda }`$, and zero otherwise, where $`\widehat{𝐤}`$ is the unit vector pointing along the direction of $`𝐤`$, $`\mathrm{\Omega }`$ denotes the volume of the system, and $`\mathrm{\Lambda }`$ is an energy cut-off, characterized by the nature of the interaction. In the case of conventional superconductivity ($`s`$-wave pairing, or $`\mathrm{}=0`$), it would be natural to identify such an energy scale with the Debye energy, as in BCS theory.

The use of spherical harmonics to expand the electron-electron interaction in anisotropic pairing superconductors characterized by a spherical Fermi surface was earlier considered by Markowitz and Kadanoff  in the weak-coupling regime (see also Refs. ), and later by Combescot  in the strong and intermediate coupling regime. A related model of exotic Cooper pairing with finite angular momentum has been discussed also in Ref. , where rotational symmetry breaking ($`\mathrm{}>0`$) is due to the interplay between the finite range of the attractive potential and the interelectronic average distance. Within such model, the authors of Ref.  also derive the $`\mathrm{}`$-dependence of the critical temperature in the Bose-Einstein limit. In the more realistic case of a non-spherical Fermi surface, spherical harmonics are naturally replaced by Allen’s Fermi surface harmonics , which have been used by Whitmore *et al.*  in extending the results of Ref.  within the framework of Eliashberg equations for $`T_c`$. Moreover, spherical harmonics afford a natural classification of anisotropic superconductors (such as heavy Fermion compounds  as well as high-$`T_c`$ cuprates ) in terms of their pairing symmetry.

Within the weak-coupling approximation, the largest attractive coupling constant $`V_\mathrm{}mV>0`$ in Eq. (2) gives rise to a pairing instability in the $`(\mathrm{},m)`$ channel. In the following, we shall neglect other instabilities, which may arise well below $`T_c`$, corresponding to mixed symmetry pairing. This amounts to retaining the $`(\mathrm{},m)`$ term only in the expansion Eq. (2). In this case, Eq. (1) has a solution given by

$$\psi _𝐤=\alpha \frac{Y_\mathrm{}m(\widehat{𝐤})}{\epsilon 2\xi _𝐤}$$
(3)
belonging to the eigenvalue

$$\epsilon 2\mathrm{\Lambda }\mathrm{exp}\left(\frac{2}{VN(0)}\right),$$
(4)
where $`N(0)`$ is the density of states at the Fermi level, and $`\alpha `$ is a normalization constant. Such a solution corresponds to a bound state ($`\epsilon <0`$), and a further mean-field analysis of the many-electron problem (see also Sec. 3) shows that it indeed corresponds to a superconducting state, characterized by a gap function at $`T=0`$:

$$\mathrm{\Delta }_{\widehat{𝐤}}=2\mathrm{\Lambda }\mathrm{\Gamma }e^{1/N(0)V}Y_\mathrm{}m(\widehat{𝐤}),$$
(5)
having the symmetry of the attractive channel under consideration, where:

$$\mathrm{ln}\mathrm{\Gamma }=𝑑\mathrm{\Omega }_𝐤|Y_\mathrm{}m(\widehat{𝐤})|^2\mathrm{ln}|Y_\mathrm{}m(\widehat{𝐤})|,$$
(6)
the integration being carried over the unit sphere . The anisotropic $`𝐤`$-dependence of the pair wave-function $`\psi _𝐤`$ in Eq. (3) in the case $`\mathrm{}>0`$ provides interesting information on the internal structure of a Cooper pair and its connection with the overall symmetry of the many-body gap function.

Given the pair wave-function $`\psi _𝐤`$, the coherence length $`\xi `$ is naturally defined by

$$\xi ^2=\frac{_𝐤|_𝐤\psi _𝐤|^2}{_𝐤|\psi _𝐤|^2}.$$
(7)
Passing to the continuous limit, in the isotropic, $`s`$-wave case ($`\mathrm{}=0`$), one obtains

$$\xi ^2=\frac{2\mathrm{}^2}{3m^{}\epsilon _\mathrm{F}}\frac{1}{x^2},$$
(8)
where $`x=|\epsilon |/(2\epsilon _\mathrm{F})`$ measures the binding energy of a pair in units of the energy of an unbound pair. Apart from a numerical factor, from Eq. (8) one thus recovers the correct order of magnitude relation between the critical temperature, the Fermi velocity $`v_\mathrm{F}=\mathrm{}^1d\xi _𝐤/dk`$, and the coherence length:

$$k_\mathrm{B}T_c|\epsilon |\frac{\mathrm{}v_\mathrm{F}}{\xi }.$$
(9)
In the anisotropic case, for a pairing instability in the $`(\mathrm{},m)`$ channel, with the pair wave-function given by Eq. (3), the denominator in Eq. (7) is easily integrated in the continuous limit as

$`{\displaystyle \underset{𝐤}{}}|\psi _𝐤|^2`$ $``$ $`\mathrm{\Omega }{\displaystyle d^3𝐤|\psi _𝐤|^2}`$ (10)
$`={\displaystyle \frac{\mathrm{\Omega }\alpha ^2}{4\pi }}{\displaystyle _0^{\mathrm{}}}𝑑E{\displaystyle \frac{N(E)}{(\epsilon 2E)^2}}{\displaystyle 𝑑\mathrm{\Omega }_𝐤|Y_\mathrm{}m(\widehat{𝐤})|^2}`$
$`\mathrm{\Omega }\alpha ^2{\displaystyle \frac{N(0)}{8\pi \epsilon }},`$
where $`N(E)/4\pi `$ is the density of states per unit solid angle, and use has been made of the normalization condition of the spherical harmonics. In the case $`\mathrm{}0`$, the anisotropic $`𝐤`$ dependence of the pair wave-function Eq. (3) gives rise to two contributions in the numerator of Eq. (7), according to the chain rule

$$|_𝐤\psi _𝐤|^2=\frac{\alpha ^2}{k^2}\frac{|_{\widehat{𝐤}}Y_\mathrm{}m(\widehat{𝐤})|^2}{(\epsilon 2\xi _𝐤)^2}+4\alpha ^2\mathrm{}^2v_\mathrm{F}^2\frac{|Y_\mathrm{}m(\widehat{𝐤})|^2}{(\epsilon 2\xi _𝐤)^4}.$$
(11)
Here, we have made use of the decomposition $`_𝐤=k^1_{\widehat{𝐤}}+\widehat{𝐤}/k`$, where $`_{\widehat{𝐤}}`$ denotes the angular part of the gradient operator in momentum space, and of the fact that $`\xi _𝐤`$ depends only on $`k=|𝐤|`$. Passing to the continuous limit, integrating separately over the angles as in Eq. (10), and making use of the identity

$`{\displaystyle 𝑑\mathrm{\Omega }_𝐤|_{\widehat{𝐤}}Y_\mathrm{}m(\widehat{𝐤})|^2}`$ $`=`$ $`{\displaystyle 𝑑\mathrm{\Omega }_𝐤Y_\mathrm{}m^{}(\widehat{𝐤})_{\widehat{𝐤}}^2Y_\mathrm{}m(\widehat{𝐤})}`$ (12)
$`=\mathrm{}(\mathrm{}+1),`$
which follows from a variant of the Green’s formula over the unit sphere, one eventually obtains

$$\xi ^2=\frac{\mathrm{}^2}{2m^{}\epsilon _\mathrm{F}}\left[\frac{4}{3}\frac{1}{x^2}+\frac{\mathrm{}(\mathrm{}+1)}{1+x}\left(1\frac{x\mathrm{ln}x}{1+x}\right)\right].$$
(13)
Equation (13) implicitly relates the binding energy $`|\epsilon |`$ of a Cooper pair to its characteristic size $`\xi `$ for anisotropic pairing superconductors.

In the case of non–$`s`$-wave superconductors, Equation (13) manifestly includes an ‘orbital’ contribution, proportional to $`\mathrm{}(\mathrm{}+1)`$, as earlier surmised in Ref.  on the basis of phenomenological considerations. There, we proposed the existence of a phenomenological relation linking $`k_\mathrm{B}T_c`$ to the characteristic energy

$$\epsilon _c=\frac{\mathrm{}^2}{2m^{}\xi ^2},$$
(14)
in the case of anisotropic superconductors. That the effective mass $`m^{}`$ should enter inversely in determining the scale of $`k_\mathrm{B}T_c`$ was earlier recognized by Uemura *et al.* , who did not, however, include the coherence length in their analysis . By comparing the experimental values for several heavy Fermion compounds ($`p`$-wave superconductors, $`\mathrm{}=1`$) as well as high-$`T_c`$ superconductors ($`d`$-wave superconductors, $`\mathrm{}=2`$), in the anisotropic case we found that $`k_\mathrm{B}T_c=f(\epsilon _c)`$ deviates from the square-root behavior, $`k_\mathrm{B}T_c\sqrt{\epsilon _c}`$, that is easily derived from Eqs. (9) and (14) for isotropic, $`s`$-wave superconductors. In particular, in the heavy Fermion case, we noticed a large initial slope in $`f(\epsilon _c)`$ for $`\epsilon _c=0`$, and a tendency of such function to approach saturation, as $`\epsilon _c`$ increases.

In view of the fact that $`k_\mathrm{B}T_c|\epsilon |x\epsilon _\mathrm{F}`$ , Eq. (13) implicitly defines the generalization to the anisotropic case of the relationship between $`k_\mathrm{B}T_c`$ and $`\epsilon _c`$ we were looking for. Equation (13) correctly reduces to Eq. (8) in the isotropic case ($`s`$-wave superconductors, $`\mathrm{}=0`$). In the case of non–$`s`$-wave superconductors ($`\mathrm{}>0`$), Eq. (13) indeed increases with increasing $`\epsilon _c`$, starting with a logarithmically infinite slope at $`\epsilon _c=0`$, and tending to a saturation value as $`\epsilon _c\mathrm{}`$. Such a limit is *e.g.* approached for $`k_\mathrm{F}\xi 1`$. Setting $`a=4/[3\mathrm{}(\mathrm{}+1)]`$ and performing an asymptotic analysis of Eq. (13), in the limit $`k_\mathrm{F}\xi 1`$ ($`\mathrm{}>0`$) one obtains

$$x=\frac{1+a}{W[(1+a)/e]},$$
(15)
where $`W(z)`$ is Lambert’s function , and in the limit of very large anisotropy ($`\mathrm{}1`$):

$$\underset{\mathrm{}\mathrm{}}{lim}x=W^1(e^1)=3.591\mathrm{}.$$
(16)
Figure 1 shows the dependence of the dimensionless measure of the pair binding energy $`x=|\epsilon |/2\epsilon _\mathrm{F}`$ on $`\epsilon _c/\epsilon _\mathrm{F}=(k_\mathrm{F}\xi )^2`$, implicitly defined by Eq. (13) in the cases $`\mathrm{}=0,1,2`$, corresponding to $`s`$-, $`p`$-, and $`d`$-wave superconductivity, respectively. Keeping fixed all other variables, one notices that for $`\xi >\xi _0`$, where

$$\frac{\mathrm{}^2}{2m^{}\xi _0^2\epsilon _\mathrm{F}}=\frac{3}{4W^2(e^1)}=9.672\mathrm{},$$
(17)
$`|\epsilon |`$ increases as $`\mathrm{}`$ increases while, for $`\xi <\xi _0`$, $`|\epsilon |`$ decreases as $`\mathrm{}`$ increases, although without large deviations from the limiting value, Eq. (16). Therefore, at least for sufficiently weakly coupled superconductors ($`\xi >\xi _0`$), anisotropy enhances superconductivity, the degree of anisotropy being here parametrized by the order $`\mathrm{}`$ of the spherical harmonic modelling the $`𝐤`$ dependence of the order parameter.

## 3 Mean-field analysis of the many-body problem and the crossover from weak to strong coupling

The definition of the Cooper pair size $`\xi `$ we employ in Eq. (7) makes use of the pair wave-function $`\psi _𝐤`$ for the Cooper problem, in which the only many-body effect is that of forbidding the occupancy of states below the Fermi level, $`\epsilon _\mathrm{F}`$. Within BCS theory, *i.e.* at the mean field level, a self-consistent treatment of the superconducting instability affords a better definition of the pair wave-function, which in momentum space is given by  $`\psi _𝐤=\mathrm{\Delta }_𝐤/(2E_𝐤)`$, where $`\mathrm{\Delta }_𝐤`$ is the gap function, now depending on $`𝐤`$ as a vector, and $`E_𝐤=(\xi _𝐤^2+|\mathrm{\Delta }_𝐤|^2)^{1/2}`$ the upper branch of the quasiparticle dispersion relation. With this definition of $`\psi _𝐤`$, when $`E_𝐤`$ is allowed to vanish locally on the Fermi surface, as is the case for $`p`$\- and $`d`$-wave superconductors, it has been recently shown that nodal quasiparticles give rise to a logarithmically divergent contribution to $`\xi `$, as defined by Eq. (7) . While such a drawback does not arise in $`s`$-wave superconductors , in the case of anisotropic superconductors the coherence length must be defined in terms of the range in real space of the static correlation function for the modulus of the order parameter . The fact that our result Eq. (13) provides a finite estimate for $`\xi `$ also in the anisotropic case clearly depends on our approximate choice for $`\psi _𝐤`$, which solves only the two-body Cooper problem.

In order to discuss the limits of validity of the approximate results derived in Sec. 2, we will then work out numerically the mean-field solution to the corresponding BCS Hamiltonian, both in the weak and in the strong-coupling limit. We start by reviewing the model and notations set out in Refs.  (see also Ref. ).

We shall consider the following Hamiltonian for a superconducting system in three dimensions :

$$H=\underset{𝐤\sigma }{}\xi _𝐤c_{𝐤\sigma }^{}c_{𝐤\sigma }+\underset{\mathrm{𝐤𝐤}^{}𝐪}{}V_{\mathrm{𝐤𝐤}^{}}c_𝐤^{}c_{𝐤+𝐪}^{}c_{𝐤^{}+𝐪}c_𝐤^{},$$
(18)
where $`c_{𝐤\sigma }^{}`$ \[$`c_{𝐤\sigma }`$\] is the creation \[destruction\] operator for an electron with wave-vector $`𝐤`$ and spin projection $`\sigma \{,\}`$ along a specified direction, $`\xi _𝐤=\mathrm{}^2k^2/(2m^{})\mu `$ is the dispersion relation for free electrons with effective mass $`m^{}`$, measured with respect to the chemical potential $`\mu `$, and

$$V_{\mathrm{𝐤𝐤}^{}}=\frac{V}{\mathrm{\Omega }}Y_\mathrm{}m(\widehat{𝐤})Y_\mathrm{}m(\widehat{𝐤}^{})$$
(19)
is the projection of the electron-electron interaction along the $`(\mathrm{},m)`$ channel, as in Eq. (2), which we assume to be attractive ($`V>0`$) . Equation (19) generalizes to the non–$`s`$-wave case the contact potential in real space discussed *e.g.* by Marini *et al.* . Standard diagonalization techniques then lead to the mean-field coupled equations for the gap energy $`\mathrm{\Delta }_𝐤=\mathrm{\Delta }_0Y_\mathrm{}m(\widehat{𝐤})`$ and the particle density $`n`$ at $`T=0`$ (see Ref.  for the case $`T0`$):

$`{\displaystyle \frac{1}{V}}`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Omega }}}{\displaystyle \underset{𝐤}{}}{\displaystyle \frac{|Y_\mathrm{}m(\widehat{𝐤})|^2}{2E_𝐤}},`$ (20a)
$`n`$ $`=`$ $`{\displaystyle \frac{2}{\mathrm{\Omega }}}{\displaystyle \underset{𝐤}{}}v_𝐤^2,`$ (20b)
where $`v_𝐤^2=\frac{1}{2}(1\xi _𝐤/E_𝐤)`$, together with $`u_𝐤^2=1v_𝐤^2`$, are the usual coherence factors of BCS theory.

Owing to our choice of a free particle dispersion relation (all band effects are embedded in a single parameter, namely the effective mass $`m^{}`$), and of a contact potential in real space, the sums over $`𝐤`$ in Eq. (20a) above give rise to an ultraviolet divergence, which requires a suitable regularization. In three dimensions, it is customary to intoduce the scattering amplitude $`a_s`$ , which for our anisotropic interaction reads:

$$\frac{m^{}}{4\pi a_s}=\frac{1}{V}+\frac{1}{\mathrm{\Omega }}\underset{𝐤}{}\frac{m^{}}{k^2}|Y_\mathrm{}m(\widehat{𝐤})|^2.$$
(21)
Subtracting Eq. (21) from Eq. (20a) one has:

$$\frac{m^{}}{4\pi a_s}=\frac{1}{\mathrm{\Omega }}\underset{𝐤}{}\left(\frac{1}{2E_𝐤}\frac{m^{}}{k^2}\right)|Y_\mathrm{}m(\widehat{𝐤})|^2.$$
(22)
Following Ref. , we render Eqs. (20) and Eq. (22) dimensionless, by introducing the dimensionless quantities

$`x^2`$ $`={\displaystyle \frac{\mathrm{}^2k^2}{2m^{}}}{\displaystyle \frac{1}{\mathrm{\Delta }_0}},`$ $`x_0`$ $`={\displaystyle \frac{\mu }{\mathrm{\Delta }_0}},`$
$`\xi _x`$ $`={\displaystyle \frac{\xi _𝐤}{\mathrm{\Delta }_0}}=x^2x_0,`$ $`E_x`$ $`={\displaystyle \frac{E_𝐤}{\mathrm{\Delta }_0}}=\sqrt{\xi _x^2+|Y_\mathrm{}m(\widehat{𝐱})|^2},`$
and the Fermi energy $`\epsilon _\mathrm{F}=\mathrm{}^2k_\mathrm{F}^2/(2m^{})=\mathrm{}^2(3\pi ^2n)^{2/3}/(2m^{})`$. In particular, Marini *et al.*  observe that $`x_0=\mu /\mathrm{\Delta }_0`$ can be used as a parameter for the crossover between the strong-coupling, Bose-Einstein (BE) limit ($`x_00`$) and the weak-coupling, BCS limit ($`x_00`$). Moreover, on one hand, the expression resulting from the dimensionless version of Eq. (22) for $`(k_\mathrm{F}a_s)^1`$ can be inverted to obtain $`x_0`$ as a function of the dimensionless scattering length $`k_\mathrm{F}a_s`$. On the other hand, one may alternatively use $`k_\mathrm{F}\xi `$ as the independent variable in place of $`k_\mathrm{F}a_s`$ . Indeed, it was earlier recognized by Pistolesi and Strinati  (following the seminal work of Nozières and Schmitt-Rink ) that a natural variable which can be used to follow the crossover from strong-coupling to weak-coupling superconductivity is the product $`k_\mathrm{F}\xi `$ of Fermi wave-vector $`k_\mathrm{F}`$ times the coherence length $`\xi `$ for two-electron correlation. In the BE limit, electrons are expected to bind in quasi-bound pairs localized in real space (Schafroth pairs ), thus realizing the condition $`k_\mathrm{F}\xi 1`$, while the BCS limit corresponds to loosely coupled pairs, with $`k_\mathrm{F}\xi 1`$, localized in momentum space close to the Fermi energy.

A zero-temperature calculation of the coherence length for a three-dimensional, $`s`$-wave superconductor along the crossover between the weak- and the strong-coupling limits has been performed both numerically  and analitically . The case of a two-dimensional, $`d`$-wave superconductor has been discussed in Ref. , where a dispersion relation typical of the cuprate superconductors has been explicitly considered. As anticipated above, in the non–$`s`$-wave case ($`\mathrm{}0`$), the standard definition, Eq. (7), of the coherence length leads to unphysical divergences. One may still conveniently define a coherence length as the range in real space of the static correlation function $`X_\mathrm{\Delta }(𝐫)`$ for the modulus of the order parameter  as

$$\xi ^1=\underset{r\mathrm{}}{lim}\frac{\mathrm{log}X_\mathrm{\Delta }(r)}{r}.$$
(23)
The Fourier transform $`X_\mathrm{\Delta }(𝐪)`$ of such a function has been derived for an $`s`$\- and a $`d`$-wave superconductor in Refs.  and , respectively. In the case of our anisotropic interaction, Eq. (19), it reads:

$`X_\mathrm{\Delta }(𝐪)^1`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle \frac{1}{2\mathrm{\Omega }}}{\displaystyle \underset{𝐤}{}}|Y_\mathrm{}m(\widehat{𝐤})|^2{\displaystyle \frac{(u_{𝐤+𝐪/2}u_{𝐤𝐪/2}v_{𝐤+𝐪/2}v_{𝐤𝐪/2})^2}{E_{𝐤+𝐪/2}+E_{𝐤𝐪/2}}}`$ (24)
$`=`$ $`{\displaystyle \frac{1}{2\mathrm{\Omega }}}{\displaystyle \underset{𝐤}{}}|Y|^2\left[{\displaystyle \frac{1}{E}}{\displaystyle \frac{1}{E_++E_{}}}\left(1+{\displaystyle \frac{\xi _+\xi _{}|\mathrm{\Delta }_+||\mathrm{\Delta }_{}|}{E_+E_{}}}\right)\right],`$
where $`EE_𝐤`$, $`YY_\mathrm{}m(\widehat{𝐤})`$, $`\xi _\pm `$, $`\mathrm{\Delta }_\pm `$, $`E_\pm `$ are calculated at momenta $`𝐤\pm 𝐪/2`$, respectively, and use has been made of the gap equation, Eq. (20a), in going into the second line.

First of all, since the summand in Eq. (24) depends only on $`k=|𝐤|`$, $`q=|𝐪|`$, and on the relative angle between $`𝐤`$ and $`𝐪`$, passing to the continuum limit and transforming the sums over wave-vectors into an integral, one has $`X_\mathrm{\Delta }(𝐪)X_\mathrm{\Delta }(q)`$. Moreover, $`X_\mathrm{\Delta }(q)`$ is an even function of its argument. Back to real space, one analogously has $`X_\mathrm{\Delta }(𝐫)X_\mathrm{\Delta }(r)`$. Then, the asymptotic behaviour of $`X_\mathrm{\Delta }(r)`$ at large distance $`r`$ will be governed by the behaviour of its Fourier transform $`X_\mathrm{\Delta }(q)`$ at small wave-vector $`q`$. In particular, assuming the expansion :

$$X_\mathrm{\Delta }^1(q)=a+bq^2+O(q^4),$$
(25)
it is straightfoward to show  that

$$\xi ^2=\frac{b}{a}.$$
(26)
Such a definition of the coherence length is now consistent also for nodal superconductors . It reduces to Eq. (7) in the weak-coupling limit for the $`s`$-wave case, and applies also to the strong-coupling regime, provided $`b>0`$ (given that $`a>0`$, identically, provided $`\mathrm{\Delta }_00`$), *i.e.* provided that $`X_\mathrm{\Delta }(q)`$ has its absolute minimum at $`q=0`$ .

We have numerically solved the gap equations, Eqs. (20), and evaluated the dimensionless coherence length $`k_\mathrm{F}\xi `$ according to this more general definition, Eq. (26), for the anisotropic potential Eq. (19), with $`\mathrm{}=0,1,2`$, corresponding to $`s`$-, $`p`$-, and $`d`$-wave pairing, respectively. Here, a convenient measure of the characteristic energy for superconductivity may be taken as the maximum gap $`\mathrm{\Delta }_{\mathrm{max}}`$, where $`\mathrm{\Delta }_{\mathrm{max}}k_\mathrm{B}T_c`$ holds also for anisotropic superconductors . Such a quantity is readily extracted from the solution of the gap equations, Eqs. (20). In Figure 2, we plot the dimensionless characteristic energy $`2\mathrm{\Delta }_{\mathrm{max}}/\epsilon _\mathrm{F}`$ versus $`\epsilon _c/\epsilon _\mathrm{F}=(k_\mathrm{F}\xi )^2`$ in the cases $`\mathrm{}=0,1,2`$. As $`k_\mathrm{F}\xi `$ decreases \[$`(k_\mathrm{F}\xi )^2`$ increases\], one crosses over from the weak-coupling, BCS limit into the strong-coupling, BE limit . Despite a mean-field solution of the self-consistent gap equations has been now taken into account in the calculation, and a more general and consistent definition of the coherence length has been employed, our numerical results are in good qualitative agreement with the approximate results of Sec. 2, with $`2\mathrm{\Delta }_{\mathrm{max}}/\epsilon _\mathrm{F}`$ increasing as a function of $`\epsilon _c/\epsilon _\mathrm{F}`$ with a steeper slope, as $`\mathrm{}`$ increases from $`\mathrm{}=0`$ ($`s`$-wave) to $`\mathrm{}=2`$ ($`d`$-wave), all curves tending to saturation, at least within the numerically accessible range of the crossover parameter $`k_\mathrm{F}\xi `$. These results are in agreement with the phenomenological plots correlating characteristic energies for both the heavy Fermion and cuprate superconductors in Refs. .

## 4 Summary and future directions

We have studied the Cooper problem for an anisotropic superconductor characterized by an electron-electron interaction expanded in terms of simple spherical harmonics over the Fermi sphere. In the weak-coupling limit of a superconducting instability in the $`(\mathrm{},m)`$ channel, we have derived an analytical expression for the relation between the pair binding energy and the correlation length, for arbitrary relative angular momentum quantum number $`\mathrm{}`$. While such an expression correctly reduces to the standard one in the $`s`$-wave case ($`\mathrm{}=0`$), in view of the fact that $`k_\mathrm{B}T_c`$ scales with such a pair binding energy, in the case of non–$`s`$-wave superconductors ($`\mathrm{}>0`$) our expression agrees qualitatively with the phenomenological correlation between $`k_\mathrm{B}T_c`$ and the characteristic energy $`\mathrm{}^2/(2m^{}\xi ^2)`$ earlier found in Ref. . These results have been confirmed by a numerical solution of the self-consistent gap equations in the crossover between the weak- and the strong-coupling limits, where a more general definition for the coherence length has been employed, applying to anisotropic superconductors.

It may, in the future, prove of significance to include the effect of Coulomb interaction on the formation of Cooper pairs in superconducting assemblies along lines such as those laid down in Ref.  (see also Ref. ). The question arises then as to whether, in the future, it will be of significance in making the present study the basis of fully quantitative calculations, to generalize beyond the result here, and of course, of many earlier workers, that the Cooper pairs formed as a consequence of an attractive interaction correspond to a single value of the binding energy at a given temperature. There is some evidence (see, *e.g.,* Refs. ) that some physical properties, and in particular specific heat and tunneling spectra of cuprate materials, may require generalization to Cooper pairs with a finite range of binding energies. We do not anticipate that the qualitative trends proposed in the present paper will be grossly affected.

G.G.N.A. gratefully acknowledges Dr. N. Andrenacci for useful and stimulating discussions. F.E.L. thanks Professor K. Van Alsenoy for much motivation and encouragement. F.E.L. also acknowledges financial support from the “Program for the Stimulation of Participation in Research Programs of the European Union” of the University of Antwerp. N.H.M. and F.E.L. thank the Department of Physics and Astronomy, University of Catania, for the stimulating environment and much hospitality.


format: text

text:
# THE EFFECTS OF PERIODICALLY GAPPED TIME SERIES ON CROSS-CORRELATION LAG DETERMINATIONS

## 1 Introduction

One of the main advances from recent X-ray observations of blazars is the discovery of the energy-dependent time lags in the X-ray emission of the three bright TeV-emitting blazars Mrk 421, Mrk 501 and PKS 2155–304. The overall spectral energy distributions (SEDs) show that the synchrotron emission component from these sources peaks at high-energy (UV/soft X-ray) band. This indicates that the X-ray emission from these sources is the high energy tail of the synchrotron component, where the most violent variability is expected. The inter-band time lag is one of the important variability parameters. The cross-correlation function (CCF) technique is the standard tool for lag determinations. ASCA and BeppoSAX discovered that in these sources the lower energy X-ray photons may lag or lead the higher energy ones, i.e., the so-called soft or hard lag, respectively. The signs and values of the lags may well depend on either energy or on the single flare analyzed. The typical lags (either soft or hard lag) range from $``$ zero to $`10^4`$ s (see Zhang et al. 2002 for a review), and appear to be correlated to the flare duration: the shorter the flare duration, the smaller the lag (Zhang et al 2002; Brinkmann et al. 2003). Zhang (2002) also found evidence for a dependence of lags on timescales in Mrk 421: lags appear to be larger on longer timescales. The observed lags have been interpreted as evidence for the interplay of the acceleration and cooling timescales of relativistic electrons responsible for the observed X-rays taking place in the jets (e.g., Kirk, Rieger & Mastichiadis 1998), and used to constrain the main parameters of the emitting region on the basis of variability (e.g., Zhang et al. 2002) rather than of SEDs, as commonly used.

However, the latest lag searches with XMM-Newton have put the discovery of ASCA and BeppoSAX into question. Edelson et al. (2001) and Sembay et al. (2002) reported that there is no evidence for measurable inter-band lags in the X-rays, with upper limits of $``$ 0.3 (PKS 2155–304) and 0.08 hr (Mrk 421). It is commonly believed that uninterrupted time series, due to the highly eccentric $``$ 48 hours orbit of XMM-Newton, allow to detect more reliably the lags. This has led these authors to propose that the lags measured with ASCA and BeppoSAX could well be an artifact of the periodic interruptions by the Earth occultation related to the short ASCA and BeppoSAX orbital period ($`1.6`$ hr). Maraschi et al. (2002) also reported no measurable lags with a full orbit XMM-Newton observation of PKS 2155–304. Brinkmann et al. (2003) re-analyzed all the available XMM-Newton observations of Mrk 421, working with light curves that have a small bin-size of 8 s, and found typical lags of $`5`$ min.

In order to address the issue of the reliability and significance of the lags discovered by ASCA and BeppoSAX, in this paper we perform Monte Carlo simulations to investigate the effects of periodically gapped time series on CCF lag determinations. A number of simulations using inverse Fourier transformation, which is model-dependent, have been carried out in the literature to evaluate the significance and uncertainty of the CCF lag determinations either for poorly sampled time series (e.g., Gaskell & Peterson 1987; Maoz & Netzer 1989; White & Peterson 1994) or for evenly sampled time series (e.g., Welsh 1999). The issue of poorly irregular sampling and noisy data was also previously investigated (see the summary in Koen 1994). There are a number of issues that can significantly bias the reliability of CCF lag determinations. We refer to Welsh (1999) for a review on such issue, where the author also showed the effects of power spectral density (PSD) slope, duration, signal-to-noise ratio, de-trending and tapering of light curves. We note, however, that the previous simulations mainly dealt with the relationship between the UV/optical emission line flux and the continuum variations in Seyfert galaxies, which applied to the reverberation mapping (to measure the size of broad-line region) and typical ground-based irregular sampling patterns. In the blazar context, the simulations by Litchfield, Robson, & Hughes (1995) were performed on the radio sampling patterns. Our simulations are tailored to match the sampling characteristics of space-based observations introduced by low earth orbit satellites, which to our knowledge have not been explicitly addressed in astronomical literature. Another distinction is that we discuss the time lags between the variations of the synchrotron emission in different X-ray bands. In particular, it is worth noting that there are obvious difficulties in accurately determining lags in the synchrotron X-ray emission of TeV blazars as (1) the lags are short compared to the bin-sizes of the available data; (2) the data are usually not equally sampled – in particular periodical interruptions of space-based observations are unavoidable for the low Earth orbit satellites; and (3) there are ambiguities in interpreting the complexities of the CCF results.

In this work we also make use of the method introduced by Peterson et al. (1998), namely a model-independent Monte Carlo method to assess the uncertainties in the lag measurements obtained. This method does make use of real data, and is known as flux redistribution/random subset selection (FR/RSS) that is based on the “bootstrap” method. Here we adopt it to study the effects of periodic gaps in lag determinations. To do so, we use a real XMM-Newton observation without interruptions, and re-sample it with the typical BeppoSAX sampling windows. The same CCF analysis is then performed on both the real and fake time series.

In §2 we conduct model-dependent simulations: §2.1 illustrates the assumptions for the simulations to be performed; and the results are presented in §2.2. The model-independent simulations are performed in §3. We discuss the significance of our results in §4.

## 2 Monte Carlo Simulations

### 2.1 Assumptions

The inverse Fourier transformation from frequency to time domain is usually performed to simulate light curves by assuming a PSD model. We use the algorithm of Timmer & König (1995), which randomizes both the amplitude and the phase of a red noise process at each Fourier frequency. In exploring the full variety of possible light curves showing the same PSD, this algorithm is superior to the commonly used one that randomizes phases only (Benlloch et al. 2001; Uttley, McHardy, & Papadakis 2002). In order to mimic a real situation as closely as possible, we normalize fake light curves on the basis of the real variability behavior of one observation of Mrk 421 obtained with BeppoSAX (Fossati et al. 2000; Zhang 2002) as this observation has a good signal-to-noise ratio that will lay stress on the key point of this investigation. The assumptions and procedure of our simulations are as follows.

A simple power-law PSD, i.e., $`P(f)f^\alpha `$ with slope $`\alpha =2.5`$ is assumed (Zhang 2002) to represent the red noise variability of the bright TeV blazars (Kataoka et al. 2001; Zhang et al. 1999; 2002), and from this the light curves are recovered on the basis of Monte Carlo technique. With one set of Gaussian distributed random numbers, we construct a fake pair of light curves that evenly sample 171 points with bin-size $`\mathrm{\Delta }t=512`$ s. This choice resembles the binned light curves of Mrk 421 used in Zhang (2002) to perform the CCF analysis. At the same time, in order to mimic the discovered time lags, we delay the phase of the second light curve of the fake pair on the basis of the relationship $`\mathrm{\Delta }\varphi (f_i)=2\pi f_i\tau (f_i)`$, where $`f_i`$ is the Fourier frequency ($`i=1,\mathrm{},N/2`$; $`N=171`$ is the total number of the evenly sampled light curve points), $`\tau (f_i)`$ is the time lag and $`\mathrm{\Delta }\varphi (f_i)`$ the phase lag at $`f_i`$. For simplicity, we assume that $`\tau `$ is frequency-independent, i.e., $`\tau =`$ constant (see, however, Zhang 2002 for the evidence of the dependence of $`\tau `$ on frequency), and in turn $`\mathrm{\Delta }\varphi (f_i)=2\pi \tau f_i`$. An assumed (true) lag $`\tau `$ therefore corresponds to the lag determined with the CCF methods. The aim of our simulations is to recover the assumed $`\tau `$ by applying the CCF methods to the fake pair with specific sampling windows. The two light curves of the fake pair are then scaled to have the same mean and variance as the real 0.1–2 keV and 2–10 keV light curves of Mrk 421 (Zhang 2002), respectively. In order to mimic photon counting (Poisson) white noise, the two fake light curves are further Gaussian randomly redistributed on the basis of the average errors on the real 0.1–2 keV and 2–10 keV light curves of Mrk 421. In order to simulate the periodic gaps of the light curves obtained with BeppoSAX, we re-sample the two fake light curves by applying the real observing windows on the 0.1–2 keV and 2–10 keV light curves of Mrk 421. Three CCF methods, i.e., interpolation cross correlation function (ICCF, White & Peterson 1994), Discrete Correlation Function (DCF, Edelson & Krolik 1998), and Fisher’s z-transformed DCF (ZDCF, Alexander 1997) that is based on the DCF, are then used to cross-correlate the two fake light curves before and after applying the real observing windows. All CCFs are normalized by the mean and standard deviation of the two cross-correlated light curves using only the data points that actually contribute to the calculation of each lag (White & Peterson 1994), since the light curves in our cases are not stationary. A simulation is deemed to have succeeded if $`r_{\mathrm{max}}`$ (the maximum value of the CCF) between the two fake light curves is significant at a level of confidence greater than $`95\%`$ (Peterson et al. 1998). For each succeeded trial, we record the lags using three techniques to interpret the CCF results: (1) using the lag corresponding to $`r_{\mathrm{max}}`$ of the CCF, $`\tau _{\mathrm{peak}}`$; (2) computing the centroid of the CCF over time lags bracketing $`r_{\mathrm{max}}`$, $`\tau _{\mathrm{cent}}`$ – all the CCF points with $`r`$ in excess of 0.8$`r_{\mathrm{max}}`$ are used (Peterson et al. 1998); (3) fitting the CCF with a Gaussian function to find the location of the CCF peak, $`\tau _{\mathrm{fit}}`$. We repeat this procedure 2000 times to construct probability distributions of the lags (i.e., the cross-correlation peak distribution, CCPD; Maoz & Netzer 1989) for $`\tau _{\mathrm{peak}}`$, $`\tau _{\mathrm{cent}}`$, and $`\tau _{\mathrm{fit}}`$.

### 2.2 Results

We simulate time lags covering the range discovered by ASCA and BeppoSAX in TeV blazars. Therefore, we arbitrarily assume $`\tau =0,300,1400,3000,5400,7100`$ s, respectively. We create fake light curves using the same sets of random number for each case except for different $`\tau `$.

Figure 1a, as an example, shows a pair of fake light curves that are evenly sampled, and the second light curve (solid circles) is delayed by $`\tau =3000`$ s with respect to the first one (open circles). The lag is clearly visible by comparing both the peaks and the troughs of the two light curves. In this case, the ratio of $`\tau `$ to $`\mathrm{\Delta }t`$ of the light curves is $`6`$. The ICCF, DCF and ZDCF of the two light curves are shown in Figure 1b-d (solid lines). We use $`\tau _{\mathrm{peak}}`$, $`\tau _{\mathrm{cent}}`$, and $`\tau _{\mathrm{fit}}`$ to measure the lag, respectively. The results are reported in Table 1. One can see that the true lag is properly recovered by the different CCF methods and different techniques used to interpret the CCFs.

We then re-sample these two light curves using the real BeppoSAX observing windows as mentioned in §2.1. The resulted light curves thus are affected by periodic gaps resembling the light curves obtained with BeppoSAX. Note that the number of points (63) in the first re-sampled light curve is smaller than that in the second one (76) because on-board BeppoSAX the LECS detector is less exposed than the MECS one. The ICCF, DCF and ZDCF of the two light curves with periodic interruptions are shown in Figure 1b-d (dotted line or open circles with error bars). The measured lags are also reported in Table 1. In this case, both $`\tau _{\mathrm{peak}}`$ and $`\tau _{\mathrm{cent}}`$ underestimate the true lag while $`\tau _{\mathrm{fit}}`$ of the DCF and ZDCF still recover it, although with larger uncertainties. Figure 1b-d also show that the overall characteristics of the CCFs obtained from the evenly sampled light curves are almost identical to those obtained from the corresponding light curves with periodic gaps. However, due to the missing of a large number of data points in the latter case, the CCF peaks somewhat shift to smaller lags. Note that $`\tau _{\mathrm{fit}}`$ is not calculated for the ICCF case because there are no estimates of the ICCF errors, and the DCF errors are overestimated with respect to the ZDCF ones.

The statistical significance of the effects of periodic gaps on the CCF lag determinations is deduced from the probability distributions of a number of simulations. We show in Figure 2 (ICCF), Figure 3 (DCF), and Figure 4 (ZDCF) the CCPDs of the simulated lags for the case $`\tau =0`$ and $`3000`$ s. It is worth noting that the CCPDs are almost always non-normal distributions, in particular for $`\tau _{\mathrm{peak}}`$ and $`\tau _{\mathrm{cent}}`$. Therefore, we use the median and $`68\%`$ confidence level (with respect to the median) to statistically characterize the CCPDs. The statistical results for all the simulations performed are tabulated in Table 2. The first column is the assumed true lag, and columns 2–4 and 5–7 give the CCPD median and $`68\%`$ confidence range for evenly sampled and periodically gapped light curves, respectively. As shown in Figures 2-4 and Table 2, the main results of our simulations can be summarized as follows: (1) the assumed true lags are recovered in all cases in terms of the CCPD medians with $`68\%`$ confidence errors; (2) the main effect of periodic gaps is to broaden the CCPDs, thus to increase the uncertainty of the CCF lag determinations. Most importantly, periodic gaps do not produce artificial CCF lags. More specifically, from the simulations with true small lags (e.g, the cases of 0 and 300 s lags), one can see that periodic interruptions in the light curves definitely do not produce spurious lags on timescales of hours; (3) in some cases, the three CCF methods and the three techniques used to interpret the CCF results do not give rise to completely consistent results. For example, $`\tau _{\mathrm{peak}}`$ strongly depends on the lag steps used to calculate the CCF; the errors of the DCF are overestimated, thus producing the broadest CCPDs.

These results imply that (1) evenly sampled light curves with large lag to bin-size ratios give more reliable CCF lag determinations; (2) more importantly, the lags on timescales of hours discovered by ASCA and BeppoSAX can not be an artifact of periodic gaps in the light curves that have intrinsically small lag; (3) the only effect of periodic gaps is to introduce uncertainty and to increase the variance on the CCF lag determinations.

## 3 A Specific Case: XMM-Newton Observations of PKS 2155–304

The simulations presented in the previous section showed that periodical gaps in the light curves do not produce spurious lags between them, but only increase the uncertainty and variance in lag determinations. Uninterrupted data with high temporal resolution are available from XMM-Newton observations of Mrk 421 and PKS 2155–304. All of the CCF analysis, performed by Brinkmann et al. (2001; 2003), Edelson et al. (2001), Sembay et al. (2002), and Maraschi et al. (2002) showed that the inter-band lags between the soft and hard energy band are close to zero, with upper limits of about 1000 s. As we pointed out in the Introduction, these results led Edelson et al. (2001) to suggest that previous claims of time lags on time scales of hours might be an artifact of the periodic interruptions every $``$1.6 hours due to the low-Earth orbits of satellites such as ASCA and BeppoSAX. Therefore, in addition to the simulations presented in the previous section arguing against the above suggestion, an important test to assess the role of periodic gaps is to consider the XMM-Newton data and re-sample them according to the ASCA or BeppoSAX observing windows, and then repeat the CCF analysis on the fake ASCA or BeppoSAX data. Whether or not inter-band lags would be detected from the fake data would be a strong argument in favor or against the claim by Edelson et al. (2001).

In order to perform such a test, we take the first part of the XMM-Newton observation of PKS 2155–304 (Maraschi et al. 2002). The details of the data reduction and CCF analysis will be presented in Maraschi et al. (in preparation). For our purposes we just extracted the light curves in two energy bands, i.e., 0.7-1 and 1-2 keV. The light curves (without interruptions) are shown in Figure 5a. First, we performed the CCF analysis on them: the results are shown in Figure 5b-d (solid line) for the ICCF, DCF and ZDCF, respectively. We then re-sampled the light curves with the typical BeppoSAX sampling windows, obtaining two light curves with periodic interruptions resembling real BeppoSAX observations. We performed the same CCF analysis on the fake light curves with periodic gaps just as we did on light curves without gaps. The results are also shown in Figure 5b-d (dotted line or open circles with error bars). In all cases the two light curves are highly correlated near zero lag, and the CCFs calculated from the real and fake data overlap near their peaks, suggesting that the periodic interruptions do not change the CCF character near the peak. We also measured the lags with the three techniques. The results are tabulated in Table 3. Due to the complexities of the CCFs, the measured lags depend on the CCF methods and the techniques used to quantify the CCFs. However, the results show in general that the real lag may be close to zero, with 0.7-1 keV photons very marginally leading the 1-2 keV photons.

Finally, we performed FR simulations (i.e., Gaussian randomly redistributing the light curves on the basis of the quoted errors). The resulting CCPDs are shown in Figure 6. We used the same statistical method as we did in § 2.2 to characterize the CCPDs. The statistical results are tabulated in Table 4. Within the $`68\%`$ confidence errors, $`\tau _{\mathrm{peak}}`$ gives lags consistent with zero, while $`\tau _{\mathrm{cent}}`$ and $`\tau _{\mathrm{fit}}`$ suggest positive lags of $``$1000–2000 s, but at a low confidence ($``$ 2$`\sigma `$). Note that the low confidence of the lag detections may be caused by the complexities of the variability behavior and the dependence of variations on the energy band considered. This also explains the fact that different methods and techniques may give rise to inconsistent results. In any case the comparison of the results obtained – using the same method and technique – from the light curves with periodic gaps and without gaps does not favor the suggestion that the periodic interruptions may produce spurious lags.

## 4 Discussion and Conclusions

We performed two sets of simulations to investigate the effects of space-based observations with short orbital period ($`1.6`$ hr) on the reliability and significance of CCF lag determinations, specifically for the energy-dependent variations of synchrotron X-ray emission of TeV blazars. The first set of simulations (§2) make use of fake light curves generated with the Fourier transformation method. We investigated two main issues: (1) the effects of periodic data gaps in the light curves; and (2) the effects of different lag to bin-size ratios of light curves. The simulations showed that evenly sampled light curves indeed yield more reliable CCF lags, with smaller variance, than the periodically gapped light curves do. However, the CCPD analysis clearly showed that the light curves with periodic gaps still preserve the nature of the true lags (regardless of the values) even though they introduce larger lag variances than the light curves without gaps. Moreover, larger ratio of lag to light curve bin-size can improve the significance of CCF lag determination. The second set of simulations is based on a real XMM-Newton observation without interruptions. We re-sampled it with the typical BeppoSAX sampling windows in order to study whether or not periodic interruptions may give rise to strong biases in CCF lag determinations. The complex nature of the variability in TeV blazars results in obvious difficulties when quantifying lags, and it is likely to be the source of discrepancies between the results quantified with different techniques. However, the comparison of the results derived with the same CCF method and the same quantifying technique, shows that the light curves with periodic interruptions do not produce spurious lags on timescales of hours if their intrinsic lag is indeed small. Therefore, our investigations argue against the proposal by Edelson et al. (2001) that the lags of about an hour discovered by ASCA and BeppoSAX are an artifact of periodic gaps introduced by low Earth orbit satellites. We thus conclude that the lags discovered by ASCA and BeppoSAX are most likely due to intrinsic variability properties of the sources, and not artificially produced from an intrinsic zero lag by periodic gaps. However, due to the complexities of variability that produces complicated CCF – irregularities and complexities of CCF could well be caused by different properties of variability on different timescales (Zhang 2002) – it is not easy to quantify real lags with the CCF method and different CCF methods and interpreting techniques may work for different cases. Our simulations also confirm that uninterrupted light curves with large lag to bin-size ratios can improve accuracy of lag determinations, in particular for small lags that require high sampling rates of light curves.

Welsh (1999) showed that the reality of CCF lag determinations also depend (1) on the light curve auto-correlation function (ACF) sharpness (the sharper the ACFs, the narrower the CCF peak and the smaller the lag bias) and (2) on the ratio of the intrinsic lag to the duration of the light curves. The first dependence can be easily explained by recalling that a CCF is a convolution of two ACFs. The ACF sharpness is determined by the PSD steepness (the steeper the PSD, the broader the ACF). The X-ray PSDs of the three TeV blazars are steep (with slopes of $``$ 2–3; Kataoka et al. 2001; Zhang et al. 1999; 2002), the ACF and the CCF peaks are therefore broad. De-trending light curves might remove such bias, but it also removes the low-(Fourier) frequency variability of the sources. This can introduce serious errors into time series and needs to be done very carefully. The second dependence becomes important only for the lengths of light curves shorter than $`4`$ times the lag, which is never the case for the TeV blazars observed in the X-rays (typical ratios of the lengths of light curves to lags are 10–100).

Finally, we stress that the real light curves of the bright TeV blazars are very complex (e.g., the relationship between light curves at different energies may not be represented by just one “fixed” lag). Such complexities definitely result in irregular CCF, e.g., the CCF peaks at zero lag but shows asymmetry, which makes the CCF methods less straightforward to lag determinations. The cross-spectral technique, a more complex tool used to determine the Fourier frequency-dependent lags, may have the advantage of avoiding such ambiguities at least over long timescales (see Zhang 2002 for details). However, because this method relies on the Fourier transformation, it is not applicable to unevenly spaced data. On the contrary, in time domain one can use the DCF to substitute the classical CCF when dealing with irregular data.

We thank the anonymous referee for constructive comments. This research has made use of the standard on-line archive provided by the BeppoSAX Science Data Center (SDC). This work is partly based on an observation with XMM-Newton, an ESA science mission. The Italian MIUR is acknowledged for financial support. The project is sponsored by the Scientific Research Foundation for the Returned Overseas Chinese Scholars, State Education Ministry.


text:
# Measurement of the Negative Muon Anomalous Magnetic Moment to 0.7 ppm

## Abstract

The anomalous magnetic moment of the negative muon has been measured to a precision of 0.7 parts per million (ppm) at the Brookhaven Alternating Gradient Synchrotron. This result is based on data collected in 2001, and is over an order of magnitude more precise than the previous measurement of the negative muon. The result $`a_\mu ^{}=\mathrm{11\hspace{0.17em}659\hspace{0.17em}214}(8)(3)\times 10^{10}`$ (0.7 ppm), where the first uncertainty is statistical and the second is systematic, is consistent with previous measurements of the anomaly for the positive and negative muon. The average for the muon anomaly is $`a_\mu (\text{exp})=\mathrm{11\hspace{0.17em}659\hspace{0.17em}208}(6)\times 10^{10}`$ (0.5 ppm).

The anomalous magnetic moments of the muon and the electron have played an important role in the development of the standard model. Compared to the electron, the muon anomaly has a relative sensitivity to heavier mass scales which typically is proportional to $`(m_\mu /m_e)^2`$. At the present level of accuracy the muon anomaly gives an experimental sensitivity to virtual $`W`$ and $`Z`$ gauge bosons as well as a potential sensitivity to other, as yet unobserved, particles in the few hundred GeV/$`c^2`$ mass range kh .

We report our result for the negative muon anomalous magnetic moment $`a_\mu ^{}=(g2)/2`$ from data collected in early 2001. The measurement is based on muon spin precession in a magnetic storage ring with electrostatic focusing. The same experimental technique was used as in our most recent measurements of $`a_{\mu ^+}`$ muplus ; g2\_2000 , and a similar precision of 0.7 ppm was achieved. Detailed descriptions of the apparatus may be found elsewhere nimpapers ; fei ; nimpapers2 ; kicker ; quads .

For polarized muons moving in a uniform magnetic field $`\stackrel{}{B}`$ perpendicular to the muon spin and to the plane of the orbit and in an electric quadrupole field $`\stackrel{}{E}`$, which is used for vertical focusing quads , the angular frequency difference, $`\omega _a`$ between the spin precession frequency and the cyclotron frequency, is given by

$$\stackrel{}{\omega }_a=\frac{e}{mc}\left[a_\mu \stackrel{}{B}\left(a_\mu \frac{1}{\gamma ^21}\right)\stackrel{}{\beta }\times \stackrel{}{E}\right].$$
(1)
The dependence of $`\omega _a`$ on the electric field is eliminated by storing muons with the “magic” $`\gamma =29.3`$ cern , which corresponds to a muon momentum $`p=3.09`$ GeV/$`c`$. Hence measurement of $`\omega _a`$ and of $`B`$, in terms of the free proton NMR frequency $`\omega _p`$ and the ratio of muon to proton magnetic moments $`\lambda `$, determines $`a_\mu `$. At the magic $`\gamma `$, the muon lifetime is approximately $`64.4\mu \mathrm{s}`$ and the $`(g2)`$ precession period is $`4.37\mu \mathrm{s}`$. With a field of 1.45 T in our storage ring nimpapers , the central orbit radius is 7.11 m.

The difference frequency $`\omega _a`$ was determined by counting the number $`N(t)`$ of decay electrons above an energy threshold. The time spectrum of decay electrons is then given by

$$N(t)=N_0e^{t/(\gamma \tau )}\left[1+A\mathrm{sin}(\omega _at+\varphi _a)\right].$$
(2)
The normalization $`N_0`$, asymmetry $`A`$, and phase $`\varphi _a`$ vary with the chosen energy threshold.

The measurement of the magnetic field frequency $`\omega _p`$ is based on proton NMR in water. A trolley with 17 NMR probes was moved typically every three days throughout the entire muon storage region. About 150 fixed NMR probes distributed around the ring in the top and bottom walls of the vacuum chamber were used to interpolate the field between trolley measurements. The system was calibrated with respect to a standard probe with a spherical $`\text{H}_2\text{O}`$ sample fei . The homogeneity of the field in 2001 (Figure 1) was similar to that achieved for the opposite polarity field in 2000 g2\_2000 .

The field $`B`$ weighted with the analyzed event sample was obtained from two largely independent analyses, whose results were found to agree to within 0.05 ppm. Its final value is expressed in terms of the free proton resonance frequency and is given by $`\omega _p/(2\pi )=\mathrm{61\hspace{0.17em}791\hspace{0.17em}400}(11)\text{Hz}`$ (0.2 ppm). Table 1 lists the uncertainties. The improved 2001 uncertainties resulted from refinements in the calibration measurements, and from an upgraded system to determine the azimuthal trolley position in the storage ring.

The 2001 $`\omega _a`$ data taking was similar to that in 2000. However, the hardware energy threshold of the detectors was kept lower and equal for all counters at 0.9 GeV compared to 1.0-1.4 GeV in 2000. This was made possible by reducing the intensity of the injected beam, which in turn reduced the light flash in the detectors muplus ; g2\_2000 . These factors allowed all the detectors to be turned on and be stable by 32 $`\mu `$s after beam injection, as opposed to 50 $`\mu `$s in 2000. As a result of the reduced rates, the fraction of overlapping signals (pileup) after 32 $`\mu `$s in 2001 was comparable to the pileup fraction after 50 $`\mu `$s in 2000. In 2000 the field focusing index $`n`$, which is proportional to the electric field gradient, was $`n=0.137`$, corresponding to a horizontal coherent betatron oscillation frequency (CBO) of 466 kHz g2\_2000 . This frequency was close to twice the $`(g2)`$ frequency of 229 kHz, which resulted in a sizable uncertainty in the fitted $`\omega _a`$ value g2\_2000 . In 2001 we used two different $`n`$-values, $`n=0.122`$ and $`n=0.142`$, which resulted in CBO frequencies, 419 kHz and 491 kHz that are further from twice the $`(g2)`$ frequency (see Figure 2). Consequently, the uncertainty caused by CBO is smaller. Furthermore it also reduced the correlation between the CBO and detector gain effects in the fits to the time spectrum.

Two independent implementations of the algorithm to reconstruct the electron times and energies from the calorimeter signals were used. The frequency $`\omega _a`$ was determined by fitting the time distribution of decay electrons. Five independent analyses were performed in order to probe the systematic uncertainties and, of course, to protect against mistakes. All five results agreed within the expected statistical deviations due to different data selection and weightings. These analyses are described below.

Two of the analyses used slightly different parametrizations muplus ; g2\_2000  that included CBO modulations and fitted the combined electron spectrum in the energy range 1.8-3.4 GeV. In the third analysis, the counts were weighted with the experimentally determined energy-dependent modulation asymmetry, which optimized the statistical power of the data. This method permitted the analyzed energy range to be extended. We used an energy range of 1.5 to 3.4 GeV, which together with the asymmetry weighting resulted in a 10% improvement of the statistical uncertainty. As in the first two analyses, the resulting spectrum of weighted counts was fitted to a function that parametrized all known and statistically significant perturbations.

The remaining analyses fit the ratio muplus ; g2\_2000  formed by randomly assigning the data to four statistically independent subsets $`n_1`$ to $`n_4`$. The subsets were rejoined in $`u(t)=n_1(t)+n_2(t)`$ and $`v(t)=n_3(t\tau _a/2)+n_4(t+\tau _a/2)`$, where $`\tau _a`$ is an estimate of the $`(g2)`$ period, and then combined to form the time spectrum $`r(t)=[u(t)v(t)]/[u(t)+v(t)]`$. The $`(g2)`$ rate modulation of $`v`$ is 180 degrees out of phase compared to that of $`u`$, and to sufficient precision $`r(t)`$ can be described by $`A\mathrm{sin}(\omega _at+\varphi _a)`$. The ratio $`r(t)`$ is largely insensitive to changes of observed counts on time scales larger than $`\tau _a=2\pi /\omega _a4\mu \mathrm{s}`$.

In one of the ratio analyses, the sensitivity to CBO was reduced by combining the data from both $`n`$-values and all detectors prior to fitting. The data were fitted from 32 $`\mu `$s after injection when all detectors were on. In the second ratio analysis the data were fitted separately for each calorimeter and $`n`$value. The fits began between 24 $`\mu `$s and 32 $`\mu `$s, and required the parametrization of the CBO effects in the fit function.

Changes in the radial and vertical muon distributions with time were quantified, and were found to have negligible effect on $`\omega _a`$. A small reduction in the pulsed electrostatic quadrupole voltages quads  during the measurement period could change the vertical muon distribution. Analysis of the data from scintillator counter hodoscopes placed in front of the calorimeters combined with a beam tracking calculation and a GEANT based simulation set a systematic error limit of 0.03 ppm. The muon radial distribution is determined by the magnetic field and the momentum distribution muplus ; yuri . The magnetic field does not change with time after injection, except due to the field from eddy currents induced by the fast kicker kicker . This was measured, and found to have a negligible effect on the muon radial distribution. Muons of lower momenta decay earlier in the laboratory frame than muons of higher momenta. The momentum distribution of the stored beam thus changes during the $`600\mu \mathrm{s}`$ measurement period. The effect on $`\omega _a`$ due to this change was studied in simulation, and was found to be 0.03 ppm.

The results for $`\omega _a`$ for the two $`n`$-values are consistent, see Figure 3, and were combined for each of the analyses. The values for $`\omega _a`$ from the five analyses are in agreement to within variations expected from the differences in the analyzed event samples and the treatment of the data. The analysis techniques are expected to have somewhat different sensitivities to different systematic effects. Detailed comparisons of the results, using all analyzed data as well as only the data in overlap, showed no evidence for unaccounted systematic differences. The five resulting values for $`\omega _a`$ were combined in a simple arithmetic mean to obtain a single value for $`\omega _a`$.

The resulting frequency value is $`\omega _a/(2\pi )=\mathrm{229\hspace{0.17em}073.59}(15)(5)\mathrm{Hz}`$ (0.7 ppm), which includes a correction of $`+0.77(6)`$ ppm for contributions to Eq. 1 caused by vertical oscillations (0.30 ppm) and for the effect of the horizontal electric fields on muons with $`\gamma 29.3`$ (0.47 ppm). The stated uncertainties account for strong correlations among the individual results, both statistical and systematic. Table 2 lists the systematic uncertainties in the combined result with these correlations taken into account.

After the $`\omega _p`$ and $`\omega _a`$ analyses were finalized separately and independently, $`a_\mu `$ was evaluated. The result is

$$a_\mu ^{}=\frac{R}{\lambda R}=\mathrm{11\hspace{0.17em}659\hspace{0.17em}214}(8)(3)\times 10^{10}\text{(0.7\hspace{0.17em}ppm)},$$
(3)
where $`R_\mu ^{}\omega _a/\omega _p=\mathrm{0.003\hspace{0.17em}707\hspace{0.17em}208\hspace{0.17em}3}(\mathrm{2\hspace{0.17em}6})`$ and $`\lambda =\mu _\mu /\mu _p=\mathrm{3.183\hspace{0.17em}345\hspace{0.17em}39}(10)`$ liu . This new result is in good agreement with the average of $`R_{\mu ^+}=\mathrm{0.003\hspace{0.17em}707\hspace{0.17em}204\hspace{0.17em}8}(\mathrm{2\hspace{0.17em}5})`$ muplus  as predicted by the CPT theorem. The difference $`\mathrm{\Delta }R=R_\mu ^{}R_{\mu ^+}=(3.5\pm 3.4)\times 10^9`$. The new average is $`R_\mu =\mathrm{0.003\hspace{0.17em}707\hspace{0.17em}206\hspace{0.17em}3}(\mathrm{2\hspace{0.17em}0})`$ and

$$a_\mu (\mathrm{exp})=\mathrm{11\hspace{0.17em}659\hspace{0.17em}208}(6)\times 10^{10}\text{(0.5\hspace{0.17em}ppm)},$$
(4)
in which the total uncertainty consists of $`5\times 10^{10}`$ (0.4 ppm) statistical uncertainty and $`4\times 10^{10}`$ (0.3 ppm) systematic uncertainty. The correlation of systematic uncertainties between the data sets has been taken into account. The combined result for the positive muon g2\_2000 , $`a_{\mu ^+}(\mathrm{exp})=\mathrm{11\hspace{0.17em}659\hspace{0.17em}203}(8)\times 10^{10}`$ (0.7 ppm) has a statistical uncertainty of $`6\times 10^{10}`$ (0.6 ppm) and a systematic uncertainty of $`5\times 10^{10}`$ (0.4 ppm). It is shown in Figure 4 together with the new result for the negative muon and their average.

The standard model prediction for $`a_\mu `$ consists of QED, hadronic and weak contributions. The uncertainty on the standard model value is dominated by the uncertainty on the lowest-order hadronic vacuum polarization. This contribution can be determined directly from the annihilation of $`e^+e^{}`$ to hadrons through a dispersion integral disp . The indirect determination using data from hadronic $`\tau `$ decays, the conserved vector current hypothesis, plus the appropriate isospin corrections, could in principle improve the precision of $`a_\mu (\mathrm{had})`$. However, discrepancies between the $`\tau `$ and the $`e^+e^{}`$ results exist dehz03 ; jeg . The two data sets do not give consistent results for the pion form factor. Using $`e^+e^{}`$ annihilation data the corresponding theoretical value is $`a_\mu (\mathrm{SM})=\mathrm{11\hspace{0.17em}659\hspace{0.17em}181}(8)\times \mathrm{\hspace{0.17em}10}^{10}(0.7\mathrm{ppm})`$. The value deduced from $`\tau `$ decay is larger by $`15\times 10^{10}`$ and has a stated uncertainty of $`7\times \mathrm{\hspace{0.17em}10}^{10}(0.6\mathrm{ppm})`$. The difference between the experimental determination of $`a_\mu `$ and the standard model theory using the $`e^+e^{}`$ or $`\tau `$ data for the calculation of the hadronic vacuum polarization is 2.7 $`\sigma `$ and 1.4 $`\sigma `$, respectively.

This is the final analysis of the anomalous magnetic moment from experiment E821 at the Brookhaven Alternating Gradient Synchrotron. We aim to substantially improve our result in a new measurement and look forward to continued efforts to improve the theoretical evaluation.

We thank T. Kirk, D.I. Lowenstein, P. Pile, and the staff of the BNL AGS for the strong support they have given this experiment. This work was supported in part by the U.S. Department of Energy, the U.S. National Science Foundation, the U.S. National Computational Science Alliance, the German Bundesminister für Bildung und Forschung, the Russian Ministry of Science, and the U.S.-Japan Agreement in High Energy Physics.


text:
# Error-correcting codes on scale-free networks

## Abstract

We investigate the potential of scale-free networks as error-correcting codes. We find that irregular low-density parity-check codes with highest performance known to date have degree distributions well fitted by a power-law function $`p(k)k^\gamma `$ with $`\gamma `$ close to 2, which suggests that codes built on scale-free networks with appropriate power exponents can be good error-correcting codes, with performance possibly approaching the Shannon limit. We demonstrate for an erasure channel that codes with power-law degree distribution of the form $`p(k)=C(k+\alpha )^\gamma `$, with $`k2`$ and suitable selection of the parameters $`\alpha `$ and $`\gamma `$, indeed have very good error-correction capabilities.

A variety of complex networks Alb1  exhibit a topological structure in which the connectivity between their constituent nodes follows a simple power law. Examples of such scale-free networks include the Internet Fal ; Cal , the World Wide Web (WWW) Alb2 ; Hub , social networks New1 , metabolic networks Jeong , etc. Extensive studies have been made to understand the topological features and evolving dynamics Bara ; Krap ; Doro  of these networks. While many intriguing properties concerning the structural aspect of complex networks have been revealed thanks to these efforts, relatively little has been known about the effects of specific connectivity structures on networks’ functional behavior Stro . In order to properly operate under a certain environment or, in a more active sense, to successfully accomplish a given task, a complex network may favor one particular structure over another. For practical applications, it now appears that more attention need be paid to the functional aspect of these complex networks viewed as whole systems or organisms working for particular purposes.

Recent advances in channel coding theory have led to the perception that the state-of-the-art capacity-approaching codes, such as Turbo codes Berr  and low-density parity-check (LDPC) codes Gall ; Mac1 ; Rich , can be understood in terms of graphs (or networks) consisting of nodes and edges Forney . The function of these graphs is to carry out error correction, i.e. to recover original data transmitted over noisy channels, by iteratively passing certain messages through edges connecting neighboring nodes. The art of developing a high-performance error-correcting code lies in designing a connectivity structure of a graph in such a way as to make the code built on it perform a desired function. One very important issue concerns finding the connectivity distribution that achieves the Shannon’s capacity limit. Most attempts, however, have been limited to numerical optimization, and a complete understanding of the connectivity structure specific to capacity-achieving codes is still lacking. Inspired by the ubiquitous nature of scale-free networks, one may ask whether the connectivity structure of scale-free networks could offer any insight into seeking good graph-based codes.

In this paper we address the question whether scale-free networks whose connectivity distribution follows a power law can function effectively as good error-correcting codes. The codes built on scale-free networks considered here are basically LDPC codes in that the associated parity-check matrices are sparse and that the belief propagation algorithm Gall ; Mac1 ; Rich  is employed for decoding. We first show that the degree distributions of LDPC codes with highest performance known to date are well fitted by power-law functions. Motivated by this finding, we generate a degree distribution according to the function $`p(k)=C(k+\alpha )^\gamma `$ and fine tune the parameters $`\alpha `$ and $`\gamma `$ to maximize the code’s performance. We investigate error-correction capability of these codes over a binary erasure channel and compare them with the Tornado code Luby , the first commercialized LDPC code.

An LDPC code can be represented by a bipartite graph in which there are two different types of nodes: variable nodes and check nodes. Nodes of one type are connected by edges only to nodes of the other type. Variable nodes are associated with data bits, and check nodes examine whether the variable nodes connected to them satisfy parity-check equations. Error correction of corrupted data bits is performed by passing certain messages, e.g. likelihood ratios, through edges back and forth between variable and check nodes. It is known from the density evolution analysis Rich  that, under the assumption of a tree-structured random graph with no closed loops, the error-correction capability of a code is solely determined by the degree distribution.

We begin by inspecting the degree distributions of some high-performance LDPC codes. Figure 1(a) shows the variable-node degree distribution of the LDPC code designed by Chung et al. Chung , which has been optimized for an additive white Gaussian noise channel and approaches the Shannon limit within 0.0045 dB, presently the world record. Here, in order to obtain a meaningful distribution from the irregularly spaced data $`\lambda (k_i)`$ in Table 2 of Chung et al. Chung , we took a local average over a bin of length $`(k_{i+1}k_{i1})/2`$:

$$p(k_i)=\frac{P(k_i)}{(k_{i+1}k_{i1})/2},$$
(1)
where $`P(k_i)`$ is the fraction of nodes with degree $`k_i`$ and is given by $`P(k_i)=C\lambda (k_i)/k_i`$, in which $`\lambda (k_i)`$ is the fraction of edges connected to a variable node of degree $`k_i`$ and $`C`$ is a normalization constant. It can be seen from Fig. 1(a) that the degree distribution is well fitted by a power-law function $`p(k)k^\gamma `$ with $`\gamma 2.14`$. A more dramatic correspondence is observed for the variable-node degree distribution of the Tornado code Luby ; Shok , as apparently seen in Fig. 1(b). The Tornado code has been optimized for an erasure channel and has a Poisson distribution for its check-node degree distribution. The power-law function that best fits the variable-node degree distribution is found to have an exponent $`\gamma 2.02`$, and the fitting appears to be nearly perfect for large $`k`$. We also find that the right-regular sequence of Shokrollahi Shok  that slightly beats the Tornado sequence is well fitted by a similar power-law function.

The fact that many high-quality LDPC codes have degree distributions well fitted by power-law functions stimulates us to test scale-free networks as error-correcting codes. Following the approach of Newman et al. New2 , we write the generating function for a general scale-free network in the form:

$$G(x)=\underset{k=0}{\overset{d_l1}{}}a_kx^k+\underset{k=d_l}{\overset{d_{max}}{}}p(k)x^k,$$
(2)
where the fraction of nodes with degree $`kd_l`$ follows a power law $`p(k)=Ck^\gamma `$, and the terms $`a_k`$ with low degree $`k<d_l`$ are separated from the second sum to allow for possible deviation from the power law for small $`k`$, as is often the case for general scale-free networks. The generating function in Eq. (2), however, contains too many parameters to be amenable to numerical optimization unless $`d_l`$ is sufficiently small. To reduce the number of parameters while still retaining the possibility that the distribution for low degrees may not obey an exact power law, we instead choose the following generating function:

$$G(x)=\underset{k=0}{\overset{d_{max}}{}}p(k)x^k,$$
(3)
where $`p(k)=C(k+\alpha )^\gamma `$. If $`\alpha <0`$ ($`\alpha >0`$), the degree distribution for small $`k`$ lies above (below) the power-law function $`k^\gamma `$. We henceforth use Eq. (3) to generate a variable-node degree distribution of our code and optimize the parameters $`\alpha `$ and $`\gamma `$ to achieve the best performance.

Some empirical results known about LDPC codes help us to further refine our code. The most well known findings related to features of good LDPC codes may be that the variable nodes of degree one should be removed since they do not contribute to error correction and that the codes with almost uniform check-node degree yield good performance Rich ; Chung ; Mac2 . Taking these into account, we let the sum in Eq. (3) start from $`k=2`$, and restrict the check-node degree to two consecutive integers: the generating function for the check-node degree is written as $`F(x)=bx^i+(1b)x^{i+1}`$, where the parameters $`b`$ and $`i`$ are easily determined once a variable-node degree distribution is selected. This choice of the check-node degree distribution enables us to design a code without restrictions on $`d_{max}`$ for any given code rate; this property, however, is not shared by the right-regular sequence Shok  for which $`d_{max}`$ is allowed to have only a special set of values.

The performance of an LDPC code over a binary erasure channel can be evaluated by the density evolution method Rich  as follows. Let $`\delta `$ be the erasure probability of a given channel, and consider a code with a degree distribution pair $`\lambda (x)=\lambda _kx^{k1}`$ and $`\rho (x)=\rho _kx^{k1}`$, where $`\lambda _k`$ ($`\rho _k`$) is the fraction of edges connected to a variable (check) node of degree $`k`$. Note that the distribution here is defined in terms of the fraction of edges, not the fraction of nodes as before. Then, if the belief propagation algorithm is used for decoding, the messages passed between the variable and check nodes are known to evolve as Rich ; Luby

$$x_l=x_0\lambda (1\rho (1x_{l1})),$$
(4)
where $`x_l`$ denotes the expected fraction of erasure messages at the $`l`$th iteration and $`x_0`$ is its initial value given by $`x_0=\delta `$. The recovery of original data is successfully done if $`x_l`$ converges to zero. The threshold $`\delta ^{}`$, defined by the supremum of all $`\delta `$ that result in successful decoding, tells the code’s performance. For a given code rate $`R`$, the threshold is upper bounded by $`1R`$ Luby .

With the help of the above density evolution method, we calculate the error-correction capability of the scale-free networks given in the form of Eq. (3). The results are shown in Fig. 2 as a function of the maximum variable-node degree $`d_{max}`$, where the code rate is fixed at $`R=0.5`$. It is seen that the threshold erasure probability $`\delta ^{}`$ increases as the maximum variable-node degree increases. For large $`d_{max}`$, the threshold almost reaches the theoretical upper bound $`1R`$, indicating that the error-correction capability of our code is very good. For comparison, we have also studied the error-correction capability of codes that have degree distributions other than the power-law distribution, namely the exponential distribution of the form $`p(k)e^{\beta (k+\alpha )}`$ and the Gaussian distribution of the form $`p(k)e^{\beta (k+\alpha )^2}`$. We find that the threshold for the exponential distribution rapidly increases with $`d_{max}`$ and converges to 0.465, a value much lower than the threshold for the power-law distribution. The case for the Gaussian distribution is observed to exhibit a similar behavior with a similar, low convergence limit.

To more clearly demonstrate the high performance of codes on scale-free networks, we compare them with the Tornado code Luby ; Shok . The threshold of our code is presented in Table 1 along with the parameters $`\alpha `$ and $`\gamma `$ that maximize the code’s performance. Table 1 shows that our code yields better performance than the Tornado code for $`d_{max}`$ smaller than about 1000. Also shown in Table 1 is the average variable-node degree $`k`$ of the two codes. From a practical viewpoint, it is important to design a code that yields good performance for small $`k`$, since the physical complexity of a code, which grows with increasing $`k`$, limits the hardware implementation of the code. For this reason, our code seems to be better suited to applications than the Tornado code.

Another merit of our code is that the iteration number required for convergence of decoding is very small. The iteration numbers of our code and the Tornado code are compared in Fig. 3(a) as a function of the erasure probability for the case of $`d_{max}=610`$, which clearly shows that our code has a smaller iteration number than the Tornado code in the whole region of $`\delta `$. Even for the case of $`d_{max}>1000`$ where the Tornado code has a little higher threshold than our code, the iteration number is smaller for our code than for the Tornado code over a broad region of $`\delta `$, except near the threshold \[Fig. 3(b)\]. For an early convergence of decoding processes, each node needs to gather messages from other nodes quickly. This implies that graphs with smaller diameter may be more advantageous to reducing the iteration number. This in turn suggests that scale-free networks, which are known to have a very small diameter $`d\mathrm{ln}\mathrm{ln}N`$ Cohen  where $`N`$ is the number of nodes, may require a smaller iteration number than regular random networks or small-world networks Watt .

The error-correction capability of codes on scale-free networks can be further enhanced by adjusting the degree distribution, especially in the low degree region, so that it more closely models realistic scale-free networks whose degree distribution does not necessarily follow a power law for small $`k`$. While doing this, we try to keep as small as possible the number of parameters added to the generating function. After a number of numerical simulations we have found the following generating function adequate enough for this purpose:

$$G(x)=C\left[w_2p(2)x^2+w_3p(3)x^3+\underset{k=4}{\overset{d_{max}}{}}p(k)x^k\right],$$
(5)
where $`p(k)=(k+\alpha )^\gamma `$. The performance of this code is displayed in Table 2, which shows that by adding two new parameters, $`w_2`$ and $`w_3`$, which permit the two lowest degrees to vary from the power law, the performance of the code is increased. Addition of more parameters is expected to give rise to an increased performance, but at the expense of rendering the optimization process more time consuming.

In summary, we have found that many high-performance LDPC codes possess degree distributions well fitted by power-law functions with exponents close to 2. Based on this finding, we have developed codes on scale-free networks that have very good error-correction capabilities. The codes with power-law degree distribution yield better performance than those with exponential and Gaussian degree distributions that have fast decreasing tails. It also would be interesting to study the effect of degree correlations on the performance of a code, which is left as future work. As good error-correcting codes, the codes on scale-free networks could find lucrative applications in areas as diverse as wireless communication, media and data transfer over the Internet, and storage.


text:
# References and links

## 1 Introduction

The experimental investigation of nonlinear matter wave dynamics is feasible since the realization of Bose-Einstein-condensation of dilute gases . The combination of this new matter wave source with periodic potentials allows for the realization of many nonlinear propagation phenomena. The dynamics depends critically on the modulation depth of the potential. For deep periodic potentials the physics is described locally taking into account mean field effects and tunneling between adjacent potential wells. In this context wave packet dynamics in Josephson junction arrays have been demonstrated experimentally  and nonlinear self trapping has been predicted theoretically . In the limit of weak periodic potentials and moderate nonlinearity rich wave packet dynamics result due to the modification of dispersion which can be described applying band structure theory . Especially matter wave packets subjected to anomalous dispersion (negative effective mass) or vanishing dispersion (diverging mass) are of great interest. In the negative mass regime gap solitons have been predicted theoretically  and have been observed recently . Also modulation instabilities can occur .

The experiments described in this work reveal wave dynamics in the linear and nonlinear regime for weak periodic potentials. The observed behavior is a consequence of the special preparation of the wave packet leading to a continuous change of the effective mass and thus the dispersion during the propagation. The initial propagation is dominated by the atom-atom interaction leading to complex wave dynamics. After a certain time of propagation slowly spreading atomic wave packets are formed which are well described by linear theory. In this work we focus on the mechanisms governing the initial stage of propagation.

The paper is organized as follows: in section 2 we describe the effective mass and dispersion concept. In section 3 we present our experimental setup and in section 4 the employed wave packet preparation schemes are discussed in detail. In section 5 the experimental results are compared with numerical simulations. We show that some features of the complex dynamics can be identified with well known nonlinear mechanisms. We conclude in section 6.

## 2 Effective mass and dispersion concept

In our experiments we employ a weak periodic potential which leads to a dispersion relation $`E_n(q)`$ shown in Fig. 1(a). This relation is well known in the context of electrons in crystals  and exhibits a band structure. It shows the eigenenergies of the Bloch states as a function of the quasi-momentum $`q`$. The modified dispersion relation leads to a change of wavepacket dynamics due to the change in group velocity $`v_g(q)=1/\mathrm{}E/q`$ (see Fig. 1(b)), and the group velocity dispersion described by the effective mass $`m_{eff}=\mathrm{}^2(^2E/q^2)^1`$ (see Fig. 1(c)), which is equivalent to the effective diffraction introduced in the context of light beam propagation in optically-induced photonic lattices . In our experiment only the lowest band is populated, which is characterized by two dispersion regimes, normal and anomalous dispersion, corresponding to positive and negative effective mass. A pathological situation arises at the quasimomentum $`q_{\mathrm{}}^\pm `$, where the group velocity $`v_g(q)`$ is extremal, $`|m_{eff}|`$ diverges and thus the dispersion vanishes.

In the following we will show that the two preparation schemes employed in the experiment lead to a continuous change of the quasimomentum distribution, and thus to a continuous change of dispersion. One of the preparation schemes allows to switch periodically from positive to negative mass values and thus a slowly spreading wave packet is formed. This is an extension of the experiment reporting on dispersion management . The second preparation gives further insight into the ongoing nonlinear dynamics for the initial propagation.

## 3 Experimental Setup

The wave packets in our experiments have been realized with a 87Rb Bose-Einstein condensate (BEC). The atoms are collected in a magneto-optical trap and subsequently loaded into a magnetic time-orbiting potential trap. By evaporative cooling we produce a cold atomic cloud which is then transferred into an optical dipole trap realized by two focused Nd:YAG laser beams with $`60\mu m`$ waist crossing at the center of the magnetic trap (see Fig.2(a)). Further evaporative cooling is achieved by lowering the optical potential leading to pure Bose-Einstein condensates with $`110^4`$ atoms in the $`|F=2,m_F=+2`$ state. By switching off one dipole trap beam the atomic matter wave is released into a trap acting as a one-dimensional waveguide with radial trapping frequency $`\omega _{}=2\pi 100Hz`$ and longitudinal trapping frequency $`\omega _{}=2\pi 1.5Hz`$. It is important to note that the dipole trap allows to release the BEC in a very controlled way leading to an initial mean velocity uncertainty smaller than 1/10 of the photon recoil velocity.

The periodic potential is realized by a far off-resonant standing light wave with a single beam peak intensity of up to $`1W/cm^2`$. The chosen detuning of 2 nm to the blue off the D2 line leads to a spontaneous emission rate below $`1Hz`$. The standing light wave and the waveguide enclose an angle of $`\theta =21^{}`$ (see Fig. 2(b)). The frequency and phase of the individual laser beams are controlled by acousto-optic modulators driven by a two channel arbitrary waveform generator allowing for full control of the velocity and amplitude of the periodic potential. The light intensity and thus the absolute value of the potential depth was calibrated independently by analyzing results on Bragg scattering  and Landau Zener tunneling .

The wave packet evolution inside the combined potential of the waveguide and the lattice is studied by taking absorption images of the atomic density distribution after a variable time delay. The density profiles along the waveguide, $`n(x,t)`$, are obtained by integrating the absorption images over the transverse dimension.

## 4 Dynamics in reciprocal space

In our experimental situation an acceleration of the periodic potential to a constant velocity leads to a collective transverse excitation as indicated in Fig. 2(d). Since the transverse motion in the waveguide has a non vanishing component in the direction of the periodic potential due to the angle $`\theta `$, a change of the transverse velocity leads to a shift of the central quasimomentum of the wave packet. The coupling between the transverse motion in the waveguide and the motion along the standing light wave gives rise to a nontrivial motion in reciprocal (see Fig. 2(e,f)) and real space.

The appropriate theoretical description of the presented experimental situation requires the solution of the three dimensional nonlinear Schrödinger equation (NLSE) and thus requires long computation times. In order to understand the basic physics we follow a simple approach which solves the problem approximately and explains all the features observed in the experiment. For that purpose we first solve the semiclassical equations of motion of a particle which obeys the equation $`\stackrel{}{F}=M^{}\ddot{\stackrel{}{x}}`$ where $`M^{}`$ is a mass tensor describing the directionality of the effective mass. We deduce the time dependent quasimomentum $`q_c(t)`$ in the direction of the periodic potential by identifying $`\mathrm{}\dot{q_c}=F_{\widehat{x}}`$ and $`\dot{\widehat{x}}=v_g(q_c)`$ (definition of $`\widehat{x}`$ see Fig. 2(b)). Subsequently we can solve the one dimensional NPSE (non-polynomial nonlinear Schrödinger equation) where the momentum distribution is shifted in each integration step according to the calculated $`q_c(t)`$. Thus the transverse motion is taken into account properly for narrow momentum distributions. We use a split step Fourier method to integrate the NPSE where the kinetic energy contribution is described by the numerically obtained energy dispersion relation of the lowest band $`E_0(q)`$. It is important to note, that this description includes all higher derivatives of $`E_0(q)`$, and thus goes beyond the effective mass approximation.

In the following we discuss in detail the employed preparation schemes:

Acceleration scheme I: After the periodic potential is adiabatically ramped up to $`V_0=6E_{rec}`$ it is accelerated within $`3`$ms to a velocity $`v_{pot}=\mathrm{cos}^2(\theta )1.5v_{rec}`$. Then the potential depth is lowered to $`V_0=0.52E_{rec}`$ within $`1.5`$ms and the periodic potential is decelerated within $`3`$ms to $`v_{pot}=\mathrm{cos}^2(\theta )v_{rec}`$ subsequently. $`V_0`$ and $`v_{pot}`$ are kept constant during the following propagation. The calculated motion in reciprocal space $`q_c(t)`$ is shown in Fig. 2(e).

Acceleration scheme II: The periodic potential is ramped up adiabatically to $`V_0=0.37E_{rec}`$ and is subsequently accelerated within $`3`$ms to a final velocity $`v_{pot}=\mathrm{cos}^2(\theta )\times 1.05v_{rec}`$. The potential depth is kept constant throughout the whole experiment. Fig. 2(f) reveals that in contrast to the former acceleration scheme the quasimomentum for the initial propagation is mainly in the negative effective mass regime.

## 5 Experimental and Numerical Results

In this section we compare the experimental results with the predictions of our simple theoretical model discussed above. The numerical simulation reveal all the experimentally observed features of the dynamics such as linear slowly spreading oscillating wave packets, nonlinear wave packet compression and splitting of wave packets. The observed nonlinear phenomena can be understood by realizing that in the negative effective mass regime the repulsive atom-atom interaction leads to compression of the wave packet in real space and to a broadening of the momentum distribution. An equivalent picture borrowed from nonlinear photon optics  is the transient formation of higher order solitons, which show periodic compression in real space with an increase in momentum width and vice versa.

### 5.1 Preparation I

The experimental results for the first acceleration scheme discussed in section 4 are shown in Fig. 3. Clearly we observe that a wave packet with reduced density is formed which spreads out slowly and reveals oscillations in real space. This wave packet results from the initial dynamics characterized by two stages of compression which lead to radiation of atoms . The observed behavior is well described by our numerical simulation which allows further insight into the ongoing physics.

In Fig. 3(c,d) we show the calculated momentum and real space distribution for the first 14ms of propagation. As can be seen the acceleration of the standing light wave leads to a oscillatory behavior in momentum space. For the chosen parameters the wave packet is initially dragged with a tight binding potential ($`V_0=6E_{rec}`$) over the critical negative mass regime. While the real space distribution does not change during this process, the momentum distribution broadens due to self phase modulation . The subsequent propagation in the positive mass regime leads to a further broadening in momentum space and real space (t=4-9ms).

The dynamics changes drastically as soon as a significant part of the momentum distribution populates quasimomenta in the negative mass regime (t=10ms). There the real space distribution reveals nonlinear compression as known from the initial dynamics of higher order solitons. This compression leads to a significant further broadening in momentum space and thus to population of quasimomenta corresponding to positive mass. This results in a spreading in real space due to the different group velocities involved and leads to the observed background. The change of the quasimomentum due to the transverse motion prohibits a further significant increase in momentum width, since the whole momentum distribution is shifted out of the critical negative mass regime at t=14ms.

The long time dynamics of the slowly spreading wave packet is mainly given by the momentum distribution marked with the shaded area for t=14ms in Fig. 3(c). The subsequent motion is dominated by the change of the quasimomentum due to the transverse motion. This leads to a periodic change from normal to anomalous dispersion and thus the linear spreading is suppressed. This is an extension of our previous work on dispersion management for matter waves - continuous dispersion management.

### 5.2 Preparation II

This preparation scheme reveals in more detail the transient solitonic propagation leading to the significant spreading in momentum space. This results in a splitting of the wave packet which cannot be understood within a linear theory. The results are shown in Fig. 4 and the observed splitting is confirmed by our numerical simulations.

In contrast to the former preparation scheme the momentum distribution is prepared as a whole in the critical negative mass regime. Our numerical simulations reveal that the wave packet compresses quickly in real space after t=4ms which is accompanied by an expansion in momentum space. The momentum distribution which stays localized in the negative mass regime reveals further solitonic propagation characterized by an expansion in real space and narrowing of the momentum distribution (t=5-10ms). The transverse motion shifts this momentum distribution into the normal dispersion regime after 11ms of propagation resulting in a wave packet moving with positive group velocity (i.e. moving to the right in fig. 4(b)). The initial compression at t=4ms even produces a significant population of atoms in the normal mass regime which subsequently move with negative group velocity showing up as a wave packet moving to the left in Fig. 4(b). Thus the splitting in real space is a consequence of the significant nonlinear broadening in momentum space.

## 6 Conclusion

In this paper we report on experimental observations of nonlinear wave packet dynamics in the regime of positive and negative effective mass. Our experimental setup realizing a BEC in a quasi-one dimensional situation allows the observation of wave dynamics for short times, where the nonlinearity due to the atom-atom interaction dominates and also for long times, where linear wave propagation is revealed.

We have shown that a slowly spreading wave packet can be realized by changing the quasimomentum periodically from the normal to anomalous dispersion regime. This can be viewed as an implementation of continuous dispersion management. We further investigate in detail the formation process of these packets, which are a result of the initial spreading in momentum space due to nonlinear compression. A second experiment investigates in more detail the nonlinear dynamics in the negative mass regime where the solitonic propagation leads to a significant broadening in momentum space. This shows up in the experiment as splitting of the condensate into two wave packets which propagate in opposite directions.

The developed theoretical description utilizing the effective mass tensor models the experimental system in one dimension and can explain all main features observed in the experiment.

## 7 Acknowledgment

This work was supported by Deutsche Forschungsgemeinschaft, Emmy Noether Program, by the European Union, Contract No. HPRN-CT-2000-00125, and the Optik Zentrum University of Konstanz.