text:
# The Absolute Abundance of Iron in the Solar Corona

## 1 INTRODUCTION

It seems remarkable that the absolute abundances (here taken to mean abundances relative to hydrogen) of trace elements in the solar corona are still a matter for debate when many measurements for much more distant stars exist. The reason for this situation is easy to understand. Classic spectroscopic techniques for measuring relative abundances require that we compare the amount of radiation emitted in spectral lines whose atomic physics we understand and which arise from the two different elements in question. In the solar corona H and He are fully ionized, and hence produce no lines. In stars with hotter coronae, thermal bremsstrahlung produces a continuum emission which can be used instead of a spectral line as a measure of the amount of H present, but the Sun’s quiescent corona is so cool that the continuum present at X–ray or EUV wavelengths is relatively weak and difficult to measure.

It has been known for some time that there is a difference in the relative abundances of elements in the photosphere and the corona organized according to the first ionization potential (FIP) of the element: high–FIP elements such as O are underabundant in the corona relative to low–FIP elements such as Fe by a factor of 4 (Meyer, 1985). However, it has not been clear whether low–FIP elements are enhanced in the corona in an absolute sense (i.e., as a fraction of the total ion number density, dominated by H and He) or the high–FIP elements are depleted relative to H. Two main techniques have been used to infer absolute abundances in the solar corona (Fludra et al., 1999). One can observe solar flares, which produce plasma hot enough for the continuum to be measured, and compare the continuum with a suitable set of lines; or one can measure particles directly in the solar wind, and assume that the abundances there truly reflect the abundances in the solar corona. Unfortunately, these two techniques have resulted in conflicting measurements for, e.g., the absolute abundance of iron. The flare measurements suggest that it is close to the absolute abundance in the photosphere (Veck & Parkinson, 1981) or perhaps factors of 1.5 – 2 larger (Fludra & Schmelz, 1999), while the solar wind particle measurements suggest that it is 4 times the photospheric ratio (Reames, 1992).

In this paper we use a third technique to measure the absolute abundance $`N_{\mathrm{Fe}}/N_\mathrm{H}`$ which has several desirable features. We use EUV observations of lines of Fe in order to determine the amount of Fe present in the corona, but to measure the amount of H present we use radio continuum measurements of thermal bremsstrahlung from coronal plasma. This measurement therefore pertains to the quiescent solar corona, including non-flaring weak active regions, rather than flare or solar wind plasma. Like classic spectroscopic techniques, this technique should be independent of the filling factor of plasma in the volume studied since both types of emission have the same dependence on number density. We describe the measurement technique and the observations used for this measurement in the following section. We then discuss the implications of our measurement.

## 2 THE MEASUREMENT TECHNIQUE

The method proposed here to measure the abundances of elements relative to H rests on the fact that the EUV line emission and radio bremsstrahlung both depend on the product of the electron number density and the number density of an element. In the case of radio bremsstrahlung the element is effectively hydrogen; in the case of an EUV line, it is the element responsible for the line. Therefore comparison of the two yields the abundance of the latter element relative to hydrogen.

Specifically, since most EUV lines produced in the solar corona are collisionally excited and optically thin, the line intensity at a photon energy $`E`$ corresponding to a transition in ionization state Z of element X is the integral over temperature of the product $`Q_{X,Z,E}DEM(T)N_X/N_\mathrm{H}`$ (Fludra et al., 1999), where $`Q_{X,Z,E}`$ is the “contribution function”: it includes the intrinsic spontaneous emission rate (the Einstein “A” coefficient) as well as the ratio $`N_{X,Z}/N_X`$ of the number density of the particular charge state $`Z`$ of the element $`X`$ responsible for the line to the total number density of the element. This ratio is usually determined by an ionization balance calculation appropriate to the temperature of the plasma, and $`Q`$ is a strong function of temperature for this reason. However, for most lines $`Q`$ is effectively independent of $`N_X`$ and $`N_e`$. The (hydrogen column) differential emission measure DEM(T) is the integral along the line of sight of $`d(N_eN_\mathrm{H})/dT`$.

Bremsstrahlung is one of the two main mechanisms for radio emission in the solar corona and also has a significant dependence on temperature: in the optically thin limit, the flux at a frequency $`f`$ is

$$S=9.78\times \mathrm{\hspace{0.17em}10}^3\frac{2k_B}{c^2}(1+4\frac{N_{\mathrm{He}}}{N_\mathrm{H}})T^{0.5}DEM(T)G(T)𝑑T𝑑\mathrm{\Omega }$$
(1)
in ergs cm-2 s-1 Hz-1 (Dulk, 1985), where $`G(T)`$ = 24.5 - ln($`T/f`$) is the Gaunt factor, $`k_B`$ is the Boltzmann constant, $`c`$ the speed of light and the second integral is over the solid angle subtended by the source. We are assuming that in a plasma at the temperature of the solar corona the ion component is completely dominated by protons and fully–ionized helium, and we adopt $`N_{\mathrm{He}}/N_\mathrm{H}=\mathrm{\hspace{0.17em}0.1}`$. It is clear from (1) that we must know the DEM in order to calculate the thermal bremsstrahlung radio flux accurately. Note that as long as $`N_X/N_\mathrm{H}`$ does not vary from place to place, (1) has exactly the same dependence on density as the EUV line emissivity and thus filling factor has exactly the same effect on the EUV line and radio continuum fluxes.

We determine the abundance $`N_{\mathrm{Fe}}/N_\mathrm{H}`$ with the following steps: (i) We measure the radio flux from a source in the solar corona. This source must be optically thin and its radio flux must be due only to bremsstrahlung. The latter requirement can be tested by observing the radio spectrum of the source, which should be essentially flat for optically thin bremsstrahlung. (ii) We observe the same region with the Coronal Diagnostic Spectrometer (CDS) on SOHO in order to measure the intensities of lines from a range of charge states of Fe sufficient to give a reliable determination of the DEM. (iii) We use the CHIANTI atomic data package (Dere et al., 1997; Landi et al., 1999) to fit the Fe line fluxes to a DEM. (iv) We assume that the shape of the DEM for electrons and protons is the same as that for Fe, and integrate over temperature in (1) to predict the optically–thin bremsstrahlung radio flux we expect to observe, for a given absolute abundance of Fe. Comparison with the radio flux measured by the VLA yields $`N_{\mathrm{Fe}}/N_\mathrm{H}`$.

Coordinated observations with CDS and the VLA111The Very Large Array is run by the National Radio Astronomy Observatory which is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. were carried out on 1997 July 27, Aug. 3 and Nov. 11. We will present an analysis of all the data elsewhere; here we discuss the best target observed, a small active region AR 8105 found at about S35W15 on 1997 Nov. 11 (with the VLA in “D” configuration). This was the best target in the sense that it was compact: the VLA is an interferometer and cannot measure the flux of large sources, so radio flux measurements are more reliable for compact sources. The region was one of four observed on this day with the VLA at 1.4, 4.8 and 8.4 GHz. CDS also observed the same 4 regions alternately using the sequence FE\_IONS which reads nine 21–pixel–wide windows from the NIS1 detector of CDS containing 14 identified lines including those of Fe X, XI, XII, XIII, XIV, XV, XVI and XVII. One run of the sequence took 32 minutes, scanning with the 4″ slit across a 240″ wide field of view. Five sequences were acquired on AR 8105 centered at about 13:14, 16:28, 18:42, 20:56 and 23:10 UT, respectively. The VLA observations covered the period from 14:30 UT to 23:10 UT, cycling between the different frequencies rapidly in order to achieve optimal $`uv`$ coverage. Images of AR 8105 acquired by CDS, the VLA and EIT are shown in Figure 1.

## 3 THE ABUNDANCE OF IRON

It is essential for this measurement that the VLA see the same plasma that is seen in the EUV, and Fig. 1 shows that this is indeed the case: the EUV and radio images match each other extremely well, indicating that gyroresonance opacity does not contribute significantly to the radio emission and we conclude that it is due to bremsstrahlung. The higher–resolution EIT images shows that two east–west loops dominated AR 8105.

The VLA observations resulted in images with resolutions of 30″ at 1.4 GHz, 13″ at 4.8 GHz and 7″ at 8.4 GHz. Using the 4.8 and 8.4 GHz images we identified an area of dimension 69″ $`\times `$ 97″ around AR 8105 which includes essentially all the flux from the region. Radio fluxes were derived by summing over this area at 4.8 and 8.4 GHz, while due to the poorer resolution at 1.4 GHz we instead fitted a Gaussian model to the feature associated with AR 8105 and determined its flux. At each frequency we have independent measurements in two separate 50 MHz bands and in the two circular polarizations, all of which are consistent with one another. The resulting fluxes are 0.25 $`\pm `$ 0.01 sfu at 1.4 GHz, 0.207 $`\pm `$ 0.004 sfu at 4.8 GHz and 0.197 $`\pm `$ 0.006 sfu at 8.4 GHz, which are consistent with the $`f^{0.1}`$ spectrum expected for optically thin bremsstrahlung due to the Gaunt factor in (1). The quoted uncertainties are the map rms per beam times the square root of the number of beams over which integration was carried out, and are much larger than the formal uncertainties in the fluxes. Polarizations were also low, consistent with the properties of bremsstrahlung.

The CDS data were calibrated using the standard reduction package in SOLARSOFT (results presented here use the version current in 1999 December). Line profiles were fitted and fluxes extracted, with uncertainties calculated according to Thompson (1998). We used both a recent radiometric CDS calibration derived from comparison with the SERTS rocket flight in 1997 (Thomas et al., 2000) and the “version 2” CDS calibration which was in the CDS software as of 1999 December. One important step is that background subtraction was carried out on the CDS data using a linear fit to pixels to the north and south of the region studied, well away from AR 8105. This step is important because of the fact that the VLA does not measure large–scale flux and hence is insensitive to the quiet–Sun contribution. In effect we assume that a smooth quiet–Sun contribution has been subtracted from the radio fluxes, and hence it should be subtracted from the CDS data as well. This predominantly affects the cooler lines, amounting to up to 20% of the photons in the original spectrum, but being fairly uniform across the CDS field of view.

The following lines (with typical line fluxes averaged over the region shown in parentheses, in ergs cm-2 s-1 st-1) were then used to calculate the DEM: Fe X at 345.7 Å (55 $`\pm `$ 3), Fe XI at 352.7 Å (140 $`\pm `$ 4), Fe XII at 352.1 Å (145 $`\pm `$ 4) and 364.4 Å (230 $`\pm `$ 4), Fe XIII at 348.3 Å (95 $`\pm `$ 4), Fe XIV at 334.2 Å (260 $`\pm `$ 4), Fe XVI at 360.7 Å (380 $`\pm `$ 4) and Fe XVII (undetected here) at 350.5 Å. These lines were chosen to be as density–insensitive as possible while still covering the temperature range expected for the non–flaring solar corona, log T = 5.8 up to 6.8. The Fe XIII and Fe XII lines do have significant sensitivity to density in the likely range of coronal densities, 109 – 1010 cm-3. In practice, AR 8105 was quiet and there is no significant plasma outside the temperature range 5.9 – 6.4. The fitting procedure uses bins of width 0.1 in log T; in effect, we find that only 4 bins contain significant emission measure and the adjacent two bins are also significant, so we are fitting to 6 parameters using 8 lines from 7 different charge states. The Fe XV line at 327.1 Å was not used to determine the DEM because it was always found to be $``$ 30% lower than the models predicted; the atomic data for it are clearly not consistent with the strong Fe XIV and XVI lines observed. Additional details will be given in a companion paper.

We adopted the ionization equilibrium calculations of Arnaud & Raymond (1992) and assume a constant pressure $`nT`$ of 1 $`\times `$ 1016 cm-3 K (e.g., $`N`$ = 5 $`\times `$ 109 cm-3 at $`T`$ = 2 $`\times `$ 106 K) in determining the DEM. The line fluxes used for the fit were generally reproduced by the resulting model DEM to better than 10% (SERTS) or 5% (“version 2”). Figure 2 shows the geometric mean DEM over the 5 CDS sequences with error bars showing the standard deviation in each temperature bin. Essentially no plasma is present below log T = 6.0 (controlled by the ratio of Fe X to Fe XI) while the amount of plasma above 2 $`\times `$ 106 K is small (controlled by the ratio of Fe XVI to Fe XIV). These properties are consistent with those of DEMs derived for other regions of similarly low activity (Brosius et al., 1996). Integrating over each of the 5 DEMs resulted in a predicted radio flux of 0.829 $`\pm `$ 0.023 sfu in the 6700 arcsec2 area at 4.8 GHz if we use the CDS calibration based on the SERTS–97 comparison and adopt $`N_{\mathrm{Fe}}/N_\mathrm{H}`$ = 3.9 $`\times `$ 10-5 as in Meyer (1985), which is close to the photospheric value. The observed radio flux of 0.207 sfu is a factor of 4.0 smaller than this. Consequently we derive a value of 1.56 $`\times `$ 10-4 for $`N_{\mathrm{Fe}}/N_\mathrm{H}`$, or $`\mathrm{log}_{10}A_{\mathrm{Fe}}`$ = 8.19 on the commonly used scale where $`\mathrm{log}_{10}A_\mathrm{H}`$ = 12.00. Using the “version 2” CDS calibration we predict 0.716 $`\pm `$ 0.013 sfu and $`N_{\mathrm{Fe}}/N_\mathrm{H}`$ = 1.35 $`\times `$ 10-4 ($`\mathrm{log}_{10}A_{\mathrm{Fe}}`$ = 8.13). We adopt the more recent SERTS–97 calibration for our final result. The uncertainty of 3% for the predicted radio flux is the standard error in the mean of the results of the five distinct CDS observations of AR 8105 and thus incorporates the effects of both temporal variability over the 12 hours of the experiment as well as the stochastic nature of the process of fitting the line fluxes to a DEM. We adopt an uncertainty of 15% for the CDS radiometric calibration (Thomas et al., 2000).

It turns out that the measured value for $`N_{\mathrm{Fe}}/N_\mathrm{H}`$ is reasonably robust to several of the assumptions made here, mostly because it depends on an integration over the DEM and is thus not very sensitive to details of the DEM. One assumption is that the result quoted uses the version of CHIANTI (0.9) used by the CDS team at NASA/GSFC as of 1999 December, which is not the most current version (2.0). For the Fe lines used here the only differences between the two versions appear to be in Fe X and Fe XI, and the newer version results in significantly poorer fits for the DEM because of these two lines: the resulting DEMs have less emission measure at low temperatures and more at high temperatures. Notwithstanding this fact, the result for the abundance differs by less than 3% with the newer version of CHIANTI. If we use the ionization equilibrium calculations of Arnaud & Rothenflug (1985), the DEMs show a pronounced peak at $`\mathrm{log}T`$ = 6.3 and the resulting abundance is about 10% larger. A change in the assumed pressure by a factor of 2 affects the results by a similar amount: it would be 8% smaller at 5.0 $`\times `$ 1015 cm-3 K and 10% larger at 2.0 $`\times `$ 1016 cm-3 K. We attempt to minimize the uncertainty due to variation in pressure across the active region by choosing a value representative of the source region for the bulk of the EUV line emission. We note also that the background subtraction carried out on the CDS data introduces some uncertainty in the final results, as does solar calibration at the VLA, which uses dedicated noise sources for solar calibration which need to be measured routinely to maintain calibration accuracy. However, the calibration at each of the three frequencies observed is independent of the others, so we take the fact that the measured spectrum is flat to be an indication that solar calibration was not in error by more than the 20% levels which all the other factors may introduce.

## 4 CONCLUSIONS

Our measurement of an absolute abundance for Fe which is 4 times photospheric is clearly consistent with the interpretation that low–FIP elements such as Fe are enhanced relative to their photospheric abundances. Of the commonly–used solar abundance tables, Feldman’s (1992) is closest to our measurement with $`\mathrm{log}_{10}A_{\mathrm{Fe}}`$ = 8.10; his value is based more on the argument that H should be treated as a high–FIP element rather than any specific measurement of $`N_{\mathrm{Fe}}/N_\mathrm{H}`$. Other commonly used values for $`\mathrm{log}_{10}A_{\mathrm{Fe}}`$ are the photospheric results 7.60 (Allen, 1973), 7.51 (Anders & Grevesse, 1989) and 7.67 (Grevesse & Anders, 1991), the coronal results 7.59 (Meyer, 1985), 7.83 (Fludra & Schmelz, 1999) and 8.50 (Waljeski et al., 1994), and the solar wind measurements 7.93 (Reames, 1995) and 8.51 (Reames, 1999). We note in passing that spectroscopic measurements of the Fe abundance in active stars with coronae much hotter than the Sun’s have tended to suggest abundances smaller than photospheric (Mewe et al., 1996; Drake, 1996), in contrast to the solar behaviour we find here.

We have reported here a single measurement of the Fe abundance. Spectroscopic studies of relative abundances suggest that they may vary considerably from flare to flare and from one region to the next (Sylwester et al., 1984; Saba, 1995; Sylwester et al., 1998), so it is clearly important to confirm our measurement. We have additional data sets which will be published in a longer paper elsewhere, but with larger uncertainties expected in the radio fluxes, and expect to carry out further studies of this type in the future.

## Acknowledgments

This research at the University of Maryland was carried out primarily under a SOHO GI grant from NASA, NAG 5–4954. Solar radiophysics at the University of Maryland is also supported by NSF grant ATM 96-12738 and NASA grants NAG 5-7370, NAG 5-6257 and NAG 5-7901. We thank Richard Harrison, Andrzej Fludra and the CDS operators for their assistance with the project, and Barry Clark of NRAO for scheduling the VLA observations to match those of SOHO. We thank the CHIANTI consortium for their efforts to make atomic data readily accessible.


format: text

text:
# FERMIONS ON THE LIGHT-FRONT

## 1 Spontaneous Symmetry Breaking

The first of the above issues has been addressed very often in the past and we will restrict ourselves here to a brief summary.222See for example Ref.  for a more extensive discussion on this question. In the LF framework, non-trivial vacuum structure can reside only in zero-modes ($`k^+=0`$ modes). Since these are high-energy modes (actually infinite energy in the continuum) one often does not include them as explicit degrees of freedom but assumes they have been integrated out, leaving behind an effective LF-Hamiltonian. The important points here are the following. If the zero-mode sector involves spontaneous symmetry breaking, this manifests itself as explicit symmetry breaking for the effective Hamiltonian. In general, these effective LF Hamiltonians thus have a much richer operator structure than the canonical Hamiltonian. Therefore, compared to a conventional Hamiltonian framework, the question of the vacuum has been shifted from the states to the operators and it should thus be clear that the issues of renormalization and the vacuum are deeply entangled in the LF framework.

## 2 Current Conservation

Despite widespread confusion on this subject, vector current conservation (VCC) is actually not a problem in the LF framework. Many researchers avoid the subject of current conservation because the divergence of the vector current

$$q_\mu j^\mu (q)=q^{}j^++q^+j^{}\stackrel{}{q}_{}\stackrel{}{j}_{}$$
(1)
involves $`j^{}`$, which is quadratic in the constrained fermion spinor component

$$j^{}=\psi _{()}^{}\psi _{()},$$
(2)
and thus $`j^{}`$ contains quartic interactions making it at least as difficult to renormalize as the Hamiltonian.

We know already from renormalizing $`P^{}`$ that the canonical relation between $`\psi _{()}`$ and $`\psi _{(+)}`$ is in general not preserved in composite operators (such as $`\overline{\psi }\psi `$). It is therefore clear that a canonical definition for $`j^{}`$ will in general violate current conservation since it does not take into account integrating out zero-modes and other high-energy degrees of freedom.

However, in LF gauge $`A^+=0`$, $`j^{}`$ does not enter the Hamiltonian and therefore one can address its definition separately from the construction of the Hamiltonian. In fact, it is very easy to find a pragmatic definition which obviously guarantees manifest current conservation, namely

$`j^{}(q^+)`$ $`=`$ $`{\displaystyle \frac{1}{q^+}}[P^{},j^+(q^+)](1+1)`$ (3)
$`j^{}(q^+,\stackrel{}{q}_{})`$ $`=`$ $`{\displaystyle \frac{1}{q^+}}\left\{[P^{},j^+(q^+,\stackrel{}{q}_{})]\stackrel{}{q}_{}\stackrel{}{j}_{}(q^+,\stackrel{}{q}_{})\right\}(3+1)`$
in 1+1 and 3+1 dimensions respectively. The corresponding expressions in coordinate space $`(x^{},\stackrel{}{x}_{})`$ can be obtained by Fourier transform. In summary,

* Most importantly, VCC is no problem in LF quantization and is manifest at the operator level, provided $`j^{}`$ is defined using Eq. (3).
* Since VCC can easily be made manifest (by using the above construction!), there is no point in testing its validity and it cannot be used as a non-perturbative renormalization condition either.
* For a non-interacting theory, Eq. (3) reduces to the canonical definition of $`j^{}`$, but in general this is not the case when $`P^{}`$ contains interactions or even non-canonical terms.

Note that (as so often on the LF) $`q^+`$ appears in the denominator of Eq. (3). Therefore, as usual, one should be very careful while taking the $`q^+0`$ limit and while drawing any conclusions about this limit. An example of this kind are the pair creation terms in $`j^{}`$. Naively they do not contribute for $`q^+0`$, since the $`q\overline{q}`$ pair emanating from $`j^{}`$ necessarily carries positive $`q^+`$. However, since $`j^{}`$ often has very singular matrix elements for $`q^+`$, such seemingly vanishing terms nevertheless survive the $`q^+0`$ limit, which often leads to confusion. For an early example of this kind see Ref. . A more recent discussion can be found in Ref. .

## 3 Parity Invariance

### 3.1 General Remarks

A parity transformation, $`x^0\stackrel{P}{}x^0`$, $`\stackrel{}{x}\stackrel{P}{}\stackrel{}{x}`$ leaves the quantization hyperplane ($`x^0=0`$) in equal time (ET) quantization invariant and therefore the parity operator is a kinematic operator in such a framework. It is thus very easy to ensure that parity is a manifest symmetry in ET quantization by tracking parity at each step in a calculation. The situation is completely different in the LF framework, where the same parity transformation exchanges LF-‘time’ ($`x^+x^0+x^3`$) and space ($`x^{}x^0x^3`$) directions, i.e.

$$x^+\stackrel{P}{}x^{}$$
(4)
and therefore the quantization hyperplane $`x^+=0`$ is not invariant. Hence, the parity operator is a dynamical operator on the LF and, except for a free field theory, it is probably impossible to write down a simple expression for it in terms of quark and gluon field operators. Thus, parity invariance is not a manifest symmetry in this framework. Note that the situation is the other way round for the boost operator (kinematic and manifest on the LF, dynamical and non-manifest in ET). It thus depends on the physics application one is interested in and the symmetries that one considers the most important ones for that particular physics problem, which framework is preferable.

For most applications of LF quantization, the lack of manifest parity actually does not constitute a problem — in fact, one can view it as an opportunity rather than a problem. The important point here is the following: due to the lack of manifest covariance in a Hamiltonian formulation, LF Hamiltonians in general contain more parameters than the corresponding Lagrangian. Parity invariance may be very sensitive to some of these parameters.

In order to illustrate this important point, let us consider the example of a 1+1 dimensional Yukawa model

$$=\overline{\psi }\left(i\overline{)}m_Fg\varphi \gamma _5\right)\psi \frac{1}{2}\varphi \left(\mathrm{}+m_B^2\right)\varphi .$$
(5)
This model actually has a lot in common with the kind of interactions that appear when one formulates QCD (with fermions) on a transverse lattice .

The main difference between scalar and Dirac fields in the LF formulation is that not all components of the Dirac field are dynamical: multiplying the Dirac equation

$$\left(i\overline{)}m_Fg\varphi \gamma _5\right)\psi =0$$
(6)
by $`\frac{1}{2}\gamma ^+`$ yields a constraint equation (i.e. an “equation of motion” without a time derivative)

$$i_{}\psi _{()}=\left(m_F+g\varphi \gamma _5\right)\gamma ^+\psi _{(+)},$$
(7)
where $`\psi _\pm \frac{1}{2}\gamma ^{}\gamma ^\pm \psi .`$ For the quantization procedure, it is convenient to eliminate $`\psi _{()}`$, using

$$\psi _{()}=\frac{\gamma ^+}{2i_{}}\left(m_F+g\varphi \gamma _5\right)\psi _{(+)}$$
(8)
from the classical Lagrangian before imposing quantization conditions, yielding

$``$ $`=`$ $`\sqrt{2}\psi _{(+)}^{}i_+\psi _{(+)}\varphi \left(\mathrm{}+m_B^2\right)\varphi \psi _{(+)}^{}{\displaystyle \frac{m_F^2}{\sqrt{2}i_{}}}\psi _{(+)}`$
$``$ $`\psi _{(+)}^{}\left(g\varphi {\displaystyle \frac{m_F\gamma _5}{\sqrt{2}i_{}}}+{\displaystyle \frac{m_F\gamma _5}{\sqrt{2}i_{}}}g\varphi \right)\psi _{(+)}\psi _{(+)}^{}g\varphi {\displaystyle \frac{1}{\sqrt{2}i_{}}}g\varphi \psi _{(+)}.`$
The rest of the quantization procedure very much resembles the procedure for self-interacting scalar fields.

The above canonical Hamiltonian contains a kinetic term for the fermions, a fermion boson vertex and a fermion 2-boson vertex. While the couplings of these three terms in the canonical Hamiltonian depend only on two independent parameters ($`m`$ and $`g`$), it turns out that these terms are renormalized independently from each other once zero-mode and other high-energy degrees of freedom are integrated out. More explicitly this means that one should make an ansatz for the renormalized LF Hamiltonian density of the form

$`𝒫^{}`$ $`=`$ $`{\displaystyle \frac{m_B^2}{2}}\varphi ^2+\psi _{(+)}^{}{\displaystyle \frac{c_2}{\sqrt{2}i_{}}}\psi _{(+)}+c_3\psi _{(+)}^{}\left(\varphi {\displaystyle \frac{\gamma _5}{\sqrt{2}i_{}}}+{\displaystyle \frac{\gamma _5}{\sqrt{2}i_{}}}g\varphi \right)\psi _{(+)}`$ (10)
$`+`$ $`c_4\psi _{(+)}^{}\varphi {\displaystyle \frac{1}{\sqrt{2}i_{}}}\varphi \psi _{(+)},`$
where the $`c_i`$ do not necessarily satisfy the canonical relation $`c_3^2=c_2c_4`$. However, this does not mean that the $`c_i`$ are completely independent from each other. In fact, Eq.(10) will describe the Yukawa model only for specific combinations of $`c_i`$. It is only that we do not know the relation between the $`c_i`$. 333Coupling coherence  does not help to determine the finite part (integration constant) for the mass if the original Lagrangian contains a nonzero bare mass!

Thus the bad new is that the number of parameters in the LF Hamiltonian has increased by one (compared to the Lagrangian). The good news is that a wrong combination of $`c_i`$ will in general give rise to a parity violating theory: formally this can be seen in the weak coupling limit, where the correct relation ($`c_3^2=c_2c_4`$) follows from a covariant Lagrangian. Any deviation from this relation can be described on the level of the Lagrangian (for free massive fields, equivalence between LF and covariant formulation is not an issue) by addition of a term of the form $`\delta =\overline{\psi }\frac{\gamma ^+}{i^+}\psi `$, which is obviously parity violating, since parity transformations result in $`A^\pm \stackrel{P}{}A^{}`$ for Lorentz vectors $`A^\mu `$; i.e. $`\delta \stackrel{P}{}\overline{\psi }\frac{\gamma ^{}}{i^{}}\psi \delta `$. This also affects physical observables, as can be seen by considering boson fermion scattering in the weak coupling limit of the Yukawa model. At the tree level, there is an instantaneous contact interaction, which is proportional to $`\frac{1}{q^+}`$. The (unphysical) singularity at $`q^+=0`$ is canceled by a term with fermion intermediate states, which contributes (near the pole) with an amplitude $`\frac{m_v^2}{m_{kin}^2}\frac{1}{q^+}`$. Obviously, the singularity cancels iff $`m_V=m_{kin}`$. Since the singularity involves the LF component $`q^+`$, this singular piece obviously changes under parity. This result is consistent with the fact that there is no zero-mode induced renormalization of $`m_{kin}`$ and thus $`m_V=m_{kin}`$ at the tree level. This simple example clearly demonstrates that a ‘false’ combination of $`m_V`$ and $`m_{kin}`$ leads to violations of parity for a physical observable, which is why imposing parity invariance as a renormalization condition may help reduce the dimensionality of coupling constant space.

### 3.2 Parity Sensitive Observables that “don’t work”

Of course there are an infinite number of parity sensitive observables, but not all of them are easy to evaluate non-perturbatively in the Hamiltonian LF framework. Furthermore, we will see below that some relations among observables, which seem to be sensitive to parity violations, are actually ‘protected’ by manifest symmetries, such as charge conjugation, or by VCC.

From the brief discussion of the (to-be-canceled) singularity above it seems that the most sensitive observable to look for violations of parity in QCD would be Compton scattering cross sections between quarks and gluons (or the corresponding fermions and bosons in other field theories) because there one could tune the external momenta such that the potential singularity enters with maximum strength. However, this is not a very good choice: first of all non-perturbative scattering amplitudes are somewhat complicated to construct on the LF. Secondly, quarks and gluons are confined particles, which makes $`qg`$ Compton scattering an unphysical process.

A much better choice are any matrix elements between bound states. Bound states are non-perturbative and all possible momentum transfers occur in their time evolution. Therefore, any parity violating sub-amplitude would contribute at some point and would therefore affect physical observables. Secondly, since one of the primary goal of LFQCD is to explore the non-perturbative spectra and structure of hadrons, matrix elements in bound states are the kind of observables for which the whole framework has been tailored.

One conceivable set of matrix elements are those of the vector current operator $`j^\pm `$. For example, consider the vacuum to meson matrix elements

$`0|j^+|n,p`$ $`=`$ $`p^+f_n^{(+)}`$
$`0|j^{}|n,p`$ $`=`$ $`p^{}f_n^{()}.`$ (11)
By boost invariance, the couplings defined in Eq. (11) must be independent of the momenta. Obviously, parity invariance requires $`|f_n^{(+)}|=|f_n^{()}|`$. However, this relation also follows from current conservation $`0=p^{}0|j^+|n,p+p^+0|j^{}|n,p`$, which makes it a useless relation for the purpose of parity tests.

Similar statements can be made about elastic formfactors but we will omit the details here. The basic upshot is that the same relations between the matrix elements of $`j^\pm `$ that arise from parity invariance can often also be derived using only VCC, i.e. such relations are in general ‘protected’ by VCC

One may also consider non-conserved currents, such as the axial vector current. However, there one would have to face the issue of defining the ‘minus’ component before one can test any parity relations. Parity constraints are probably very helpful in this case when constructing the minus components, but then those relations can no longer be used to help constrain the coefficients in the LF-Hamiltonian.

Another class of potentially useful operators consists of the scalar and pseudoscalar densities $`\overline{\psi }\psi `$ and $`\overline{\psi }\gamma _5\psi `$. Obviously, if parity is conserved then at most one of the two couplings

$`f_S`$ $`=`$ $`0|\overline{\psi }\psi |n`$
$`f_P`$ $`=`$ $`0|\overline{\psi }\gamma _5\psi |n`$ (12)
can be nonzero at the same time, since a state $`|n`$ cannot be both scalar and pseudoscalar. What restricts the usefulness of this criterion is the fact that the same ‘selection rule’ also follows from charge conjugation invariance ($`\overline{\psi }`$ and $`\overline{\psi }\gamma _5\psi `$ have opposite charge parity!) which is a manifest symmetry on the LF. Therefore, only for theories with two or more flavors, where one can consider operators such as $`\overline{u}s`$ and $`\overline{u}\gamma _5s`$, does one obtain parity constraints that are not protected by charge conjugation invariance. But even then one may have to face the issue of how to define these operators. Parity may of course be used in this process, but then it has again been ‘used up’ and one can no longer use these selection rules to test the Hamiltonian.

In summary, many parity relations are probably very useful to determine the LF representation of the operators involved, but may not be very useful to determine the Hamiltonian.

### 3.3 A useful observable to test parity

We have seen above that relations between $`j^+`$ and $`j^{}`$ are often protected by vector current conservation and are therefore not very useful as parity tests. A much more useful parity test can actually be obtained by considering the ‘plus’-component only: Let us now consider a vector form factor,

$$p^{},n|j^\mu |p,m\stackrel{!}{=}\epsilon ^{\mu \nu }q_\nu F_{mn}(q^2),$$
(13)
where $`q=p^{}p`$, between states of opposite parity. When writing the r.h.s. in terms of one invariant form factor, use was made of both vector current conservation and parity invariance. A term proportional to $`p^\mu +p_{}^{}{}_{}{}^{\mu }`$ would also satisfy current conservation, but has the wrong parity. A term proportional to $`\epsilon ^{\mu \nu }\left(p_\nu +p_\nu ^{}\right)`$ has the right parity, but is not conserved and a term proportional to $`q^\mu `$ is both not conserved and violates parity. Other vectors do not exist for this example. The Lorentz structure in Eq. (13) has nontrivial implications even if we consider only the “plus” component, yielding

$$\frac{1}{q^+}p^{},n|j^+|p,m=F_{mn}(q^2).$$
(14)
That this equation implies nontrivial constraints can be seen as follows: as a function of the longitudinal momentum transfer fraction $`xq^+/p^+`$, the invariant momentum transfer reads ($`M_m^2`$ and $`M_n^2`$ are the invariant masses of the in and outgoing meson)

$$q^2=x\left(M_m^2\frac{M_n^2}{1x}\right).$$
(15)
This quadratic equation has in general, for a given value of $`q^2`$, two solutions for $`x`$ which physically correspond to hitting the meson from the left and right respectively. The important point is that it is not manifestly true that these two values of $`x`$ in Eq. (14) will give the same value for the form factor, which makes this an excellent parity test.

In Ref. , the coupling as well as the physical masses of both the fermion and the lightest boson where kept fixed, while the “vertex mass” was tuned (note that this required re-adjusting the bare kinetic masses). Figure 1 shows a typical example, where the calculation of the form factor was repeated for three values of the DLCQ parameter K (24, 32 and 40) in order to make sure that numerical approximations did not introduce parity violating artifacts.

For a given physical mass and boson-fermion coupling, there exists a “magic value” of the vertex mass and only for this value one finds that the parity condition (14) is satisfied over the whole range of $`q^2`$ considered. This provides a strong self-consistency check, since there is only one free parameter, but the parity condition is not just one condition but a condition for every single value of $`q^2`$ (i.e. an infinite number of conditions). In other words, keeping the vertex mass independent from the kinetic mass is not only necessary, but also seems sufficient in order to properly renormalize Yukawa1+1.

In the above calculation, it was sufficient to work with a vertex mass that was just a constant. However, depending on the interactions and the cutoffs employed, it ma be necessary to introduce counter-term functions. But even in cases where counter-term functions need to be introduced, parity constraints may be very helpful in determining the coefficient-functions for those more complex effective LF Hamiltonians non-perturbatively.

## References


text:
# Weak Gravitational Lensing Bispectrum

## 1. Introduction

Weak gravitational lensing of faint galaxies probes the distribution of matter along the line of sight. Lensing by large-scale structure (LSS) induces correlation in the galaxy ellipticities at the percent level (e.g., Blandford et al  1991; Miralda-Escudé  1991; Kaiser  1992). Though challenging to measure, these correlations provide important cosmological information that is complementary to that supplied by the cosmic microwave background and potentially as precise (e.g., Jain & Seljak  1997; Bernardeau et al  1997; Kaiser  1998; Schneider et al  1998; Hu & Tegmark  1999; Cooray  1999; Van Waerbeke et al  1999; see Bartelmann & Schneider  2000 for a recent review). Indeed several recent studies have provided the first clear evidence for weak lensing in so-called blank fields (e.g., Van Waerbeke et al  2000; Bacon et al  2000; Wittman et al  2000; Kaiser et al  2000), though more work is clearly needed to understand even the statistical errors (e.g. Cooray et al  2000b).

Given that weak gravitational lensing results from the projected mass distribution, the statistical properties of weak lensing convergence reflect those of the dark matter. Non-linearities in the mass distribution induce non-Gaussianity in the convergence distribution. With the growing observational and theoretical interest in weak gravitational lensing, statistics such as the skewness have been suggested as probes of cosmological parameters and the non-linear evolution of large scale structure (e.g., Bernardeau et al  1997; Jain et al  2000; Hui  1999; Munshi & Jain  1999; Van Waerbeke et al  1999).

Here, we extend previous studies by considering the full convergence bispectrum, the Fourier space analog of three-point function. The bispectrum contains all the information present at the three point level, whereas conventional statistics, such as skewness, do not. The calculation of the convergence bispectrum requires detailed knowledge of the dark matter density bispectrum, which can be obtain analytically through perturbation theory (e.g., Bernardeau et al  1997) or numerically through simulations (e.g., Jain et al  2000; White & Hu  1999). Perturbation theory, however, is not applicable at all scales of interest, while numerical simulations are limited by computational expense to a handful of realizations of cosmological models with modest dynamical range. Here, we use a new approach to obtain the density field bispectrum analytically by describing the underlying three point correlations as due to contributions from (and correlations between) individual dark matter halos.

Techniques for studying the dark matter density field through halo contributions have recently been developed (Seljak  2000; Ma & Fry  2000b; Scoccimarro et al.  2000) and applied to two-point lensing statistics (Cooray et al  2000b). The critical ingredients are: the Press-Schechter formalism (PS; Press & Schechter  1974) for the mass function; the NFW profile of Navarro et al  (1996), and the halo bias model of Mo et al.  (1997). The dark matter halo approach provides a physically motivated method to calculate the bispectrum. By calibrating the halo parameters with N-body simulations, it can be made accurate across the scales of interest. Since lensing probes scales ranging from linear to deeply non-linear, this is an important advantage over perturbation-theory calculations.

Throughout this paper, we will take $`\mathrm{\Lambda }`$CDM as our fiducial cosmology with parameters $`\mathrm{\Omega }_c=0.30`$ for the CDM density, $`\mathrm{\Omega }_b=0.05`$ for the baryon density, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$ for the cosmological constant, $`h=0.65`$ for the dimensionless Hubble constant and a scale invariant spectrum of primordial fluctuations, normalized to galaxy cluster abundances ($`\sigma _8=0.9`$ see Viana & Liddle  1999) and consistent with COBE (Bunn & White  1997). For the linear power spectrum, we take the fitting formula for the transfer function given in Eisenstein & Hu  (1999).

In §2, we review the dark matter halo approach to modeling the density field. In §3 we apply the formalism to the convergence power spectrum, skewness, and bispectrum. We summarize our results in §4.

## 2. Density Power Spectrum and Bispectrum

### 2.1. General Definitions

Underlying the halo approach is the assertion that dark matter halos of virial mass $`M`$ are locally biased tracers of density perturbations in the linear regime. In this case, functional relationship between the overdensity of halos and mass can be expanded in a Taylor series

$$\delta ^\mathrm{h}(𝐱,M;z)=b_1(M;z)\delta (𝐱;z)+\frac{1}{2}b_2(M;z)\delta ^2(𝐱;z)+\mathrm{}$$
(1)
Mo et al.  (1997) give the following analytic predictions for the bias parameters which agree well with simulations:

$$b_1(M;z)=1+\frac{\nu ^2(M;z)1}{\delta _c},$$
(2)
and

$`b_2(M;z)`$ $`=`$ $`{\displaystyle \frac{8}{21}}[b_1(M;z)1]+{\displaystyle \frac{\nu ^2(M;z)3}{\sigma ^2(M;z)}}.`$ (3)
Here $`\nu (M,z)=\delta _c/\sigma (M,z)`$, where $`\sigma (M,z)`$ is the rms fluctuation within a top-hat filter at the virial radius corresponding to mass $`M`$, and $`\delta _c`$ is the threshold overdensity of spherical collapse (see Henry  2000 for useful fitting functions).

Roughly speaking, the perturbative aspect of the clustering of the dark matter is described by the correlations between halos, whereas the nonlinear aspect is described by the correlations within halos, i.e. the halo profiles. We will consider the Fourier analogies of the 2 and 3 point correlations of the density field defined in the usual way

$$\delta ^{}(𝐤)\delta (𝐤^{})=(2\pi )^3\delta (𝐤𝐤^{})P^\mathrm{t}(k),$$
(4)
$$\delta (𝐤_1)\delta (𝐤_2)\delta (𝐤_3)=(2\pi )^3\delta (𝐤_1+𝐤_2+𝐤_3)B^\mathrm{t}(k_1,k_2,k_3).$$

Here and throughout, we occasionally suppress the redshift dependence where no confusion will arise.

As we shall see, these spectra are related to the linear density power spectrum $`P(k)`$ through the bias parameters and the normalized 3d Fourier transform of the halo density profile $`\rho (r,M)`$

$$y(k,M)=\frac{1}{M}_0^{r_v}𝑑r\mathrm{\hspace{0.17em}4}\pi r^2\rho (r,M)\frac{\mathrm{sin}(kr)}{kr}.$$
(5)
Note that $`y(k,M)`$ can be written as a combination of sine and cosine integrals for computational purposes and $`y(k,M)0`$ as $`k0`$.

It is convenient then to define a general integral over the halo mass function $`dn/dM`$,

$`I_\mu ^\beta (k_1,\mathrm{},k_\mu ;z)`$ $``$ $`{\displaystyle 𝑑M\left(\frac{M}{\rho _b}\right)^\mu \frac{dn}{dM}(M,z)b_\beta (M)}`$ (6)
$`\times y(k_1,M)\mathrm{}y(k_\mu ,M),`$
where $`b_01`$. We use the Press-Schechter (PS; Press & Schechter  1974) mass function to describe $`dn/dM`$. We take the minimum mass to be $`10^3`$ M while the maximum mass is varied to study the effect of massive halos on lensing convergence statistics. In general, masses above $`10^{16}`$ M do not contribute to low order statistics due to the exponential decrease in the number density of such massive halos.

### 2.2. Power Spectrum and Bispectrum

Following Seljak  (2000), we can decompose the density power spectrum, as a function of redshift, into contributions from single halos (shot noise or “Poisson” contributions),

$$P^{\mathrm{PP}}(k)=I_2^0(k,k),$$
(7)
and correlations between two halos,

$$P^{\mathrm{hh}}(k)=\left[I_1^1(k)\right]^2P(k),$$
(8)
such that

$$P^\mathrm{t}=P^{\mathrm{PP}}+P^{\mathrm{hh}}.$$
(9)
As $`k0`$, $`P^{\mathrm{hh}}P(k)`$.

Similarly, we decompose the bispectrum into terms involving one, two and three halos (see Scherrer & Bertschinger  1991; Ma & Fry  2000b):

$`B^\mathrm{t}`$ $`=`$ $`B^{\mathrm{PPP}}+B^{\mathrm{Phh}}+B^{\mathrm{hhh}},`$ (10)
where

$`B^{\mathrm{PPP}}(k_1,k_2,k_3)=I_3^0(k_1,k_2,k_3),`$ (11)
for single halo contributions,

$`B^{\mathrm{Phh}}(k_1,k_2,k_3)=I_2^1(k_1,k_2)I_1^0(k_3)P(k_3)+\mathrm{Perm}.`$ (12)
for double halo contributions, and

$`B^{\mathrm{hhh}}(k_1,k_2,k_3)`$ $`=`$ $`\left[2J(k_1,k_2,k_3)I_1^1(k_3)+I_1^2(k_3)\right]`$
$`\times I_1^1(k_1)I_1^1(k_2)P(k_1)P(k_2)+\mathrm{Perm}.`$
for triple halo contributions. Here the 2 permutations are $`k_3k_1`$, $`k_2`$. Second order perturbation theory tells us that (Fry  1984; Bouchet et al  1992; Kamionkowski & Buchalter  1999)

$`J(k_1,k_2,k_3)`$ $`=`$ $`1{\displaystyle \frac{2}{7}}\mathrm{\Omega }_m^{2/63}+\left({\displaystyle \frac{k_3^2k_1^2k_2^2}{2k_1k_2}}\right)^2`$ (14)
$`\times \left[{\displaystyle \frac{k_1^2+k_2^2}{k_3^2k_1^2k_2^2}}+{\displaystyle \frac{2}{7}}\mathrm{\Omega }_m^{2/63}\right].`$
As $`k0`$, $`B^{\mathrm{hhh}}B^{\mathrm{PT}}`$, where $`B^{\mathrm{PT}}`$ is the bispectrum predicted by second-order perturbation theory

$`B^{\mathrm{PT}}(k_1,k_2,k_3)`$ $`=`$ $`2J(k_1,k_2,k_3)P(k_1)P(k_2)+\mathrm{Perm}.,`$
with permutations following $`k_3k_1`$, $`k_2`$.

### 2.3. Halo Profiles

To apply the halo prescription for the power spectrum and bispectrum, we need to know the profile of the halos. We take the NFW profile (Navarro et al  1996) with a density distribution

$$\rho (r,M)=\frac{\rho _s}{(r/r_s)(1+r/r_s)^2}.$$
(16)
The density profile can be integrated and related to the total dark matter mass of the halo within $`r_v`$

$$M=4\pi \rho _sr_s^3\left[\mathrm{log}(1+c)\frac{c}{1+c}\right],$$
(17)
where the concentration, $`c`$, is defined as $`r_v/r_s`$. Choosing $`r_v`$ as the virial radius of the halo, spherical collapse tells us that $`M=4\pi r_v^3\mathrm{\Delta }(z)\rho _b/3`$, where $`\mathrm{\Delta }(z)`$ is the overdensity of collapse (see e.g. Henry  2000) and $`\rho _b`$ is the background matter density today. We use comoving coordinates throughout. By equating these two expressions, one can eliminate $`\rho _s`$ and describe the halo by its mass $`M`$ and concentration $`c`$.

Following Cooray et al  (2000b), we take the concentration of dark matter halos to be

$$c(M,z)=a(z)\left[\frac{M}{M_{}(z)}\right]^{b(z)},$$
(18)
where $`a(z)=10.3(1+z)^{0.3}`$ and $`b(z)=0.24(1+z)^{0.3}`$. Here $`M_{}(z)`$ is the non-linear mass scale at which the peak-height threshold, $`\nu (M,z)=1`$. The above concentration is chosen so that dark matter halos provide a reasonable match to the the non-linear density power spectrum as predicted by the Peacock & Dodds  (1996); it extends the treatment of Seljak  (2000) to the redshifts of interest for lensing. We caution the reader that eqn. (18) is only a good fit for the $`\mathrm{\Lambda }`$CDM model assumed.

### 2.4. Results

In Fig. 1(a-b), we show the density field power spectrum today ($`z=0`$), written such that $`\mathrm{\Delta }^2(k)=k^3P(k)/2\pi ^2`$ is the power per logarithmic interval in wavenumber. In Fig 1(a), we show individual contributions from the single and double halo terms and a comparison to the non-linear power spectrum as predicted by the Peacock & Dodds  (1996) fitting function. In Fig. 1(b), we show the dependence of density field power as a function of maximum mass used in the calculation.

Since the bispectrum generally scales as the square of the power spectrum, it is useful to define

$$\mathrm{\Delta }_{\mathrm{eq}}^2(k)\frac{k^3}{2\pi ^2}\sqrt{B(k,k,k)},$$
(19)
which represents equilateral triangle configurations, and its ratio to the power spectrum

$$Q_{\mathrm{eq}}(k)\frac{1}{3}\left[\frac{\mathrm{\Delta }_{\mathrm{eq}}^2(k)}{\mathrm{\Delta }^2(k)}\right]^2.$$
(20)
In second order perturbation theory,

$$Q_{\mathrm{eq}}^{\mathrm{PT}}=1\frac{3}{7}\mathrm{\Omega }_m^{2/63}$$
(21)
and under hyper-extended perturbation theory (HEPT; Scoccimarro & Frieman  1999),

$$Q_{\mathrm{eq}}^{\mathrm{HEPT}}(k)=\frac{42^{n(k)}}{1+2^{n(k)+1}},$$
(22)
which is claimed to be valid in the deeply nonlinear regime. Here, $`n(k)`$ is the linear power spectral index at $`k`$.

In Fig. 1(c-d), we show $`\mathrm{\Delta }_{\mathrm{eq}}^2(k)`$ separated into its various contributions (c) and as a function of maximum mass (d). Since the power spectra and equilateral bispectra share similar features, it is more instructive to examine $`Q_{\mathrm{eq}}(k)`$ (see Fig. 2). Here we also compare it with the second order perturbation theory (PT) and the HEPT prediction. In the halo prescription, $`Q_{\mathrm{eq}}`$ at $`k10k_{\mathrm{nonlin}}10h`$Mpc<sup>-1</sup> arises mainly from the single halo term. The HEPT prediction exceeds the halo prediction on larger scales and falls short on smaller scales.

### 2.5. Discussion

Even though the dark matter halo formalism provides a physically motivated means for calculating the statistics of the dark matter density field, there are several limitations of the approach that should be borne in mind when interpreting the results.

The approach assumes all halos to be spherical with a single profile shape. Any variations in the profile through halo mergers and resulting substructure can affect the power spectrum and higher order correlations. Also, real halos are not perfectly spherical which affects the configuration dependence of the bispectrum.

Furthermore, there are parameter degeneracies in the formalism that prevent a straightforward interpretation of observations in terms of halo properties. For example, one might think that the power spectrum and bispectrum can be used to measure any mean deviation from the assumed NFW profile form. However as pointed out by Seljak  (2000), changes in the slope of the inner profile can be compensated by changing the concentration as a function of mass; this degeneracy is also preserved in the bispectrum.

We do not expect these issues to affect our qualitative results. However, if this technique is to be used for precision studies of cosmological parameters, more work will be required in testing it quantitatively against simulations. Studies by Ma & Fry  (2000a) and Scoccimarro et al.  (2000) show that the bispectrum predictions of the halo formalism are in good agreement with simulations, at least when averaged over configurations. The replacement of individual halos found in numerical simulations with synthetic smooth halos with NFW profiles by Ma & Fry  (2000b) show that the smooth profiles can regenerate the measured power spectrum and bispectrum in simulations. This agreement, at least at scales less than 10$`k_{\mathrm{nonlin}}`$, suggests that mergers and substructures may not be important at such scales.

## 3. Convergence Power Spectrum and Bispectrum

### 3.1. Power Spectrum and Variance

The angular power spectrum of the convergence is defined in terms of the multipole moments $`\kappa _{lm}`$ as

$$\kappa _{lm}^{}\kappa _{l^{}m^{}}=C_l^\kappa \delta _{ll^{}}\delta _{mm^{}}.$$
(23)
$`C_l`$ is numerically equal to the flat-sky power spectrum in the flat sky limit. It is related to the dark matter power spectrum by (Kaiser  1992; 1998)

$$C_l^\kappa =𝑑r\frac{W(r)^2}{d_A^2}P^\mathrm{t}(\frac{l}{d_A};r),$$
(24)
where $`r`$ is the comoving distance and $`d_A`$ is the angular diameter distance. When all background sources are at a distance of $`r_s`$, the weight function becomes

$$W(r)=\frac{3}{2}\mathrm{\Omega }_m\frac{H_0^2}{c^2a}\frac{d_A(r)d_A(r_sr)}{d_A(r_s)};$$
(25)
for simplicity, we will assume $`r_s=r(z_s=1)`$ throughout. In deriving Eq. (24), we have used the Limber approximation (Limber  1954) by setting $`k=l/d_A`$ and the flat-sky approximation. In Cooray et al  (2000b), we used the projected mass of individual halos to construct the weak lensing power spectrum directly. The two approaches are essentially the same since the order in which the projection is taken does not matter.

In Fig. 3, we show the convergence power spectrum of the dark matter halos compared with that predicted by the Peacock & Dodds  (1996) power spectrum. The lensing power spectrum due to halos has the same behavior as the dark matter power spectrum. At large angles ($`l100`$), the correlations between halos dominate. The transition from linear to non-linear is at $`l500`$ where halos of mass similar to $`M_{}(z)`$ contribute. The single halo contributions start dominating at $`l>1000`$.

As shown in Fig. 3(b), and discussed in Cooray et al  (2000b), if there is a lack of massive halos in the observed fields convergence measurements will be biased low compared with the cosmic mean. The lack of massive halos affect the single halo contribution more than the halo-halo correlation term, thereby changing the shape of the total power spectrum in addition to decreasing the overall amplitude.

Similar statements apply to variance statistics (second moments) in real space. The variance of a map smoothed with a window is related to the power spectrum by

$$\kappa ^2(\sigma )=\frac{1}{4\pi }\underset{l}{}(2l+1)C_l^\kappa W_l^2(\sigma ).$$
(26)
where $`W_l`$ are the multipole moments (or Fourier transform in a flat-sky approximation) of the window. For simplicity, we will choose a window which is a two-dimensional top hat in real space with a window function in multipole space of $`W_l(\sigma )=2J_1(x)/x`$ with $`x=l\sigma `$.

In Fig. 4(a-b), we show the second moment as a function of smoothing scale $`\sigma `$. Here, we have considered angular scales ranging from 5 to 90, which are likely to be probed by ongoing and upcoming weak lensing experiments. As shown, most of the contribution to the second moment comes from the double halo correlation term and is mildly affected by a mass cut off.

### 3.2. Bispectrum and Skewness

The angular bispectrum of the convergence is defined as

$$\kappa _{l_1m_1}\kappa _{l_2m_2}\kappa _{l_3m_3}=\left(\begin{array}{ccc}l_1& l_2& l_3\\ m_1& m_2& m_3\end{array}\right)B_{l_1l_2l_3}^\kappa .$$
(27)
Extending our derivation of the SZ bispectrum in Cooray et al  (2000b), we can write the angular bispectrum of the convergence as

$`B_{l_1l_2l_3}^\kappa `$ $`=`$ $`\sqrt{{\displaystyle \frac{(2l_1+1)(2l_2+1)(2l_3+1)}{4\pi }}}\left(\begin{array}{ccc}l_1& l_2& l_3\\ 0& 0& 0\end{array}\right)`$ (31)
$`\times \left[{\displaystyle 𝑑r\frac{[W(r)]^3}{d_A^4}B^\mathrm{t}(\frac{l_1}{d_A},\frac{l_2}{d_A},\frac{l_3}{d_A};r)}\right].`$
The more familiar flat-sky bispectrum is simply the expression in brackets (Hu  2000). The basic properties of Wigner-3$`j`$ symbol introduced above can be found in Cooray et al  (2000a).

Similar to the density field bispectrum, we define

$$\mathrm{\Delta }_{\mathrm{eq}l}^2=\frac{l^2}{2\pi }\sqrt{B_{lll}^\kappa },$$
(32)
involving equilateral triangles in $`l`$-space.

In Fig. 3(c-d), we show $`\mathrm{\Delta }_{\mathrm{eq}l}^2`$. The general behavior of the lensing bispectrum can be understood through the individual contributions to the density field bispectrum: at small multipoles, the triple halo correlation term dominates, while at high multipoles, the single halo term dominates. The double halo term contributes at intermediate $`l`$’s corresponding to angular scales of a few tens of arcminutes.

In Fig. 5, we show the configuration dependence

$$R_{l_1l_2}^{l_3}=\frac{l_1l_2}{2\pi }\frac{\sqrt{B_{l_1l_2l_3}^\kappa }}{\mathrm{\Delta }_{\mathrm{eq}l}^2}$$
(33)
as a function of $`l_1`$ and $`l_2`$ when $`l_3=1000`$. The surface, and associated contour plot, shows the contribution to bispectrum from triangular configurations in $`l`$ space relative to that from the equilateral configuration. Due to the triangular conditions associated with $`l`$’s, only the upper triangular region of $`l_1`$-$`l_2`$ space contribute to the bispectrum. The symmetry about $`l_1=l_2`$ line is due to the intrinsic symmetry associated with the bispectrum. Though the weak lensing bispectrum peaks for equilateral configurations, the configuration dependence is weak.

The skewness is simply one, easily measured, aspect of the bispectrum. It is associated with the third moment of the smoothed map (c.f. eqn. )

$`\kappa ^3(\sigma )`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle \underset{l_1l_2l_3}{}}\sqrt{{\displaystyle \frac{(2l_1+1)(2l_2+1)(2l_3+1)}{4\pi }}}`$ (36)
$`\times \left(\begin{array}{ccc}l_1& l_2& l_3\\ 0& 0& 0\end{array}\right)B_{l_1l_2l_3}^\kappa W_{l_1}(\sigma )W_{l_2}(\sigma )W_{l_3}(\sigma ).`$
We then construct the skewness as

$$S_3(\sigma )=\frac{\kappa ^3(\sigma )}{\kappa ^2(\sigma )^2}.$$
(38)
The effect of the mass cut off is dramatic in the third moment. As shown in Fig 4(c-d), most of the contributions to the third moment come from the single halo term, with those involving halo correlations contributing significantly only at angular scales greater than $``$ 25. With a mass cut off, the total third moment decreases rapidly and is suppressed by more than three orders of magnitude when the maximum mass drops to $`10^{13}`$ M. The skewness only saturates when the maximum mass is raised to a few times $`10^{15}`$ M. Even though a small change in the maximum mass does not greatly change the convergence power spectrum (Fig. 3 of Cooray et al  2000b), the third moment, or the bispectrum, is strongly sensitive to the rarest or most massive dark matter halos.

In Fig. 6 we plot the skewness as a function of maximum mass, ranging from $`10^{11}`$ to $`10^{16}`$ M. Our total maximum skewness agrees with what is predicted by numerical particle mesh simulations (White & Hu  1999) and yields a value of $``$ 116 at 10. However, it is lower than predicted by HEPT arguments and simulations of Jain et al  (2000), which suggest a skewness of $``$ 140 at angular scales of 10. The skewness based on second-order PT is factor of $``$ 2 lower than the maximum skewness predicted by halo calculation. As shown, the PT skewness decreases slightly from angular scales of few arcmins to 90 and increases thereafter.

The effect of maximum mass on the skewness is interesting. When the maximum mass is decreased to $`10^{15}`$ M from the maximum mass value where skewness saturates ($`3\times 10^{15}`$ M), the skewness decreases from $``$ 116 to 98 at an angular scale of 10, though the convergence power spectrum only changes by less than few percent when the same change on the maximum mass used is made. When the maximum mass used in the calculation is $`10^{13}`$ M, the skewness at 10 is $`8`$, which is roughly a factor of 15 decrease in the skewness from the total.

The variation in skewness as a function of angular scale is due to the individual contribution to the second and third moments. The increase in the skewness at angular scales less than $``$ 30 is due to the single halo contributions for the third moment. The triple halo correlation terms dominate angular scales greater than 50, leading to a slight increase toward large angles, e.g. from $``$ 74 at 40 to $``$ 85 at 90. However, this increase is not present when the maximum mass used in the calculation is less than $`10^{14}`$ M. Even though mass cut off affects the single halo contributions more than the halo contribution, at such masses, the change in halo contribution with mass cut off prevents an increase in skewness at large angular scales.

The absence of rare and massive halos in observed fields will certainly bias the skewness measurement from the cosmological mean. One therefore needs to exercise caution in using the skewness to constrain cosmological models Hui  1999). In Cooray et al  (2000b), we suggested that lensing observations in a field of $``$ 30 deg<sup>2</sup> may be adequate for an unbiased measurement of the convergence power spectrum. For the skewness, observations within a similar area may be biased by as much as $``$ 25%. This is consistent with the sampling errors found in numerical simulations: 1$`\sigma `$ errors of 24% at $`10^{}`$ with a 36 deg<sup>2</sup> field (White & Hu  1999). To obtain the skewness within few percent of the total, one requires a fair sample of halos out to $`3\times 10^{15}`$ M, requiring observations of $``$ 1000 deg<sup>2</sup>, which is within the reach of upcoming lensing surveys involving wide-field cameras, such as the MEGACAM at Canada-France-Hawaii-Telescope (Boulade et al  1998), and proposed dedicated telescopes (e.g., Dark Matter Telescope; Tyson, private communication).

Still, this does not mean that non-Gaussianity measured in smaller fields will be useless. With this halo approach one can calculate the expected skewness if one knows that the most massive halos are not present in the observed fields. This knowledge may come from external information such as X-ray data and Sunyaev-Zel’dovich measurements or internally from the lensing data.

### 3.3. Related Statistics

The halo description in general allows one to test the effect of rare massive halos on any statistic related to the two and three point functions. In particular, it can be used to design more robust statistics.

Generalized three point statistics have been considered previously by Jain et al  (2000) following Nusser & Dekel  (1993) and Juszkiewicz et al  (1995). One such statistic is the $`\kappa |\kappa |_{\kappa >0}`$, which is expected to reduce the sampling variance from rare and massive halos (see, Jain et al  2000 for details). This statistic is proportional to $`\kappa ^3/\kappa ^2^{1/2}`$. In Fig. 7(a), we show this statistic as a function of maximum mass used in the calculation. We still find strong variations with changes to the maximum mass. Similar variations were also present in other statistics considered by Jain et al  (2000).

Consider instead the generalized statistic

$$S_3(\sigma ,m)=\kappa ^3/\kappa ^2^m$$
(39)
where $`m`$ is an arbitary index. We varied $`m`$ such that the effect of mass cuts are minimized on skewness. In Fig. 7(b), we show such an example with $`m=3.7`$. Here, the values are separated to two groups involving with most massive and rarest halos and another with halos of masses $`10^{14}`$ M or less. Though the values from the two groups agree with each other on small angular scales, they depart significantly above $`25^{}`$ reaching a difference of 2.5 at 80. Statistics involving such a high index $`m`$, weight the single halo contributions highly when the most massive halos are present, whereas they weight the halo correlation terms more strongly for $`M<10^{14}`$ M. To some extent this may be useful to identify the presence of rare halos in the observations.

However the consequence of using these generalized statistics is that one progressively loses their independence on the details of the cosmological model, e.g. the shape and amplitude of the underlying density power spectrum, as one departs from $`m=2`$, thereby contaminating the probe of dark matter and dark energy. Further work is necessary find the optimal trade off between robustness and cosmological independence of these and other generalized statistics.

Another observable statistic is the probability distribution function (pdf) of the convergence maps smoothed on the scale $`\sigma `$. This possibility has been recently studied by Jain & van Waerbeke  (1999), where the reconstruction of pdf using peak statistics were considered. Using the Edgeworth expansion to capture small deviations from Gaussianity, one can write the pdf of convergence to second order as

$`p(\kappa )`$ $`=`$ $`{\displaystyle \frac{1}{\sqrt{2\pi \kappa ^2(\sigma )}}}e^{\kappa (\sigma )^2/2\kappa ^2(\sigma )}`$
$`\times \left[1+{\displaystyle \frac{1}{6}}S_3(\sigma )\sqrt{\kappa ^2(\sigma )}H_3\left({\displaystyle \frac{\kappa (\sigma )}{\sqrt{\kappa ^2(\sigma )}}}\right)\right],`$
where $`H_3(x)=x^33x`$ is the third order Hermite polynomial (see, Juszkiewicz et al  1995 for details).

In Fig. 8, we show the pdf of convergence at 12 as a function of maximum mass used in the calculation. As shown, the greatest departures from Gaussianity begin to occur when the maximum mass included is greater than $`10^{14}`$ M. Given that we have only constructed the pdf using terms out to skewness, the presented pdfs should only be considered as approximate. With increasing non-Gaussian behavior, the approximated pdfs are likely to depart from this form especially in the tails. As studied in Jain & van Waerbeke  (1999), the measurement of the full pdf can potentially be used a a probe of cosmology. Its low order properties describe deviations from Gaussianity near the peak as opposed to the skewness which is more weighted to the tails.

## 4. Summary & Conclusions

We have presented an efficient method to calculate the non-Gaussian statistics of lensing convergence at the three point level based on a description of the underlying density field in terms of dark matter halos. The bispectrum contains all of the three point information, including the skewness. The prior attempts at calculating lensing bispectrum and skewness were limited by the accuracy of perturbative approximations and the dynamic range and sample variance of simulations.

Though the present technique provides a clear and an efficient method to calculate the statistics of the convergence field, it has its own shortcomings. Halos are not all spherical, which can to some extent affect the configuration dependence in moments higher than the two point level. Substructures due to mergers of halos can also introduce scatter. Though such effects unlikely to dominate our calculations, further work using numerical simulations will be necessary to determine to what extent present method can be used as a precise tool to study the higher order statistics associated with weak gravitational lensing.

The dark matter halo approach also allows one to study possible selection effects that may be present in weak lensing observations due to the presence or absence of rare massive halos in the small fields that are observed. We have shown that the weak lensing skewness is mostly due to the most massive and rarest dark matter halos in the universe. The effect of such halos is stronger at the three point level than the two point level. The absence of massive halos, with masses greater than $`10^{14}`$ M, leads to a strong decrease in skewness, suggesting that a straightforward use of measured skewness values as a test of cosmological models may not be appropriate unless prior observations are available on the distribution of masses in observed lensing fields.

One can correct for such biases using the halo approach, however. To implement such a correction in practice, further work will be needed to calibrate the technique precisely against simulations across a wide range of cosmologies. Efficient techniques to correct for mass biases both in the lensing power spectrum and bispectrum will be needed. Alternatively, this technique can be used to search for generalized three point statistics that are more robust to sampling issues. Given the great potential to study the dark matter distribution through weak lensing, this issues merit further study.

We acknowledge useful discussions with Lam Hui, Roman Scoccimaro, Uros Seljak and Ravi Sheth. We thank Jordi Miralda-Escudé for initial collaborative work that led to this paper. ARC acknowledges financial support from John Carsltrom and Don York. WH is supported by the Keck Foundation.


text:
# The Instability Strip of M3

### Acknowledgments.

This work was supported in part by OTKA T-024022 and T-030954.

## References

Bakos, G. Á., et al. 2000, in preparation

Castelli F., Gratton, R. G., & Kurucz, R. L. 1997, A&A, 318, 841

Jurcsik, J. & Kovács G. 1996, A&A, 312, 111

Kovács G. & Jurcsik, J. 1997, A&A, 322, 218

Stetson, P. B. 1994, PASP, 106, 250

Szabó, R., Csubry, Z., Kolláth, Z., & Buchler, R. 2000, in these proceedings, p. 374


text:
# Random spread on the family of small-world networks

## I Introduction

The small-world network exhibits unusual connection properties. On one hand it shows strong clustering like regular graphs and on the other hand it shows very small average shortest path between any two nodes like random graphs. Watts and Strogatz have proposed a simple model to describe small-world networks WattsStrogatz . The model gives a prescription to generate a one parameter family of graphs, ranging from the highly clustered (regular) graphs to the random graphs.

Various properties of this model have been studied Sag ; NewmanWatts ; Moukarzel ; MooreNewman . The spread and percolation properties investigated in Refs. Sag ; NewmanWatts ; Moukarzel ; MooreNewman , deal with the spread of information (disease) along the shortest path in the graph or the spread along the spanning tree.

In this paper, we study random walk on the family of small-world networks. Such a random walk corresponds to random spread of information on the network. In any realistic application of the spread on a graph, we expect the spread to be somewhere in between the two extremes viz. shortest path and the random walk. e.g. In Milgram’s experiment Milgram , which studies the connection properties of social networks, the path of a letter from a randomly chosen point to a fixed target is traced. The only condition imposed on the transfer of letter was that, the letter should be given to the person whom the sender knows by first name. The path followed by such a letter would have both random as well as shortest path elements in it. Another example is the path of internet protocol packet which follows a similar algorithm for forwarding the packet Stevens .

The determination of shortest path between two nodes typically requires $`O(n^3)`$ operations Aho , where $`n`$ is the number of nodes. This process becomes prohibitively expensive as $`n`$ becomes very large. Examples of such networks are social networks, telephone networks, internet etc. Another problem with the determination of shortest path is the incomplete knowledge of the network. Hence, it is clear that an alternate method of generating path (which need not be shortest) becomes necessary in these networks. Our analysis of random walk shows that average access time between two nodes goes as $`O(n)`$, for small-world geometry. It is benificial to consider a network with random routing or switching, particularly if it has small-world properties. Thus the random routing emerges to be both, practical and computationally cheaper method for large networks.

In Section II we discuss analytical and numerical results for the average access time of the random walk on one parameter family of graph ranging from regular case to the random case. We introduce an independent step approximation which allows us to get analytical expressions for the average access time. We discuss these results in Section III. It is found that the random walk results are similar to that of the shortest path results. Thus from the nature of the outcome of an experiment it may be difficult to conclude whether the spread was random or along the shortest path. An important consequence of this result can be in the routing and switching in very big networks. Random routing is a promising method, particularly if the network has small-world properties. Section IV summarizes the results.

## II Random walk on small-world graphs

The random walk on a graph is performed as follows: We start with a fixed node (say $`i`$) and at each step make a jump to a node connected to $`i`$ with uniform probability $`1/d(i)`$, where $`d(i)`$ is the degree of the node $`i`$, thus performing a random walk. Such a random walk gives a finite Markov chain Lovasz . One of the most important quantities of interest in a finite Markov chain is the average access time. Let $`D_{i,j}`$ be the average access time, defined as the first passage time to the node $`j`$ if the walk starts from the node $`i`$. We denote $`D_j=D_{0,j}`$.

We perform the random walks on the family of graphs generated by the algorithm given by Watts and Strogatz WattsStrogatz . The prescription gives a one parameter family of graphs which interpolates between the regular case and the random case. We refer to this family as the family of small-world networks or graphs. The regular graph (denoted by $`p=0`$) is a graph with $`n`$ vertices on a circle with each node connected to $`2k`$ nearest neighbors. The parameter $`k`$ is suitably chosen to keep the graph sparse but connected. The other elements of the family are obtained by random rewiring of each edge in the graph with probability $`p`$. It is seen that the small-world behavior is prominent around the parameter value $`p=0.01`$. i.e., when only $`1\%`$ edges are rewired. The case $`p=1`$ corresponds to the random case 111In the rewiring algorithm of Watts and Strogatz the random graph defined by $`p=1`$ is not truly random. However, as far as spread is concern the small-world network ($`p=0.01`$) itself is good enough approximation to the random network. Hence, $`p=1`$ case is adequate to describe the random behavior..

### II.1 The regular case $`(p=0)`$

Fig. 1 shows the results of the average access time in simulation of random walk on the regular graph with $`1000`$ nodes for several values of $`k`$. The average access time shows a linear behavior for small $`m`$ and it shows quadratic nature for large $`m`$ due to the circular topology. The lines are the analytic curves obtained as follows.

From the expectation values of conditional events, we can easily write the recursion relation for average access time for the walk starting from node $`i`$ as Feller

$`D_{i,m}={\displaystyle \frac{1}{2k}}{\displaystyle \underset{j=1}{\overset{k}{}}}(D_{i+j,m}+D_{ij,m})+1`$ (1)
This is a $`2k`$-th order difference equation for $`D_{i,m}`$. We first consider the case $`k=1`$. In this case the above recursion relation reduces to a quadratic form given by

$`D_{i,m}={\displaystyle \frac{1}{2}}(D_{i+1,m}+D_{i1,m})+1`$ (2)
The equation can be solved using standard techniques Feller  and the solution is given by

$`D_{i,m}=(mi)^2+A(mi)+B`$ (3)
where $`A`$ and $`B`$ are constants. The constants are determined using the boundary conditions $`D_{m,m}=0`$ and $`D_{mn,m}=0`$. Hence,

$`D_{i,m}=i^2+(2mn)im^2+mn`$
Without loss of generality, we assume that the walk starts from $`i=0`$. So the average access time for site $`m`$ starting from zero is

$`D_m=m^2+mn`$ (4)
Curve (a) in Fig. 1 shows both the analytical and numerical results for the case $`k=1`$. The linear and the quadratic nature is clearly seen.

For general value of $`k`$, it is not possible to solve the difference equation (1) exactly. The numerical simulations show that the nature of the curves for different $`k`$ is the same as for $`k=1`$. This suggests that there may be a scaling relation for a general $`k`$. Using numerical data fit we find that the following scaling relation fits the data reasonably well.

$`D_m^k(n)(1+\mu \mathrm{ln}(k))D_{\frac{m}{k}}^1({\displaystyle \frac{n}{k}})`$ (5)
where $`D_m^k`$ is the average access time from the site $`0`$ to site $`m`$ on the graph with $`2k`$ nearest neighbors, $`\mu =0.86`$ 222 We have also tried other types of scaling relations in particular scaling relation of the form $`D_m^k(n)k^\mu D_{\frac{m}{k}}^1(\frac{n}{k})`$. However, we find that the relation (5) gives better fit to the numerical data..

A possible explanation for such a scaling can be a quotienting of a graph by sub-graph of size $`k`$.

The Fig 2 explains the quotienting procedure clearly. The outer graph in Fig 2 is of size $`30`$ and the inner graph is the quotiented graph which has only nearest neighbor connections. While the walk on the quotiented graph is described by equation (4), the coefficient probably comes from the average time spent in each block. Note that the number of blocks that the walker has to pass is $`\frac{m}{k}`$.

From Fig. 1 we see that the scaling relation (5) shows an excellent matching with the numerical results for various values of $`k`$.

### II.2 The random case $`(p=1)`$

For the completely random case the access time becomes independent of $`m`$ and $`D_{i,j}=n1`$, $``$ $`i`$, $`j`$. This result can be obtained as follows.

We note that to calculate average access time one must consider an ensemble of graph for a given $`p`$. The average access time is obtained by first averaging over several realizations of random walk on a given graph and then it is averaged over various members of the ensemble of graphs. We now introduce an independent step approximation where we assume that the order of these two averages is interchanged. Thus in this approximation each step of the random walk is averaged over all the realizations of the graphs. This approximation is a kind of mean field approximation done in statistical mechanics.

Let $`P_{i,j}`$ be the probability of reaching site $`j`$ from site $`i`$ in one step. It is obvious that $`_iP_{i,j}=1`$. Using independent step approximation, we get

$`P_{i,j}={\displaystyle \frac{1}{n1}}`$
Thus the probability of reaching a site $`m`$ at any time step is $`\frac{1}{n1}`$. Hence, the probability of reaching site $`m`$ for the first time in $`t`$ time steps is

$`P_m(t)=\left(1{\displaystyle \frac{1}{n1}}\right)^{t1}{\displaystyle \frac{1}{n1}}`$ (6)
Thus $`D_m`$ is given by

$`D_m`$ $`=`$ $`{\displaystyle \underset{t}{}}tP_m(t)`$ (7)
$`=`$ $`n1`$
### II.3 The intermediate case $`(0<p<1)`$

Fig. 3 shows the results of average access time on the graph with $`1000`$ nodes and $`k=5`$. Two distinct behaviors of average access time can be identified from the Fig. 3. For small values of $`m`$, the behavior is similar to that of regular case, while for the larger values of $`m`$ the average access time saturates and behaves like that of a random graph (See Sec. II.2). As $`p`$ increases the saturation of average access time becomes more prominent and the crossover from regular to saturation behavior takes place at smaller and smaller values of $`m`$. For $`p=0.01`$, which corresponds to the small-world behavior WattsStrogatz , the average access time behavior is almost same as that of the random graph except for small values of $`m`$ of the order of $`k`$.

To obtain analytical estimates of average access time we again make use of the independent step approximation defined in Sec. II.2. In this approximation each step in the walk is taken with a probability $`(1p)`$ to nearest $`2k`$ sites and with probability $`p`$ to the remaining sites. In analogy with the $`p=0`$ case we write a recursion relation for $`D_{i,m}`$ (See Appendix).

$`D_{i,m}`$ $`=`$ $`(1p)\left[{\displaystyle \frac{1}{2k}}{\displaystyle \underset{l=0}{\overset{k}{}}}(D_{i+l,m}+D_{il,m})+1\right]`$ (8)
$`+p\left[{\displaystyle \frac{1}{n2k1}}{\displaystyle \underset{\stackrel{l=0}{\text{ }lik,\mathrm{},i+k}}{\overset{n1}{}}}D_{l,m}+1\right]`$
This is a $`2k`$-th order difference equation. As in the case of $`p=0`$, we first consider the case $`k=1`$. The recursion relation becomes

$`D_{i,m}`$ $`=`$ $`(1p)\left[{\displaystyle \frac{1}{2}}(D_{i+1,m}+D_{i1,m})+1\right]`$ (9)
$`+p\left[{\displaystyle \frac{1}{n3}}{\displaystyle \underset{\stackrel{l=0}{\text{ }ii,i1,i+1}}{\overset{n1}{}}}D_{l,m}+1\right]`$
which can be written as

$`D_{i,m}`$ $`=`$ $`(1p)\left[{\displaystyle \frac{1}{2}}(D_{i+1,m}+D_{i1,m})+1\right]`$
$`+p\left[{\displaystyle \frac{1}{n3}}{\displaystyle \underset{l=0}{\overset{n1}{}}}D_{l,m}+1{\displaystyle \frac{1}{n3}}(D_{i,m}+D_{i+1,m}+D_{i1,m})\right]`$
with little algebra, we finally get

$`D_{i,m}=\xi (D_{i+1,m}+D_{i1,m})+\zeta `$ (10)
where

$`\xi ={\displaystyle \frac{(n3)p(n1)}{2(n3+p)}}`$
and

$`\zeta ={\displaystyle \frac{p\underset{j=0}{\overset{n1}{}}D_{i,m}+n3}{n3+p}}`$
Note that the sum $`_{j=0}^{n1}D_{i,m}`$ is independent of the site index and can be treated as a constant to be determined self-consistently from the solution. We solve the different equation (10) using the standard methods Feller .

$`D_{i,m}=A\theta _+^{(mi)}+B\theta _{}^{(mi)}{\displaystyle \frac{\zeta }{2\xi 1}}`$
where

$`\theta _\pm ={\displaystyle \frac{1}{2\xi }}(1\pm \sqrt{14\xi ^2})`$
The constants $`A`$ and $`B`$ are determined using the boundary conditions $`D_{m,m}=0`$ and $`D_{mn,m}=0`$. The solution is given by

$`D_m={\displaystyle \frac{\zeta }{2\xi 1}}\left[{\displaystyle \frac{\theta _{}^m\theta _+^n\theta _+^m\theta _{}^n\theta _{}^m+\theta _+^m}{\theta _+^n\theta _{}^n}}1\right]`$ (11)
where without loss of generality we have put the starting point as $`i=0`$.

The sum $`_{m=0}^{n1}D_m`$ occurring in $`\zeta `$ is calculated by summing Eq. (11) for all values of $`m`$ and is given by

$`{\displaystyle \underset{m=0}{\overset{n1}{}}}D_m={\displaystyle \frac{(n3)}{p\left(1+n/_{m=1}^{n1}\left[\frac{\theta _{}^m\theta _+^n\theta _+^m\theta _{}^n\theta _{}^m+\theta _+^m}{\theta _+^n\theta _{}^n}1\right]\right)}}`$
Note that Eq. (8) and further analysis exhibit random behavior as $`p\frac{(n2k1)}{(n1)}`$, rather than $`p=1`$ (See Appendix).

For general value of $`k`$, we again write a scaling relation similar to that of the case $`p=0`$.

$`D_m^k(n,p)(1+\mu (p,k)\mathrm{ln}(k))D_{\frac{m}{k}}^1({\displaystyle \frac{n}{k}},p)`$ (12)
The Fig. 4 shows the match of the numerical data with the analytic expressions. The parameter $`\mu `$ has a weak dependence on $`k`$ for non-zero $`p`$ which we have neglected in further analysis.

Fig. 5 shows the behavior of $`\mu `$ for the various values of $`p`$. In Fig. 4 it is clearly seen that for small and large values of $`p`$ the match numerical and equation is quite good but for intermediate values there is a considerable deviation from the numerical results for small values of $`m`$. This fact can be understood as follows:

For small values of $`p`$, i.e., when $`p10^4`$, the graph is nearly regular and the blocks of size $`k`$ are nearly completely connected graphs, but as $`p`$ increases the probability of leaving the block randomly to far away point increases, giving rise to higher average access time for the nearby points in the block, than the analytical values. Again for high values of $`p`$ the expression has a good match because as far as the average access time is concern the completely connected block and random block behave in similar way.

## III Discussion

It is interesting to compare the limiting behaviors of intermediate cases with those of results for the regular and random cases. The Eq. (9) reduces to Eq. (2) as $`p0`$. However, the solution (Eq. (11)) does not smoothly reduce to Eq. (4). This is because of the degeneracy in the roots of quadratic indicial equation of Eq. (2) for $`p=0`$ which is lifted for non-zero $`p`$. Thus the solution changes from a polynomial to an exponential in site index as $`p`$ becomes non-zero. This explains the reason behind the saturation of the solution for non-zero $`p`$ for asymptotic $`m`$ which can not be obtained for $`p=0`$.

For the random case we again find a similar situation. As $`p\frac{(n3)}{(n1)}`$, which corresponds to the random case (see Appendix), the Eq. (9) reduces to the $`D_m=\zeta `$ and is consistent with the solution $`D_0=0`$ and $`D_m=n1`$, $`m0`$. However, the solution (Eq. (11)) does not smoothly reduce to a constant function because one of the roots ($`\theta _+`$) diverges.

Next, we consider the average access time for the diametrically opposite node i.e., $`G(p)=D_{\frac{n}{2}}(p)`$. This quantity is of interest as it should be correlated to the average cover time for the graph.

Fig. 6 shows the behavior of the access time of the diametrically opposite point normalized with that of $`p=0`$ case, i.e. $`G(p)/G(0)`$, as a function of $`p`$. The figure also shows the graph of characteristic length $`L(p)/L(0)`$ WattsStrogatz , which is the average shortest path between any two sites. Both the curves in the graph corresponding to the shortest path and the random walk show similar behavior. This observation has interesting consequences. The random spread on a small-world network considerably reduces the access time compared to the regular graph as in the case of the shortest path spread. This can be very useful in applications such as routing and switching where random routing is cheaper. The determination of shortest path is generally very expensive (the number of operations goes as $`n^3`$) and also in many cases the complete information of the network is not available. Examples of such networks are social networks, telephone networks, internet etc. where the number of nodes is very large. In such cases random routing will be more effective and cheaper, particularly if the graph has the small-world geometry. In this case the average access time goes as $`O(n)`$.

In the calculations we have introduced an independent step approximation, where we average over the various realizations of graphs at each time step. This is a kind of mean field approximation and has allowed us to obtain recursion relations for the average access time. We expect this approximation to be reasonably good for the random case. Even for the intermediate cases, the results obtained by this approximation are in good agreement with the numerical values.

For general values of $`k`$ the recursion relations for average access time can not be solved. However, the behavior of average access time shows an interesting scaling relation in $`k`$. As discussed in Sec. II.1 the scaling relation corresponds to a quotienting procedure where a graph of $`n`$ nodes with $`2k`$ connections is scaled to a graph of $`\frac{m}{k}`$ nodes with nearest neighbor connections and preserving the far-edges Sag . The scaling gives very good fit to the numerical data except for small values of $`m`$, when $`p0.001`$.

## IV Summary

In this paper, we have studied random walk on the family of small-world graphs. For the regular case the average access time shows linear behavior for small distances and shows quadratic nature for large distances. An interesting scaling relation in $`k`$ is observed for the average access time. For the random case the average access time is $`(n1)`$ and is independent of distance and $`k`$. An independent step approximation has enabled us to get the analytical result for the average access time. The same approximation allows us to write a recursion relation for the intermediate values of $`p`$. For intermediate cases the average access time shows a crossover from regular to random behavior with increasing distances.

The normalized average access time of the diametrically opposite nodes shows almost identical behavior as that of the characteristic length as a function $`p`$. This observation can be very important in several applications where the number of nodes in graph are very large or the complete information about the graph is not available. In these cases random routing or switching will be beneficial and cheaper.

We thank Prof. V. Balakrishnan, Prof. G. Baskaran, Dr. G. Menon and Dr. R. Ramanujam for fruitful discussions.

## V Appendix

### V.1 Derivation of equation (9)

We assume that the probability of breaking of an edge is $`p`$. We identify that the contribution to the average access time of site $`j`$ comes from three types of events.

1. All the $`2k`$ neighbors of site $`j`$ are connected to $`j`$. The probability associated with this event is $`(1p)^{2k}`$.
2. Some of the neighbors are connected to $`j`$ (say $`2kr`$) and $`r`$ are connected to the far away sites. The probability associated with this event is $`p^r(1p)^{2kr}`$.
3. None of the neighbors of $`j`$ are connected to $`j`$. The probability associated with this event is $`p^{2k}`$.

By independent step approximation the degree of each node is $`2k`$. This enables up to write a recursion relation for average access time using the properties of expectation value of conditional probability Feller . We write

$`D_{i,m}`$ $`=`$ $`{\displaystyle \underset{r=0}{\overset{2k}{}}}\text{ }^{2k}C_rp^r(1p)^{2kr}\left\{{\displaystyle \frac{2kr1}{2k}}\left[{\displaystyle \underset{l=1}{\overset{k}{}}}\left(D_{il,m}+D_{i+l,m}\right)\right]+{\displaystyle \frac{r}{n2k1}}{\displaystyle \underset{\stackrel{l=0}{\text{ }lik,\mathrm{},i+k}}{\overset{n1}{}}}D_{l,m}\right\}+1`$ (13)
$`=`$ $`(1p)\left[{\displaystyle \frac{1}{2k}}{\displaystyle \underset{l=0}{\overset{k}{}}}(D_{i+l,m}+D_{il,m})+1\right]+p\left[{\displaystyle \frac{1}{n2k1}}{\displaystyle \underset{\stackrel{l=0}{\text{ }lik,\mathrm{},i+k}}{\overset{n1}{}}}D_{l,m}+1\right]`$ (14)
In the above calculation we have assumed that the edge broken with probability $`p`$ is not rewired to one of the $`2k`$ nearest neighbors. Due to this the random graph limit corresponds to $`p=\frac{(n2k1)}{(n1)}`$ instead of $`p=1`$. This can be seen by equating the probability of a connection to a nearest neighbor to the probability of a connecting to any of the other sites

$`{\displaystyle \frac{(1p)}{2k}}={\displaystyle \frac{p}{n2k1}}.`$