text:
# Matter Power Spectrum from the Lyman-Alpha Forest: Myth or Reality?

## 1 Introduction

The Lyman-alpha forest is perceived as being simple, and in this simplicity promises to become a key to many unsolved problems in cosmology. In the most recent demonstration of the power of simplicity, Croft et al.  presented a measurement of the matter linear power spectrum on scales which are nonlinear today – a measurement not easily reproducible by other means, such as galaxy surveys at low redshift.

Not only was the phenomenon that they modeled simple, but also the method that they used to recover the linear power spectrum was simple – perhaps excessively simple. For example, Croft et al.  considered only one prior cosmological model in deriving the correction function that translates the observed flux power spectrum into a linear power spectrum.

In this paper our purpose is to check whether their method is reliable by sampling over a range of prior cosmological models. We also attempt to estimate systematic errors and their covariance. We begin in §2 by making sure that we are able to reproduce Croft et al.  results. In §3 we give our results. We close in §4 with an optimistic conclusion.

## 2 Method

Since our goal is to evaluate the accuracy of the Croft et al.  method, we apply their method to a range of cosmological models. We adopt the same linear transfer function, but we vary the amplitude and the tilt, and we experiment with different geometries and different Hubble constants. In §3.2 we also report the effect of varying the linear power spectrum over narrow intervals of wavenumber. We use a standard Particle-Mesh method to simulate the distribution the matter at $`z=2.72`$. Following Croft et al.  we assume that the underlying baryon density follows that of the dark matter. At small scales the baryon density is smoothed out compared to the dark matter density, but the scale at which this occurs, the “filtering” scale, is much smaller than the range of scales of interest here  (we discuss this in detail in § 3.3.2). All our simulations have a box size of $`20h^1\text{ Mpc}`$ (in comoving coordinates) and $`256^3`$ particles on a mesh of the same size. The gravitational force is calculated with a Green function that includes the Optimal Antialising filter , by summing over three Brilloun zones.

In calculating the synthetic Lyman-alpha absorption spectra we follow the the Croft et al.  procedure in precise detail: we compute the synthetic spectra for 1000 lines of sight which cross the computational box at random angles (to avoid bias introduced for lines of sight parallel to the box sides) for an assumed value for the spatial homogeneous photoionization rate $`\mathrm{\Gamma }`$. We then compute the mean opacity $`\overline{\tau }(\mathrm{\Gamma })`$ and find the value of $`\mathrm{\Gamma }`$ that gives the desired value for the mean opacity ($`\overline{\tau }=0.349`$ except as described in §3.1.3). We then compute the 1D flux power spectra of the normalized flux for each line of sight using a Lomb periodogram and average the results over 1000 lines of sight. We have verified that this number – 1000 lines of sight – is fully sufficient to obtain numerical results with the precision well below the observational error.

Before we can test for systematic errors in the Croft et al.  method, we must first check that we can reproduce their results accurately.

In order to reduce the uncertainty in the mean flux power spectrum, an average of several different random realizations is usually used. Normally, this procedure requires a large number of realizations, because the variance decreases only as the square root of the number of realizations. Instead, we adopt a slightly different approach. Out of all random realizations, we first choose “good” ones, i.e. those that give a measured linear power spectrum close to the input linear power spectrum. We choose 3 “best” realizations out of 100 random ones, and an average over those three recovers the input power spectrum at least as well as a plain average over 100 truly random realizations.

Figure 1 shows the flux power spectrum $`P_F(k)`$ (shown as $`\mathrm{\Delta }^2(k)k^3P(k)/2\pi ^2`$) averaged over 3 and over 12 best realizations (equivalent of 100 and 400 random realizations). The difference between the two curves provides an estimate of about 1/2 of the uncertainty in the mean flux power spectrum (because averaging over 12 realizations gives half the variance of averaging over just 3 realizations). We notice that this difference is significantly smaller than the statistical error-bars of the observational data, indicating that our flux power spectrum is calculated with sufficient precision.

Hereafter we use three “good” realizations per normalization for each cosmological model. Thus, each model requires at least six simulations (two different normalizations above and below the best-fit value), and more if our initial choice for the two normalizations does not bracket the best-fit value.

As evidence that our method reproduces the Croft et al.  results, we show in Figure 2 their fiducial model (courtesy Rupert Croft) and the flux power spectrum from our simulation of exactly the same cosmological model. Notice that there exist a difference on large scales, which is most likely due to the smaller size of our computational box, but that it is smaller than the observational errors and so is unimportant.

## 3 Results

### 3.1 Systematic errors of the recovered linear power spectrum

In this paper our main concern is the effect that an assumption of a specific cosmological model makes on recovering the linear power spectrum from the flux power spectrum. Specifically, Croft et al.  assume that the flux power spectrum $`P_F(k)`$ is proportional to the linear power spectrum $`P_L(k)`$ at the same wavenumber

$$P_F(k)=b^2(k)P_L(k)$$
(1)
the bias factor $`b(k)`$ being independent of, or at least insensitive to, the power spectrum $`P_L(k)`$.

If the “bias factor” $`b(k)`$ were independent of the linear power spectrum, this would be a direct analogy to the biased galaxy linear power spectrum. However, this is not the case \- the “bias factor” $`b(k)`$ depends on the amplitude of the linear power spectrum $`P_L(k)`$, so that the relation (1) is not a linear relation. Correctly, equation (1) should be written as

$$P_F(k)=b^2[k,P_L]P_L(k),$$
(2)
where we use the square brackets to underline that $`b`$ is actually an operator acting on the linear power spectrum $`P_L(k)`$. Thus, the recovered “observed” linear power spectrum,

$$P_L^{\mathrm{obs}}(k)=\frac{P_F^{\mathrm{obs}}(k)}{b^2[k,P_L]},$$
(3)
actually depends on the assumed linear power spectrum, and this dependence appears as a systematic error which is not included in the estimate of the error bars by Croft et al. .

For a particular choice of $`P_L`$ that fits the observed flux power spectrum, relationship (2) can be linearized. For small deviations $`\mathrm{\Delta }P_L`$ about an assumed linear power spectrum $`P_L`$ we obtain

$$\mathrm{\Delta }P_F(k)=\underset{k^{}}{}b^2(k,k^{})\mathrm{\Delta }P_L(k^{}),$$
(4)
where the matrix $`b(k,k^{})`$ can be obtained by differentiating equation (2). Even in this linearized form this equation is more complicated that the Croft et al.  assumption (1).

Throughout this paper, when we refer to the “recovered linear power spectrum”, we always mean $`P_L^{\mathrm{obs}}(k)`$ obtained via equation (3). The sole purpose of this paper is to investigate how the dependence of $`b`$ on $`P_L`$, ignored by Croft et al. , affects the recovery procedure.

#### 3.1.1 Dependence on the Hubble constant or spectral curvature

An important part of the Croft et al.  procedure is the translation from the spatial units used in a simulation to velocity units of the observational data. This requires adopting a value for the Hubble constant at $`z=2.72`$ which is used to scale spatial scales to velocity space. Specifically, following Croft et al. , we adopt

$$v_{\mathrm{tot},i}=H_3\left(r_i+v_{\mathrm{pec},i}/H_{\mathrm{EdS},3}\right),$$
(5)
where $`r_i`$ is the physical position of a pixel $`i`$ in the synthetic spectrum, $`v_{\mathrm{tot},i}`$ is the velocity position of the same pixel, $`v_{\mathrm{pec},i}`$ is the peculiar velocity of a given pixel in the simulation (in the EdS cosmology), $`H_3`$ is the adopted value of the Hubble constant at $`z=2.72`$, and $`H_{\mathrm{EdS},3}=720h\text{ km/s/Mpc}`$ is the Hubble constant in the EdS cosmology at $`z=2.72`$. This procedure is identical to the one adopted in Croft et al. : we have applied our procedure to the simulation data kindly provided to us by Rupert Croft, and we obtained numerically indistinguishable results.

It is important to emphasize here that this rescaling of the Hubble constant only affects the transformation from position to velocity space. The primordial power spectrum of fluctuations as a function of position adopted in the simulations is not changed, but rescaling of the velocity units does result in the linear power spectra being different in velocity space.

It can be shown that for pure power-law power spectra all the dependence on the Hubble constant cancels out exactly. However, this is not true for a general power spectrum with varying slope. In order to test this dependence we have run a low density flat cosmological model in addition to the fiducial model. Both models have identical initial conditions, and both run in the matter-dominated universe ($`\mathrm{\Omega }_{m,3}=1`$ for both models). They differ only by the values of the Hubble constant $`H=H_0\sqrt{\mathrm{\Omega }_{m,0}}(1+z)^{3/2}`$, which are different by a factor $`\sqrt{0.4}=0.6`$. Croft et al.  did not use a Hubble constant appropriate to their cosmological model. Rather, they ran the simulation in the EdS cosmology, but adopted the Hubble constant from a flat cosmology with $`\mathrm{\Omega }_{m,0}=0.4`$. We also show the so rescaled model in Fig. 3 (which is the model also shown in Fig. 2). As can be seen from Figure 3, the three recovered power spectra differ by less than the random error of the measurement, and we thus conclude that the dependence on the Hubble constant (or, equivalently, on the rate of change in the local slope - but see § 3.2) is insignificant.

#### 3.1.2 Assumed effective equation of state

Croft et al.  point out that uncertainty in the “effective equation of state of the IGM”, i.e. in the relationship between the temperature and the density of the photoionized gas, is one of the dominant systematic errors in recovering the linear power spectrum. This relation is commonly parameterized as

$$TT_0(1+\delta )^{\gamma 1},$$
(6)
where $`T_0`$ and $`\gamma `$ are parameters and $`\delta `$ is the gas overdensity. Croft et al.  adopt the values of $`T_0=15,000\text{ K}`$ and $`\gamma =1.6`$. However, recent measurements of the effective equation of state at $`z3`$  suggest somewhat different values for these parameters: $`T_0=20,000\pm 3,000\text{ K}`$ and $`\gamma =1.3\pm 0.2`$ at $`z=2.72`$ (all three measurements taken together and weighted by their uncertainties). We have run the six models with values for $`T_0`$ and $`\gamma `$ in the observed range, and the results are shown in Figure 4. We find that the effect of the slope $`\gamma `$ on the recovered power spectrum is small, below the random uncertainties of the measurement, and so we do not attempt to model it, but regard it as part of the total systematic error quoted in Table 1. The amplitude $`T_0`$ however does affect the recovered linear power spectrum, which scales approximately inversely proportionally to $`T_0`$,

$$P_L(k;T_0)=P_L(k;20,000\text{ K})\left(\frac{20,000\text{ K}}{T_0}\right).$$
(7)
Rescaling the recovered power spectra this way reduces the rms difference between various models to about 1%.

This dependence is easy to understand: for a given density and velocity distribution and a given photoionization rate, higher temperature means broader absorption lines, which in turn implies a lower mean opacity. When the mean opacity is renormalized to the fiducial value 0.349, the photoionization rate is adjusted downward, and a given value of the flux now corresponds to a lower value of the density and thus the lower amplitude of the power spectrum (since readjustment of the photoionization rate scales all densities by the same factor).

Because the value $`T_0=20,000\text{ K}`$ is favored by observations, we adopt this value as our fiducial value in the rest of this paper.

#### 3.1.3 Assumed mean optical depth

Another possible systematic error discussed by Croft et al.  is the assumed value for the mean optical depth - the value used to normalize synthetic absorption spectra. They assumed a value of $`\tau =0.349`$ based on Press, Rybicki, & Schneider  (hereafter PRS) value. The PRS measurement actually gives $`\tau _{\mathrm{PRS}}=0.349_{0.034}^{+0.051}`$. A more recent measurement by McDonald et al.  (hereafter MCD) gives a somewhat lower value, $`\tau _{\mathrm{MCD}}=0.285\pm 0.025`$. In order to test all the parameter range suggested by observations, we have run four models with the mean optical depth different from the fiducial value: $`\tau =0.400`$ (PRG+$`1\sigma `$), $`\tau =0.315`$ (PRG-$`1\sigma `$ and about MCD+$`1\sigma `$), $`\tau =0.285`$ (MCD value), and $`\tau =0.260`$ (MCD-$`1\sigma `$). The flux power spectra and the recovered linear power spectra for these four cases and the fiducial case are given in Figure 5. As one can see, the recovered power spectra scale with $`\tau `$. We find a somewhat weaker dependence than that reported by Croft et al. , mostly because we only consider a range of optical depth around the observed values. The following simple power-law scaling takes out almost all dependence of the recovered linear power spectrum on $`\tau `$ (the rms difference between the rescaled power spectra is less than 1%):

$$P_L(k;\tau )=P_L(k;0.349)\left(\frac{0.349}{\tau }\right)^{0.75}.$$
(8)
The physical origin of this dependence is also clear: lower $`\tau `$ means higher densities for the same value of the flux.

#### 3.1.4 Slope of the prior power spectrum

Another possible source of systematic error is the slope of the prior linear power spectrum $`P_L(k)`$. We test four models, with $`n=1.0`$, $`n=0.95`$ (the Croft et al. model), $`n=0.9`$, and $`n=0.7`$. Figure 6 shows the flux and recovered linear variance for these models. The $`n=0.7`$ model does not fit the flux power spectrum as well as other models, but the fit is still statistically acceptable ($`\chi ^2=11`$ with 11 degrees of freedom). In all cases however, the recovered linear power spectrum agrees well with the fiducial model, indicating that the method is insensitive to quite significant variations in the slope of the prior power spectrum. In other words, the effective bias factor $`b(k)`$ in equation (1) appears similar for all the model spectra, an encouraging result.

#### 3.1.5 Dependence on the matter density

Similarly, the dependence on the matter density parameter at z=2.72, $`\mathrm{\Omega }_{m,3}`$, was not tested by Croft et al. . This may not be an issue after all, because for any reasonable flat model with a cosmological constant, $`\mathrm{\Omega }_{m,3}=1`$ to better than 5%. Nevertheless, we have tested the matter density dependence. We have run three open models with $`n=0.7`$ and $`\mathrm{\Omega }_{m,0}=0.45`$ ($`\mathrm{\Omega }_{m,3}=0.75`$), $`\mathrm{\Omega }_{m,0}=0.2`$ ($`\mathrm{\Omega }_{m,3}=0.5`$), and $`\mathrm{\Omega }_{m,0}=0.1`$ ($`\mathrm{\Omega }_{m,3}=0.3`$) respectively, and the results of these simulations are presented in Figure 7. One can see that the recovered linear power spectrum does depend on the assumed value of the density parameter at $`z=2.72`$, even if all three models fit the observed flux power spectrum.

The dependence on the density parameter cannot be rigorously derived analytically, but the form found empirically can be understood approximately on the basis of arguments from linear theory. The Gunn-Peterson optical depth $`\tau `$  at every spatial point is inversely proportional to the $`du/dx`$ along the line of sight , where $`u`$ is the total velocity (Hubble flow plus peculiar velocity) and $`x`$ is the comoving distance. In the linear regime

$$\frac{du}{dx}\frac{\stackrel{}{u}}{\stackrel{}{x}}=aH\left(1+f\delta \right),$$
(9)
where $`fd\mathrm{ln}D_+/d\mathrm{ln}a\mathrm{\Omega }^{0.6}`$. The first term in the parenthesis is the Hubble flow, while the second one is due to peculiar velocities. In the linear regime the Hubble flow of course dominates, but at $`z=3`$ most of the Lyman-alpha forest is mildly nonlinear, with $`\mathrm{\Delta }_k^2k^3P_L(k)/(2\pi ^2)1`$, as can be seen from Fig. 7.

If we assume that different $`k`$ values are independent (which is not quite true, but good enough for our hand-waving arguments), we would expect that the power spectrum for the optical depth $`\tau `$ is approximately inversely proportional to $`(1+f\mathrm{\Delta }_k)^2`$. In other words, we hypothesize that the recovered linear power spectrum depends on cosmological parameters in the following way:

$$P_L^{\mathrm{obs}}(k)=P_L^{\mathrm{fct}}(k)\left(\frac{1+\mathrm{\Delta }_k}{1+f_3\mathrm{\Delta }_k}\right)^2,$$
(10)
where we will call $`P_L^{\mathrm{fct}}(k)`$ the “factorized” linear power spectrum, $`f_3`$ is the values of $`f`$ at $`z=2.72`$, and

$$\mathrm{\Delta }_k\left(\frac{k^3P_L^{\mathrm{obs}}(k)}{2\pi ^2}\right)^{1/2}$$
(11)
and thus depends on $`P_L^{\mathrm{obs}}(k)`$. In this form $`P_L^{\mathrm{fct}}(k)`$ should be nearly independent of the value of the matter density parameter at $`z=2.72`$ $`\mathrm{\Omega }_{m,3}`$.

Figure 8 now gives $`P_L^{\mathrm{fct}}(k)`$ for the four models shown in 7. As one can see, the factorization (10) does account for the main dependence of the recovered linear power spectrum on cosmological parameters, although some differences remain. It is however not convenient to use in practice because it is nonlinear and implicit for $`P_L^{\mathrm{obs}}(k)`$. But we can notice that for $`\mathrm{\Omega }_{m,3}`$ substantially different from 1, $`\mathrm{\Delta }_k`$ is about unity for the range of scales where the observational error-bars are small. It turns out that a comparable factorization can be obtained if we adopt a fixed value $`\mathrm{\Delta }_k=1.4`$ in equation (10),

$$P_L^{\mathrm{obs}}(k)=P_L^{\mathrm{fct}}(k)\left(\frac{2.4}{1+1.4f_3}\right)^2.$$
(12)
This fit is linear, so it is easier to use in joint parameter estimation, while it still recovers the density dependence to about 2.2% precision.

#### 3.1.6 The systematic error of the recovered linear power spectrum

We can now summarize our results. The recovered linear power spectrum can be factorized in the following way:

$$P_L^{\mathrm{obs}}(k)=P_L^{\mathrm{fct}}(k)Q_\mathrm{\Omega }Q_TQ_\tau ,$$
(13)
where

$`Q_\mathrm{\Omega }`$ $``$ $`\left({\displaystyle \frac{2.4}{1+1.4f_3}}\right)^2,`$
$`Q_T`$ $`=`$ $`20,000\text{ K}/T_0,`$
$`Q_\tau `$ $`=`$ $`(0.349/\tau )^{0.75},`$ (14)
and $`P_L^{\mathrm{fct}}(k)`$ is independent of anything else (just numbers) and is shown in Figure 9 and in Table 1 together with its “systematic” error-bars. We put the word “systematic” in quotes because we considered only a subsample of all possible cosmological models, and our sampling of the parameter space is by no means uniform. We thus suggest that errors quoted by us can be considered as an estimate for the systematic error in the Croft et al.  measurement, but the precise value of such an error and its covariance matrix may depend very much on the specifics of the prior.

The observational data  sets the value of $`Q_T`$ at

$$Q_T=1\pm 0.15$$
(15)
and the error on $`Q_T`$ is uncorrelated with the systematic error on $`P_L^{\mathrm{fct}}(k)`$.

With the mean optical depth the story is more complicated, because there is no sufficiently accurate measurement of it. The PRS data give $`Q_{\tau ,\mathrm{PRS}}=1.00_{0.08}^{+0.11}`$, while the MCD data give $`Q_{\tau ,\mathrm{MCD}}=1.18\pm 0.07`$. Combining the two values in quadrature (because the data sets they used are independent), we find

$$Q_\tau =1.11\pm 0.05.$$
(16)
And for the most conservative estimate, we can include the whole range of quoted numbers:

$$0.9<Q_\tau <1.25.$$
(17)
### 3.2 Band-Power Windows

Measurements of flux power sample power not at a single wavenumber, but rather over a finite band of wavenumbers. The band power windows $`b(k,k^{})^2`$ in equation (4) can be extracted by differentiating the flux power spectrum $`P_F(k)`$ with respect to the linear power spectrum $`P_L(k)`$. The shape of the band-powers emerges most clearly if they are scaled with the fiducial powers, so we define scaled band-power windows by

$$\frac{\mathrm{ln}P_F(k)}{\mathrm{ln}P_L(k^{})}=b^2(k,k^{})P_L(k)/P_F(k^{}).$$
(18)
In order to measure the effective band-power windows in the present case, we ran four additional models with the fiducial power spectrum of Croft et al.  reduced by 20% in four different wave-bands, corresponding to four $`k`$ values. Linear power spectra for these models are shown in Figure 10a, and their respective flux power spectra are shown in Fig. 10b. The difference between the different models is comparable to the residual uncertainty in the mean flux power spectrum (Fig. 1), and thus we are only able to measure the difference between the various models with the reduced power to an accuracy of about 30%. A more accurate determination of the band-power windows would require an implausibly large number of simulations.

As can be seen, the abrupt reductions in the linear power spectrum lead not to abrupt reductions in the flux power spectrum, but rather to broad depressions.

Figure 11 shows the inferred band-power windows themselves, the derivative of the flux power spectrum with respect to the input linear power spectrum at a given wavenumber, both for the case when the amplitudes of the input power spectra are kept fixed, and when the amplitudes are adjusted to fit the observed flux power spectrum. Evidently the effective band-power windows are quite broad.

Fig. 11 indicates that the band-power windows are broader at smaller scales. Figure 12 attempts to quantify this feature by showing the same derivatives as Fig. 11, except as a function of $`r=2\pi /k`$ rather than $`k`$. To bring out the similarity of shapes, the four band-power windows from Fig. 11 have been shifted horizontally so that their center points $`r_0=2\pi /k_0`$ coincide (the centers of the band-power windows $`k_0`$ do not necessarily fall on $`k^{}`$, but are always within one bin in $`k`$ space, from which we conclude that this discrepancy is most likely due to the finite $`k`$-space sampling), and vertically to the same maximum. The four band-power windows do not have exactly the same width in real space. Fitting to a Gaussian of width $`\sigma `$, we find that $`\sigma `$ as a function of $`k^{}`$ can be approximated as

$$\sigma (k^{})\frac{25\text{ km/s}}{\sqrt{k^{}\times 1\text{ km/s}}}.$$
(19)
In order to illustrate that it is indeed peculiar velocities that are responsible for the smoothing of the recovered linear power spectrum, we show in Fig. 12 the derivative of the flux power spectrum with respect to the prior linear power spectrum at $`k^{}=0.0096\text{ (km/s)}^1`$ (an equivalent of the long-dashed line) with all peculiar velocities set to zero. In this case the peak is much narrower, and its width is explained by finite sampling of the $`k`$-space alone.

It is important to underscore here that it is the power spectrum itself which is smoothed, and not the flux. Croft et al.  did consider the effect of smoothing the Lyman-alpha absorption spectrum (Fig. 10 in their paper), but the effect that we discuss here is different – it is a correlation of power between the neighboring $`k`$ values rather than a reduction in power on small scales.

This effect is easy to understand physically. Let us imagine that the matter power spectrum is a delta function, i.e. the power is only nonzero at a given value $`k_0`$. In the absence of peculiar velocities the three-dimensional flux power spectrum would be non-zero only in a narrow range of scales around $`k_0`$ (we cannot claim that this range is infinitesimally small because the relationship between the linear matter power spectrum and the Lyman-alpha flux power spectrum is nonlinear; however, the solid line in Fig. 12 does demonstrate that this range is less than the $`k`$-space sampling in Croft et al.). In the presence of peculiar velocities a given structure in physical space would appear shifted along the line of sight in redshift space (which is where the flux is measured), thus smearing power over a range of scales comparable to the 1D velocity dispersion.

### 3.3 Missing physics

In the previous subsections we have considered only one possible source of the systematic error, namely the choice of the prior linear power spectrum. Another possible uncertainty comes from the fact that no simulation includes all of the relevant physics. We can identify at least three major physical ingredients that are missed by pure PM simulations:

1. gas pressure;
2. inhomogeneities in the ionizing background;
3. shocks in the gas.

We discuss these effects in that order.

#### 3.3.1 Gas pressure

The role of gas pressure is well understood. In the linear regime it can simply be expressed by the filtering scale of baryonic perturbations, $`k_F`$:

$$P_{\mathrm{gas}}(k)P_{\mathrm{DM}}(k)e^{2k^2/k_F^2},$$
(20)
where $`k_F`$ is a complicated function of the entire thermal history of the universe , and equation (20) is a good approximation to the exact solution to the linear theory equations for $`k<0.8k_F`$. There are also strong reasons to believe that the filtering scale describes the effect of gas pressure even in the nonlinear regime \[2000b\].

Recent measurements of the thermal history of the universe  give $`k_F=(35\pm 5)h\text{ Mpc}^1`$, significantly smaller than the range of scales $`k3h\text{ Mpc}^1`$ measured from the Lyman-alpha forest power spectrum (the reason being that pressure filtering, always lagging behind the growth of the Jeans mass, is subdominant to the thermal broadening of the absorption spectrum). Thus, the pressure filtering produces only about 2% correction to the recovered linear power spectrum, and it is unmeasurable at the current time. This is, of course, in a full agreement with the conclusions of Croft et al. .

#### 3.3.2 Inhomogeneities in the ionizing background

Inhomogeneities in the ionizing background can come from two classes of sources: large scale inhomogeneities from quasars, and small scale inhomogeneities from galaxies. The former is known to be unimportant , but the effect of the latter can only be estimated from advanced numerical simulations which include effects of the radiative transfer. Here we use a simulation from Gnedin \[2000a\]. We follow the Croft et al.  procedure for the last output of this simulation at $`z=4`$, but we analyze it two different ways. In the first case we use the neutral hydrogen fraction from the simulation data (which reflects the non-uniformity of the ionizing radiation background) to compute a synthetic Lyman-alpha forest at $`z=4`$. In the second case we assume that the ionizing radiation is uniform and compute the neutral hydrogen fraction from ionization equilibrium. The value of the background is calculated by normalizing to the mean opacity in the simulation, so that in both cases the mean opacity is the same. In both cases we adopt the gas temperature as given in the simulation. After computing the flux power spectra in both cases, we find that they differ by less than 1%.

The Gnedin \[2000a\] simulation suffers from its small box size ($`4h^1\text{ Mpc}`$), and therefore a 1% difference is probably an underestimate, but it appears unlikely that it underestimates the role of inhomogeneities in the ionizing background by, say, a factor of 10. We therefore conclude that inhomogeneities in the ionizing background do not bias the measurement of the matter power spectrum from the Lyman-alpha forest.

Surprisingly, while the power spectra in the two calculations are very similar, the actual value of the volume averaged photoionization rate differs by 20%, being larger for the inhomogeneous ionizing background. Thus, simulations that assume homogeneous ionization and adjust the photoionization rate to fit the mean opacity of the forest actually underestimate the photoionization rate by at least 20%, or, equivalently, overestimate the baryon density by at least 10% (and probably by a significantly larger factor). This may explain larger numbers for the baryon density typically found in such simulations .

One of the possible concerns with this conclusion is that the Gnedin \[2000a\] simulation was terminated at $`z=4`$, whereas the Croft et al.  measure the linear matter power spectrum at $`z=3`$.

However, observationally we know that the following three quantities do not change by more than a factor of two between $`z=4`$ and $`z=3`$: the mean photoionization rate , the star formation rate , and the mean separation between star forming galaxies (deduced from constancy of the mean UV luminosity density and the star formation rate). This immediately implies that the proximity regions of star forming galaxies occupy approximately the same fraction of the total volume at $`z=4`$ and $`z=3`$, and thus expected fluctuations in the ionizing background are comparable at $`z=4`$ and $`z=3`$. We thus can be confident that our conclusions about the role of fluctuations in the ionizing background at $`z=4`$ also hold at $`z=3`$.

#### 3.3.3 Shocks and other mess

Obviously, non-gravitational effects like shocks can significantly affect the distribution of the Lyman-alpha forest. Fortunately, we have observational data  that suggest that non-gravitational effects are not significant in the forest. First, the IGM temperature at $`z3`$ is quite low, only $`20,000\text{ K}`$, and is entirely consistent with pure photoheating of the gas. Second, the narrowness of the cut-off in the Doppler parameter vs column density distribution of the Lyman-alpha lines suggests that the photoionization-induced “effective equation of state” has a small scatter, again fully consistent with the pure photoheating of the gas. Any realistic shock heating would destroy this tight correlation.

Thus, while we cannot rigorously prove that shocks and other non-gravitational effects can be neglected in modeling the Lyman-alpha forest, they would require a rather special fine-tuning which closely mimics the work of gravity and photoionization.

## 4 Conclusions

The two main drawbacks of the Croft et al.  paper are (a) that it misses the dependence of the recovered linear power spectrum on cosmological parameters (eq. ) and (b) that it misses the fact that the effective band-power windows of the flux power spectrum are broad. The broad band-powers must induce broad correlations between estimates of the flux power, and consequently also of the recovered matter power spectrum. But as the whole, we must acknowledge that the measurement that they present is remarkable – it currently offers the best measurement of the linear matter power spectrum over a range of scales that are nonlinear today.

The effective band-power windows can be approximated as Gaussians with a dispersion $`\mathrm{\Delta }k=k^{3/2}/(2\pi )\times 25\text{ (km/s)}^{1/2}`$ in Fourier space:

$$P_L^{\mathrm{REC}}(k)=\frac{1}{\mathrm{\Delta }k\sqrt{2\pi }}_0^{\mathrm{}}𝑑k^{}P_L^{\mathrm{TRUE}}(k^{})e^{{\displaystyle \frac{(kk^{})^2}{2\mathrm{\Delta }k^2}}},$$
(21)
where $`P_L^{\mathrm{REC}}(k)`$ is the recovered and $`P_L^{\mathrm{TRUE}}(k)`$ is the true underlying linear power spectra respectively. The smoothing is caused by peculiar velocities.

And while the resulting strong covariance between estimates of power at different wavenumbers will make careful comparison of the data and cosmological models more complicated than simple $`\chi ^2`$ fitting, this measurement will undoubtedly play an important role in securing accurate values of the cosmological parameters.

We also notice in passing that residual inhomogeneities in the ionizing background do not affect the measurement of the linear power spectrum, but may substantially bias estimates of the cosmic baryon density from the Lyman-alpha forest data.

We thank Rupert Croft and David Weinberg for valuable comments. We are also grateful to Rupert for providing the revised version of their paper prior to posting it to astro-ph. This work was partially supported by NASA ATP grant NAG5-10763 and by the National Computational Science Alliance under grant AST-960015N, and utilized the SGI/CRAY Origin 2000 array at the National Center for Supercomputing Applications (NCSA).


format: text

text:
# 1 Introduction

## 1 Introduction

The standard model (SM), as formulated for a flat four-dimensional space-time, can be easily generalized to the RS model . In the RS model, the four dimensional world around us arises from the compactification of a curved 5D geometry. The fifth dimension is a $`S_1/_2`$ orbifold and contains two 4D branes at the orbifold fixed points at $`y=0`$ and $`y=\pi R`$. The 5D geometry is a slice of AdS with the line element $`ds^2=dy^2+e^{2\sigma (y)}\eta _{\mu \nu }dx^\mu dx^\nu `$, where $`\sigma (y)=k|y|`$, $`k`$ measures the curvature along the fifth dimension, $`x`$ represents the familiar four dimensions and $`y`$ represents the fifth dimension. The brane at $`y=0`$ is called the Planck-brane and the brane at $`y=\pi R`$ is called the TeV-brane because the typical mass scale on each of the branes is of the order of the effective 4D Planck mass and a TeV respectively. This exponential hierarchy of energy scales is generated by the warp factor $`\mathrm{\Omega }=e^{\pi kR}`$, with $`kR11`$, and provides a new solution to the gauge hierarchy problem.

By compactifying the fifth dimension, a field living in the bulk, i.e. in the full 5D space-time, can be decomposed into an infinite number of 4D fields with different effective 4D masses using a method known as the Kaluza-Klein (KK) decomposition. Depending on the choice of $`k`$ and $`R`$, there is a large gap between the mass of the lightest KK mode of a bulk field and its next excited state. The mass of the first excited state is typically of the order of 10 TeV (since experimental constraints rule out excited states which are significantly lighter), while the mass of the ground state is constrained by experimental data , since such fields correspond to the ones we actually observe in accelerators.

Bulk fields allow one to address problems related to non-renormalizable operators within the RS framework, fermion masses and mixings , neutrino masses  and gauge coupling unification . With bulk gauge fields, the couplings and masses of the weak gauge bosons deviate from their SM values. In ref.  the modifications of couplings and masses were treated independently of each other. This is a good approximation if the SM fermions live close to one of the two branes. In this letter we present a combined analysis of the two effects in order to cover the case where the fermions are weakly localized or delocalized in the extra dimension as well. Weakly localized fermions are expected for instance, if the warped geometry also induces the fermion mass hierarchy .

We find that the lowest KK excitations of gauge bosons and fermions have to be heavier than about 10 TeV. There is no window for “light” KK states even for delocalized fermions, as was the case for the constraints derived in refs. . Hence, it is questionable if these excitations can be directly studied at the LHC. Indirect evidence for extra dimensions, e.g. rare processes such as $`n`$$`\overline{n}`$ oscillations or $`\mu e+\gamma `$ , are therefore very important. Nevertheless, one may speculate if the scenario presented here could account for the small deviations of the SM predictions from recent electroweak precision data . This might allow a lowering of the bound to around 7 TeV.

## 2 The Kaluza-Klein reduction

With some slight modifications, the SM is easily embedded in a warped geometry. All of the fields are assumed to live in the bulk, except the Higgs field, which is confined to the TeV-brane . (This can be considered as an approximation to the case where the Higgs field lives in the bulk but is concentrated exponentially around the TeV-brane. The ground state of a scalar field behaves as $`e^{4\sigma }`$ and for values of $`kR`$ around 11, this approximation is accurate up to a few permille .)

The masses of the weak gauge bosons are generated by spontaneous symmetry breaking arising from the Higgs mechanism, $`M^2(y)=a^2k\delta (y\pi R)`$, where $`a`$ is a dimensionless parameter which is determined by how strongly the gauge boson couples to the Higgs field. Since there is no large hierarchy between the weak and KK scales, the 5D mass term $`M^2`$ should be included in the KK reduction of the weak gauge bosons from the very beginning . The masses and 5D wave functions of a gauge field are obtained from its 5D equations of motion, $`\frac{1}{\sqrt{g}}_M(\sqrt{g}g^{MN}g^{RS}F_{NS})M^2(y)g^{RS}A_S=0`$. Here, the $`A_4=0`$ gauge is imposed and this is only possible because the Higgs field is localized on a brane. If the Higgs were to propagate in the bulk instead, it would not be possible to impose both the unitarity and the $`A_4=0`$ gauge simultaneously and we would have to contend with $`A_4`$ excitations as well. Using the method of separation of variables, a gauge field can be decomposed as follows:

$$A_\mu (x,y)=\frac{1}{\sqrt{2\pi R}}\underset{n=0}{\overset{\mathrm{}}{}}A_\mu ^{(n)}(x)f_n(y),$$
(1)
where $`A_\mu ^{(n)}`$ satisfies the field equation for a gauge boson of mass $`m_n`$. This expansion is only valid in the weak coupling limit for non-Abelian gauge fields. However, if we drop the requirement that the $`A_\mu ^{(n)}`$’s satisfy the 4D equations of motion for a Yang-Mills field, any field configuration can be decomposed according to eq. (1). The wave functions $`f_n`$ follow from $`(_y^22\sigma ^{}_yM^2+e^{2\sigma }m_n^2)f_n=0`$, with $`\sigma ^{}=d\sigma /dy`$. The $`f_n`$ are normalized by the condition $`\frac{1}{2\pi R}_{\pi R}^{\pi R}f_n(y)^2𝑑y=1`$. The solution is given by

$$f_n=\frac{e^\sigma (J_1(\frac{m_n}{k}e^\sigma )+b_nY_1(\frac{m_n}{k}e^\sigma ))}{N_n},$$
(2)
where $`N_n`$ is the normalization constant, and $`b_n`$ and $`m_n`$ are obtained by solving the following system of equations :

$`b(m,a^2)`$ $`=`$ $`{\displaystyle \frac{\frac{a^2}{2}J_1(\frac{m}{k})+\frac{m}{k}J_0(\frac{m}{k})}{\frac{a^2}{2}Y_1(\frac{m}{k})+\frac{m}{k}Y_0(\frac{m}{k})}},`$ (3)
$`b_n(e^{\pi kR}kx_n,0)`$ $`=`$ $`b(kx_n,a^2),`$
where $`x_n=e^{\pi kR}m_n/k`$. In fig. 1 we present the wave functions of the ground state and the first two excited states of the $`Z`$ boson.

Bulk fermions are described by the 5D equation of motion $`(g^{MN}\gamma _M(_M+\mathrm{\Gamma }_M)+m_\mathrm{\Psi })\mathrm{\Psi }=0`$, where $`m_\mathrm{\Psi }=c\sigma ^{}`$ and $`\mathrm{\Gamma }_M`$ is the spin connection in the tetrad formulation. The Dirac mass term $`m_\mathrm{\Psi }`$ takes on the functional form it has because of the parity restriction $`\mathrm{\Psi }(y)=\pm \gamma _5\mathrm{\Psi }(y)`$. For the same reason, the KK ground state of a fermion can only be either left-handed or right-handed. These ground states are identified with the fermion fields actually observed by experiments. Again, the 5D field is decomposed into a KK tower ,

$$\mathrm{\Psi }(x,y)=\frac{1}{\sqrt{2\pi R}}\underset{n=0}{\overset{\mathrm{}}{}}\mathrm{\Psi }^{(n)}(x)f_n^c(y),$$
(4)
where $`\mathrm{\Psi }^{(n)}(x)`$ obeys the equation of motion for a fermion of mass $`m_n`$,

$$f_0^c(y)=\frac{e^{(2c)\sigma }}{N_0},$$
(5)
and

$$f_n^c(y)=\frac{e^{\frac{5}{2}\sigma }(J_\alpha (\frac{m_n}{k}e^\sigma )+b_nY_\alpha (\frac{m}{k}e^\sigma ))}{N_n},$$
(6)
for $`n>0`$. The wave functions $`f_n^c`$ are normalized by

$$_{\pi R}^{\pi R}\frac{e^{3\sigma }f_n^c(y)^2}{2\pi kR}𝑑y=1.$$
(7)
Fig. 2 shows the wave function for the ground state fermions for different values of $`c`$. For $`c>`$1/2 ($`<`$1/2) the zero mode is localized towards the Planck-brane (TeV-brane). For $`c`$ close to 1/2 the zero mode is only weakly localized or even delocalized ($`c=1/2`$).

## 3 Measuring the weak mixing angle

Experimentally, there are several ways to measure the weak mixing angle, $`\theta _W`$. Although the equations of the SM remain unchanged in 5D, the effective 4D behavior of the fields is slightly modified. In particular, measurements of $`\theta _W`$ made on the assumption of a flat 4D space-time would result in a value different from the actual 5D weak mixing angle, $`\theta _5`$. The weak mixing angle can be calculated from the ratio of the masses of the $`W`$ and $`Z`$ bosons,

$$\mathrm{cos}^2\theta _1=\frac{m_W^2}{m_Z^2},$$
(8)
where the subscript 1 indicates that this is the first possible definition. Alternatively, $`\theta _W`$ can also be obtained from the corresponding ratio of the gauge coupling strengths. In the SM, both definitions are equivalent (at the tree level).

In the framework of warped models, deviations from standard 4D electroweak physics arise. The relationship between the (5D) gauge coupling strengths and the physical gauge boson mass is no longer linear because of eq. 3 . As a result, the $`W`$ mass gets shifted upward compared to the SM.444 Curiously, the new experimental value for $`m_W`$ is indeed slightly above the SM prediction .  Furthermore, the couplings of the weak gauge bosons to fermions are modified. In particular, they become dependent on the localization of the fermion within the extra dimension . The effective 4D coupling between the ground state of a fermion and the $`n^{th}`$ excited state of a gauge boson is given by

$$g^{(n)}=\frac{g^{(5)}}{(2\pi R)^{3/2}}_{\pi R}^{\pi R}e^{3\sigma }f_0^c(y)^2f_n(y)𝑑y,$$
(9)
where $`g^{(5)}`$ is the 5D gauge coupling.

In the approach presented above electroweak symmetry breaking is treated in the 5D framework. Corrections to weak gauge boson masses and gauge coupling strengths arise automatically in the KK reduction. One can interpret these results also as a mixing effect between a massless ($`y`$ independent) zero mode of a gauge field and its ($`y`$ dependent) KK states. The mixing between the weak gauge bosons and their KK excitations is on the order of $`m_W^2/M_{KK}^2`$. This explains why all deviations from SM physics are proportional $`1/M_{KK}^2`$. Note that somewhat similar results are obtained for the case of a flat extra dimension .

In the following, we perform a tree level comparison between the warped model with bulk fields and the ordinary SM. This procedure is justified because SM-like quantum corrections are essentially the same in both models. The tiny modifications (of order $`10^3`$) of the tree level physics will only cause sub-leading corrections in the loops. Equivalently, we can just as well work in the picture of a gauge field zero mode mixing with its KK states. In the treatment of the zero mode we then include the usual SM quantum corrections, while the mixing with the excited states is approximated at tree-level. In the literature this approach has been successfully used to investigate the impact of extra Z bosons on electroweak observables , which is very similar to what is discussed here. In the warped model (like in models with extra Z bosons), there are additional radiative corrections coming from loops involving KK states (extra Z bosons). However, their influence is small because of the large masses (about 10 TeV) of the excited states .

After electroweak symmetry breaking, the neutral gauge bosons in the bulk mix via the 5D weak mixing angle $`\mathrm{tan}\theta _5=g_1^{(5)}/g_2^{(5)}`$, where $`g_2^{(5)}`$ and $`g_1^{(5)}`$ denote the $`SU(2)`$ and $`U(1)`$ gauge coupling strengths in five dimensions respectively. The gauge coupling strength of the $`Z`$ boson in 5D is then given by $`g_Z^{(5)}=g_2^{(5)}/\mathrm{cos}\theta _5`$. Analogous relations hold for the $`W`$ boson and the photon. The effective 4D gauge couplings to fermions, $`g_W`$, $`g_Z`$ and $`g_\gamma `$ are obtained from the $`n=0`$ case of eq. 9, using the relevant 5D gauge coupling strengths. From these quantities, the weak mixing angle can be defined as

$$\mathrm{cos}^2\theta _2=\frac{g_W^2}{g_Z^2},$$
(10)
or

$$\mathrm{cos}^2\theta _3=1\frac{g_\gamma ^2}{g_W^2}.$$
(11)
Experimentally, these definitions can be regarded as comparisons of the strengths of the charged and neutral current, and of the electric and charged current respectively.

In the 4D SM, the values of the weak mixing angle inferred from eqs. 8, 10 and 11 agree, of course. Since in the warped model, the gauge coupling strengths and gauge bosons masses are obtained from eqs. 3 and 9, this is no longer the case. The differences between these angles, as measured by

$$\mathrm{\Delta }_{ij}=\mathrm{cos}^2\theta _i\mathrm{cos}^2\theta _j,$$
(12)
is too small to be measured, given the current range of uncertainty in the experimental data. However, upper bounds can be set on the $`\mathrm{\Delta }`$’s based on the range of uncertainty. This in turn provides important constraints on the parameters of the warped SM model. Rather than listing the permissible range of values for the internal parameters of the model, a more transparent way of expressing these constraints from an experimental perspective is to provide a lower bound for the masses of the first KK excitation of the weak gauge bosons.

## 4 Numerical analysis

The graphs in this section show the constraints on the mass of the first KK excitation for different values of $`k/M_5`$, where $`M_5`$ (the Planck mass in 5D) has been chosen to fit the observed value for the Planck mass in 4D, $`M_{\mathrm{Pl}}^2M_5^3/k`$ . In each of the graphs, the lower bounds for the mass of the excited $`Z`$ boson, $`M_{KK}M_Z^{(1)}`$ (the difference in mass between the excited $`Z`$ and $`W`$ bosons is insignificant) arising from each of the three constraints are plotted. The mass of the first KK state is adjusted in our model by varying the radius of the orbifold. In all the following graphs, we use the experimental constraints $`\mathrm{\Delta }_{12}1.2\times 10^3`$, $`\mathrm{\Delta }_{13}1.2\times 10^3`$ and $`\mathrm{\Delta }_{23}1.6\times 10^4`$, based on data from the Particle Data Group . The bound on $`\mathrm{\Delta }_{12,13}`$ is dominated by the experimental error of the $`W`$ mass, $`M_W=80.419\pm 0.056`$ GeV, which results in an uncertainty in $`\mathrm{cos}^2\theta _1`$ of $`1.1\times 10^3`$. For the weak mixing angle, we take $`\mathrm{sin}^2\theta _W=0.23117\pm 0.00016`$ with the error providing the constraint on $`\mathrm{\Delta }_{23}`$.

The numerical analysis is done in the following manner. We fix the values of $`kR`$, $`k/M_5`$ and $`c`$. From the equation for the KK gauge boson spectrum (3), we determine the gauge boson mass parameters, $`a^2(W)`$ and $`a^2(Z)`$ by fitting the ground state mass to the experimental $`W`$ and $`Z`$ masses. The masses of the KK excitations are then also fixed, as is the value of $`\mathrm{cos}^2\theta _1`$. The 5D weak mixing angle follows from $`\mathrm{cos}^2\theta _5=a^2(W)/a^2(Z)`$. With this information, we can calculate the ratios of the effective gauge coupling strengths (9) which enter the definitions of $`\mathrm{cos}^2\theta _{2,3}`$. The deviations from SM physics are then easily obtained from eq. (12). We stress again that for different values of the KK mass, the deviations scale as $`\mathrm{\Delta }_{ij}1/M_{KK}^2`$.

In fig. 3 we display the three different constraints on $`M_{KK}`$ as a function of the 5D fermion mass for $`k/M_5=1`$. For fermions located close to the Planck-brane, $`c>1/2`$, the strongest bound, $`M_{KK}\stackrel{>}{}11.1`$ TeV, arises from $`\mathrm{\Delta }_{12}`$, i.e. from the modification of the ratio of the weak gauge boson masses. The constraint from $`\mathrm{\Delta }_{23}`$, which only involves the effective gauge coupling strengths, gives a somewhat lower bound, $`M_{KK}\stackrel{>}{}7.8`$ TeV. This is because near the Planck-brane, the gauge boson wave function is almost flat and independent of the gauge boson mass . Therefore, $`\mathrm{\Delta }_{12}`$ and $`\mathrm{\Delta }_{13}`$ lead to almost identical results. For $`c\stackrel{<}{}1/2`$, the fermions become sensitive to the dip in the gauge boson wave function at the TeV-brane  and the gauge couplings deviate from their SM values. The constraint from $`\mathrm{\Delta }_{23}`$ takes over and pushes the value of $`M_{KK}`$ far above 10 TeV. For fermions strictly confined to the TeV-brane ($`c\mathrm{}`$), we find $`M_{KK}\stackrel{>}{}62`$ TeV. For comparison, we also include the lower bound on $`M_{KK}`$ from ref.  in fig. 3. It arises from the contribution of the excited gauge bosons to electroweak observables, like for example, effective 4-fermion interactions. Our bounds turn out to be much stronger for all values of $`c`$.555This justifies our neglect of loops containing KK excitations in the comparison with the SM. In particular, $`M_{KK}1`$ TeV is excluded even for fermions close to the Planck-brane. This is the main result of the present paper.

In fig. 4, we give the corresponding results for $`k/M_5=0.01`$. Qualitatively, the picture is unchanged, but the constraints on $`M_{KK}`$ are slightly looser by about 7 percent. For fermions localized close to the Planck-brane, the bound from $`\mathrm{\Delta }_{12}`$ on $`M_{KK}`$ is relaxed to 10.3 TeV. In general, the bound on $`M_{KK}`$ depends logarithmically on $`k/M_5`$, with the relative change given approximately by $`0.015\mathrm{ln}(k/M_5)`$. In fig. 5, this behavior is demonstrated for $`\mathrm{\Delta }_{12}`$.

So far we have only included a single 5D mass parameter $`c`$. If all SM fermions have a common location in the extra dimension, the constraint $`c\stackrel{<}{}0.3`$ applies. Otherwise the overlap of the fermion wave functions and the Higgs is insufficient to generate the observed top quark mass. Moreover, universal values of $`c\stackrel{>}{}1/2`$ are disfavored by deviations from universality . Requiring $`c\stackrel{<}{}0.3`$ for the universal $`c`$ parameter, the KK scale has to be larger than at least 25 TeV, making this scenario less attractive. However, in the SM there is no symmetry that requires the $`c`$ parameters of different fermion species to be degenerate. The $`c`$ parameters can be chosen in such a way as to generate the fermion masses and mixings without introducing hierarchies in the 5D Yukawa couplings . The leptons and light quarks then turn out to have $`c`$ parameters larger than 1/2. Since these particles are the ones used in experiments, the 11 TeV bound on $`M_{KK}`$ applies also for this case.

## 5 Discussion and summary

According to ref. , the LHC can probe $`M_{KK}`$ up to at most 7 TeV. From our results so far, it seems unlikely that KK excitations of gauge bosons can be produced directly at the LHC. Still, there may be interesting experimental signatures of the presented scenario, like for instance, exotic processes such as lepton flavor violation or $`n\overline{n}`$ oscillations .

One should keep in mind however, that the crucial input in deriving the 10 TeV constraint is the tolerated deviation in the $`W`$ mass. We took $`\mathrm{\Delta }M_W=56`$ MeV . The current experimental accuracy on the $`W`$ mass is about 34 MeV . If we were to take this value to be $`\mathrm{\Delta }M_W`$, the quantities $`\mathrm{\Delta }_{12,13}`$ would have to be smaller than $`0.86\times 10^3`$. The lower bound on $`M_{KK}`$ would increase to 13.1 TeV (for $`k=M_5`$). Most interestingly, however, recent data seem to favor a $`W`$ mass which is somewhat above the SM prediction. At the $`1\sigma `$ level, the SM prediction and the experimental result have no overlap . The central values differ by about 76 MeV for a Higgs mass of 115 GeV. Taking this value for $`\mathrm{\Delta }M_W`$ results in $`M_{KK}\stackrel{>}{}9.4`$ TeV (for $`k=M_5`$). For larger Higgs masses, the SM deviates even more from experiment. For a Higgs mass of 250 GeV, the SM prediction is about 130 MeV lower than the value found experimentally. This could be accounted for by KK gauge bosons with a mass of 7 TeV. Thus, being very optimistic, KK gauge boson excitations may be at the verge of being discovered at the LHC.

For $`M_{KK}7`$ TeV, the modifications to the effective gauge coupling strengths push $`\mathrm{\Delta }_{23}`$ close to its allowed value, as can be observed in fig. 3. Since $`\mathrm{\Delta }_{23}`$ is rather sensitive to the location of fermions, experiments with different quarks and/or leptons may actually result in different values for $`\mathrm{sin}^2\theta _W`$. This may be of some relevance to the $`3.6\sigma `$ discrepancy in measuring $`\mathrm{sin}^2\theta _{\mathrm{eff}}^{\mathrm{lept}}`$ using leptonic and hadronic asymmetries .

In more elaborate models using extra branes or non-trivial bulk mass profiles multi-localization of gauge bosons and fermions is possible. . Then exceptionally light KK states can arise which may escape the bounds we found in our analysis. To avoid conflict with experiment these states have to decouple from SM physics.

In conclusion, our calculations suggest that the KK excitations of W and Z within the RS model are expected to be rather heavy, of the order of 10 TeV or so. With some caveats, this may be brought down to about 7 TeV, which would make their discovery at the LHC much easier.

## Acknowledgements

S.J.H. thanks S. Heinemeyer for a useful discussion, the Bartol Research Institute for its hospitality during the crucial stages of this work and the Alexander von Humboldt Stiftung for a Feodor-Lynen fellowship. We thank Tony Gherghetta for useful comments on the draft and for bringing ref.  to our attention.

This work was supported in part by DOE Contract DE-FG02-91ER40626.


text:
# I Preliminaries

## I Preliminaries

The most important result of Einstein’s special relativity is probably the relationship between mass and energy (see Chapter 7 of Ref. 1). How it happened that the mass became so critical in this theory? The answer is very simple. The theory involves the square of the speed of light, $`c^2`$. What physical quantity incorporates the square of velocity? It is the energy, $`mv^2/2`$. Hence it is the mass that goes with $`v^2`$ and $`c^2`$. Thus $`m`$ is taken as an invariant under the Lorentz transformation. This becomes the rest mass $`m_0`$. The inertial mass $`m`$ follows to depend on the velocity.

What is the comparable physical quantity in cosmological special relativity? Certainly not the mass. Here we have $`\tau ^2`$, the square of the Hubble time constant. What physical quantity goes with the square of time? It is the acceleration a because $`𝐚t^2/2`$ describes the distance a particle makes at time $`t`$ when it is subject to acceleration a. Hence a should be taken as the invariant quantity under the cosmological transformation. But then the acceleration depends on the cosmic time just as the mass depends on the velocity.

In this paper these concepts are explored.

## II Velocity and acceleration four-vectors

We start our four-dimensional spacevelocity analysis by defining the velocity and acceleration (see Section 5.5 of Ref. 1 for the special relativistic treatment in spacetime).

The velocity four-vector of a particle in spacevelocity is defined as the dimensionless quantity

$$u^\mu =\frac{dx^\mu }{ds},$$
$`(1)`$
where $`\mu =0,1,2,3,`$ $`x^\mu =(x^0,x^1,x^2,x^3)=(\tau v,x,y,z)`$, and $`\tau `$ is Hubble’s time in the limit of zero gravity, a universal constant whose value is $`12.5\times 10^9`$ years . In flat spacevelocity one has for the line element,

$$\tau ^2dv^2\left(dx^2+dy^2+dz^2\right)=ds^2,$$
$`(2)`$
thus

$$\tau ^2\left(\frac{dv}{ds}\right)^2\left(1\frac{dx^2+dy^2+dz^2}{\tau ^2dv^2}\right)=1.$$
$`(3)`$
This gives

$$\tau ^2\left(\frac{dv}{ds}\right)^2\left(1\frac{t^2}{\tau ^2}\right)=1,$$
$`(4)`$
and therefore

$$\frac{dv}{ds}=\frac{1}{\tau \sqrt{1{\displaystyle \frac{t^2}{\tau ^2}}}}.$$
$`(5)`$
The velocity four-vector in spacevelocity can thus be expressed as

$$u^\mu =\frac{dx^\mu }{ds}=\frac{dx^\mu }{dv}\frac{dv}{ds}=\frac{1}{\tau \sqrt{1{\displaystyle \frac{t^2}{\tau ^2}}}}\frac{dx^\mu }{dv}.$$
$`(6)`$
The velocitylike component of $`u^\mu `$ is therefore given by

$$u^0=\gamma ,$$
$`(7)`$
whereas its spatial components

$$\text{u}=u^k=(u^1,u^2,u^3)$$
$`(8)`$
are given by

$$u^k=\frac{\gamma }{\tau }\frac{dx^k}{dv},$$
$`(9)`$
where

$$\gamma =\frac{1}{\sqrt{1{\displaystyle \frac{t^2}{\tau ^2}}}}.$$
$`(10)`$
It will be noted that, by Eq. (1),

$$u_\alpha u^\alpha =1,$$
$`(11)`$
namely, the length of $`u^\mu `$ is unity.

The acceleration four-vector of a particle in spacevelocity is defined by

$$\frac{du^\mu }{ds}=\frac{d^2x^\mu }{ds^2}.$$
$`(12)`$
By differentiating Eq. (11) we find that the acceleration four-vector satisfies the orthogonality condition

$$u_\alpha \frac{du^\alpha }{ds}=0.$$
$`(13)`$
The components of the acceleration four-vector of a particle in spacevelocity are, by Eqs. (5) and (7)–(9), then given by

$$\frac{du^0}{ds}=\frac{du^0}{dv}\frac{dv}{ds}=\frac{\gamma }{\tau }\frac{d\gamma }{dv},$$
$`(14)`$
$$\frac{du^k}{ds}=\frac{du^k}{dv}\frac{dv}{ds}=\frac{\gamma }{\tau ^2}\frac{d}{dv}\left(\gamma \frac{dx^k}{dv}\right),$$
$`(15)`$
where $`\gamma `$ is given by Eq. (10).

## III Acceleration and cosmic distances

Multiplying Eq. (4) by $`\text{a}_0^2\tau ^4`$, where $`\text{a}_0`$ is the ordinary three-vector acceleration as measured in the cosmic frame of reference (see Section 2.5 of Ref. 1) at cosmic time $`t=0`$ (i.e. now). We obtain

$$\text{a}_0^2\tau ^6\left(\frac{dv}{ds}\right)^2\left(1\frac{t^2}{\tau ^2}\right)=\text{a}_0^2\tau ^4.$$
$`(16)`$
Using now Eq. (5) in Eq. (16) the latter then yields

$$\left(1\frac{t^2}{\tau ^2}\right)^1\left(\text{a}_0^2\tau ^4\text{a}_0^2t^2\tau ^2\right)=\text{a}_0^2\tau ^4.$$
$`(17)`$
We now define the acceleration a at an arbitrary cosmic time $`t`$ by

$$\text{a}=\frac{\text{a}_0}{\sqrt{1{\displaystyle \frac{t^2}{\tau ^2}}}},$$
$`(18)`$
then Eq. (17) will have the form

$$\text{a}^2\tau ^4\text{a}^2t^2\tau ^2=\text{a}_0^2\tau ^4,$$
$`(19)`$
or

$$\text{a}^2\tau ^4\tau ^2\text{v}^2=\text{a}_0^2\tau ^4,$$
$`(20)`$
where $`\text{v}=\text{a}t`$ is the ordinary three-dimensional velocity. Equation (20) is the analog to

$$m^2c^4c^2\text{p}^2=m_0^2c^4,$$
$`(21)`$
in ESR, where $`\text{p}=m\text{v}`$ is the linear momentum (see Section 7.2 of Ref. 1). Thus a reduces to $`\text{a}_0`$ when the cosmic time $`t=0`$ (i.e. at present).

The comparable to Eq. (18) in ESR is, of course,

$$m=\frac{m_0}{\sqrt{1{\displaystyle \frac{v^2}{c^2}}}},$$
$`(22)`$
where $`m`$ and $`m_0`$ are the inertial mass and the rest mass of the particle, with $`m`$ reduces to $`m_0`$ at $`v=0`$. If we multiply Eq. (22) by $`c^2`$ and expand both sides in $`v/c`$, we obtain

$$mc^2=m_0c^2+\frac{m_0}{2}v^2+\mathrm{}.$$
$`(23)`$
Doing the same with Eq. (18) but multiplication by $`\tau ^2`$, and expanding in $`t/\tau `$, we obtain

$$\text{a}\tau ^2=\text{a}_0\tau ^2+\frac{\text{a}_0}{2}t^2+\mathrm{}.$$
$`(24)`$
Equation (24) in CSR is of course the analog to Eq. (23) in ESR.

## IV Energy in ESR versus cosmic distance in CSR

While Eq. (23) yields

$$E=E_0+\frac{m_0}{2}v^2+\mathrm{},$$
$`(25)`$
Eq. (24) gives

$$\text{S}=\text{S}_0+\frac{\text{a}_0}{2}t^2+\mathrm{}.$$
$`(26)`$
In the above equations $`E`$ is, of course, the energy of the particle, whereas S is the cosmic distance. $`E_0=m_0c^2`$ is the rest energy of the particle whereas $`m_0v^2/2`$ is the Newtonian kinetic energy. What about the terms in Eq. (26)? The term $`\text{a}_0t^2/2`$ is, of course, the Newtonian distance the particle makes due to the acceleration $`\text{a}_0`$, the term $`\text{S}_0=\text{a}_0\tau ^2`$ is unique to CSR, and it might be called the intrinsic cosmic distance of the particle.

Equation (20) can now be written as

$$\text{S}^2\tau ^2\text{v}^2=\text{S}_0^2,$$
$`(27)`$
in complete analogy to

$$E^2c^2\text{p}^2=E_0^2$$
$`(28)`$
in ESR.

## V Cosmic distance-velocity four-vector

We now define the cosmic distance-velocity four-vector. It is the analogous to the energy-momentum four-vector in ESR (see Section 7.4 of Ref. 1). It is defined by

$$\text{v}^\mu =\text{a}_0\tau ^2u^\mu ,$$
$`(29)`$
where $`u^\mu `$ has been defined in Section 2. We have

$$\text{v}^0=\text{a}_0\tau ^2u^0,$$
$`(30a)`$
$$\text{v}^k=\text{a}_0\tau ^2u^k,$$
$`(30b)`$
where $`u^0`$ and $`u^k`$ are given by Eqs. (7)–(10), with

$$\text{v}_0=\text{v}^0,\text{v}_k=\text{v}^k.$$
$`(31)`$
Accordingly we have

$$\text{v}^2=\text{v}_0\text{v}^0+\text{v}_k\text{v}^k=\text{a}_0^2\tau ^4\left(u_0u^0+u_ku^k\right)$$

$$=\text{a}_0^2\tau ^4u_\alpha u^\alpha =\text{a}_0^2\tau ^4=\text{S}_0^2.$$
$`(32)`$
But, using Eqs. (7), (9) and (18),

$$\text{v}^0=\text{a}_0\tau ^2\gamma =\text{a}\tau ^2=\text{S},$$
$`(33a)`$
$$\text{v}^k=\text{a}_0\tau \gamma \frac{dx^k}{dv}=\text{a}\tau \frac{dx^k}{dv}=\text{a}\tau \frac{dx^k}{dt}\frac{dt}{dv}=\tau \text{v},$$
$`(33b)`$
where v is the three-dimensional velocity. Hence

$$\text{v}_\alpha \text{v}^\alpha =\text{v}_0\text{v}^0+\text{v}_k\text{v}^k=\text{a}_0^2\tau ^4\gamma ^2\text{a}_0^2\tau ^2\gamma ^2\left(\frac{dx^k}{dv}\right)^2$$

$$=\text{a}^2\tau ^4\tau ^2\text{v}^2=\text{S}^2\tau ^2\text{v}^2,$$
$`(34)`$
and accordingly, using (32),

$$\text{S}^2\tau ^2\text{v}^2=\text{S}_0^2,$$
$`(35)`$
which is exactly Eq. (20) with $`\text{S}=\text{a}\tau ^2`$ and $`\text{S}_0=\text{a}_0\tau ^2`$. The above analysis also shows that Eq. (35) is covariant under spacevelocity cosmological transformation.

Equation (35) is the analog to

$$E^2c^2\text{p}^2=E_0^2$$
$`(36)`$
in ESR, where $`E`$ and $`E_0`$ are the energy and rest energy, respectively, S and $`\text{S}_0`$ are the cosmic distances at cosmic time $`t`$ and present time ($`t=0`$), respectively.

Finally, in ESR when the rest mass is zero (like the photon), one then has

$$E=cp,$$
$`(37)`$
which is valid at the light cone (see Chapter 6 of Ref. 1). In the case of CSR one also has, when $`\text{S}_0=0`$,

$$\text{S}=\tau \text{v},$$
$`(38)`$
now valid at the galaxy cone (see Section 2.14 of Ref. 1).

## VI Conclusions

In this paper it has been shown that the comparable quantity to the mass in Einstein’s special relativity is the ordinary acceleration three-vector in cosmological special relativity. They both have similar behavior, one with respect to $`v/c`$ and the other with $`t/\tau `$:

$$m=\frac{m_0}{\sqrt{1{\displaystyle \frac{v^2}{c^2}}}},$$
$`(39)`$
and

$$\text{a}=\frac{\text{a}_0}{\sqrt{1{\displaystyle \frac{t^2}{\tau ^2}}}}.$$
$`(40)`$
Furthermore, the role of the energy in Einstein’s theory is being taken over by the cosmic distance,

$$E=mc^2,$$
$`(41)`$
and

$$\text{S}=\text{a}\tau ^2.$$
$`(42)`$
Finally, the analog of the energy formula

$$E^2c^2\text{p}^2=E_0^2$$
$`(43)`$
in ordinary special relativity is

$$\text{S}^2\tau ^2\text{v}^2=\text{S}_0^2$$
$`(44)`$
in cosmological special relativity.


text:
# Geometry of frictionless and frictional sphere packings

## I Introduction

Dense amorphous packings of frictionless spheres have proven to be an extremely useful paradigm in different physical contexts, such as metallic glasses, colloidal crystals, and emulsion rheology. Granular materials are another example of a system with macroscopically large particles, with one major difference: grain-grain interactions involve frictional forces. As a result, granular packings may be quite different from frictionless sphere packings in ways which may impact significantly on their physical properties.

A common quantity of interest in packings of hard spheres is the average number of contacts per particle (coordination number) $`z`$. In order to achieve static mechanical equilibrium, each sphere in the packing needs a sufficient number of constraints that freeze out its translational and rotational degrees of freedom. These constraints are provided by contacts, and once there are a sufficient number of them, the packing can accommodate external body or boundary forces, as long as a set of contact forces satisfying mechanical equilibrium can be found for the given arrangement of such contacts. The minimal average coordination number required to obtain static packings of $`d`$-dimensional frictionless spheres that are stable against external perturbations is $`z_n=2d`$ , whereas for spheres with friction, $`z_f=d+1`$ . In three dimensions, $`z_n=6`$ and $`z_f=4`$. We call such packings “isostatic”.

In this study, we investigate whether or not sphere packings readily achieve isostaticity under generic packing conditions. This isostaticity hypothesis is important in theories focusing on the macroscopic response of such packings,and at first appears reasonable, given the strong numerical evidence from simulation  that $`z_n=6`$ for dense random packings of frictionless spheres. Simulation studies of frictional spheres compressed in a gravity-free environment have shown that $`z_f`$ is significantly less than 6, but with the lowest achieved value of around 4.5, it remains unclear whether the minimal value of 4 is reached in the limit of zero confining pressure (equivalent to the hard-sphere limit.) It is also unclear in what way the packings would change for arbitrarily small friction coefficient $`\mu `$ between the spheres in order to achieve a reduction in $`z`$ to four, if the isostaticity hypothesis were true.

To address such questions, we perform a systematic simulation study of the effect of various parameters on sphere packings. In particular, we vary the following materials properties: the sphere hardness $`k_n`$, the coefficient of restitution $`ϵ`$, and the friction coefficient $`\mu `$ between two particles. We also vary the initial conditions of the packing by varying the initial packing density $`\varphi ^i`$, as well as the initial velocities of the spheres. We investigate how the density, coordination number and the nature of the contacts change as these parameters are varied. Although frictionless hard-spheres appear to form isostatic packings regardless of construction history and restitution coefficient, frictional hard-spheres achieve a multitude of hyperstatic packings $`(z>z_f)`$ that depend on system parameters and construction history . The coordination number reduces smoothly from $`z=6`$ as the friction coefficient is increased, disagreeing with the isostaticity hypothesis.

## II Simulation Method

We present molecular dynamics (MD) simulations in three dimensions on model systems of $`N=20,000`$ monodisperse, cohesionless spheres of diameter $`d`$ and mass $`m`$. The system is spatially periodic in the $`xy`$plane, with a unit cell of size $`20d\times 20d`$, and is bounded in the $`z`$direction by a rough bed at the bottom and an open top. The starting configurations consist of randomly positioned non-overlapping spheres, with packing fractions in the range $`0.02<\varphi ^i<0.3`$, obtained by varying the overall height while keeping $`N`$ and the $`x,y`$ dimensions fixed. The system is subsequently allowed to settle under gravity on top of the rough bed (see Fig. 1). The equilibrated static packing height is about $`50d`$. This method of construction mimics the pouring of granular materials through a sieve to an area far away from side walls, without forming a conical heap.

Different ways of preparing static granular packings include compressing them, and reducing the inclination angle of gravity-driven chute flows below the angle of repose. In the frictionless case, conjugate gradient methods have been used to study dense random packings . However, for the case of particles with friction, MD simulations are more appropriate in order to properly account for both the normal and tangential forces, since the latter depend on the loading history of the contact.

The spheres interact only on contact through a linear spring–dashpot interaction law  in the normal and tangential directions to their lines of centers . Contacting spheres $`i`$ and $`j`$ positioned at $`𝐫_i`$ and $`𝐫_j`$ experience a relative normal compression $`\delta =|𝐫_{ij}d|`$, where $`𝐫_{ij}=𝐫_i𝐫_j`$, which results in a force

$$𝐅_{ij}=𝐅_n+𝐅_t.$$
(1)
The normal and tangential contact forces are given by

$$𝐅_n=k_n\delta 𝐧_{ij}\frac{m}{2}\gamma _n𝐯_n,$$
(2)
$$𝐅_t=k_t\mathrm{\Delta }𝐬_t\frac{m}{2}\gamma _t𝐯_t,$$
(3)
where $`𝐧_{ij}=𝐫_{ij}/r_{ij}`$, with $`r_{ij}=|𝐫_{ij}|`$, $`𝐯_n`$ and $`𝐯_t`$ are the normal and tangential components of the relative surface velocity, and $`k_{n,t}`$ and $`\gamma _{n,t}`$ are elastic and viscoelastic constants respectively. $`\mathrm{\Delta }𝐬_t`$ is the elastic tangential displacement between spheres, obtained by integrating surface relative velocities during elastic deformation of the contact. The magnitude of $`\mathrm{\Delta }𝐬_t`$ is truncated as necessary to satisfy a local Coulomb yield criterion $`F_t\mu F_n`$, where $`F_t|𝐅_t|`$ and $`F_n|𝐅_n|`$. Frictionless spheres can be simulated simply by setting $`\mu =0`$. Finally, the particles are moved according to the total forces and torques applied to them through contacts and the gravitational field. Additional details can be found in Ref..

The presented simulations were carried out for a range of materials parameters. The normal spring constant $`k_n`$ varied from $`210^5`$ to $`210^9mg/d`$ in order to understand how the system behaves as it approaches the hard sphere limit. In all cases, $`k_t=2k_n/7`$ . In order to study the crossover from frictionless to frictional systems, the local particle friction coefficient $`\mu `$ was varied from $`0`$ to $`10`$. Finally the coefficient of restitution $`ϵ`$, i.e. the ratio of the final to initial normal velocities in a head-on binary collision, was also varied. For a linear spring–dashpot interaction,

$$ϵ=\mathrm{exp}(\gamma _nt_{\mathrm{col}}/2),$$
(4)
where the collision time $`t_{\mathrm{col}}`$,

$$t_{\mathrm{col}}=\pi (2k_n/m\gamma _n^2/4)^{1/2}.$$
(5)
Three values of $`ϵ`$=$`0.26`$, $`0.50`$, and $`0.88`$ were used. The effect of $`ϵ`$ is like that of varying the quench rate in a thermal system: for smaller values of $`ϵ`$, collisions are more inelastic, energy is dissipated faster, and the spheres have less ability to move from the the first position where they are mechanically stable, which may result in packings that are less stable than similarly prepared packings of more elastic grains.

Most of the simulations were started from a static configuration of $`\varphi ^i0.13`$ as shown in Fig. 1(a). The timestep $`\delta t0.05\sqrt{m/k_n}`$ was chosen to accommodate the decreasing collision time as the particle hardness $`k_n`$ is increased \[cf. Eq. (5)\]. For $`k_n=210^5mg/d`$, $`\delta t10^4\sqrt{d/g}`$. Simulations were then run until the kinetic energy per particle was less than $`10^8mgd`$ for small $`k_n`$, and up to three orders of magnitude less for large $`k_n`$. This requires $`3810^6\delta t`$ for small $`k_n`$, and $`4810^7\delta t`$ for $`k_n=210^9mg/d`$.

Because of the large amount of time required to reach static equilibrium, each set of parameters was run for 1 to 5 configurations. In some cases, the particles were given an initial, random velocity with an average KE per particle of approximately $`20100mgd`$. The results were identical in all cases within sample to sample fluctuations.

## III Results

### A Coordination Number

Figure 2 shows the effect of friction coefficient $`\mu `$ on $`z`$ for $`k_n=210^5mg/d`$ and $`ϵ=0.88`$ and $`0.26`$. For both values of $`ϵ`$, $`z=6.144\pm 0.002`$ for frictionless packings. As we shall see in Sec. III C, the deviation from the isostatic value of 6 can be attributed to the finite stiffness of the spheres, the isostatic value is apparently obtained in the hard-sphere limit. However, there is no sudden drop from $`z=6`$ as friction is turned on; rather there is a gradual decrease in $`z`$ to a parameter dependent minimal value, accompanied by a similar decrease in the final volume fraction $`\varphi ^f`$ (see Fig. 3.)

As depicted in Fig. 4, the decrease in $`z`$ is primarily due to an overall shift in the distribution of coordination numbers to lower values, rather than a change in its shape and width. Consequently, the frequency of particles with eight or more neighbours reduces as $`\mu `$ increases, and particles with as few as two contacts start to appear at $`\mu =0.5`$, indicative of arching within the packing at large $`\mu `$. The saturation of $`z`$ and $`\varphi ^f`$ for $`\mu 1`$ is due to the fact that the typical tangential forces $`F_t`$ in a packing with $`\mu =\mathrm{}`$ is expected to be of order $`F_n`$, and lowering $`\mu `$ has little effect on the packing down to $`\mu 1`$. This behavior of the contact forces is further verified in Sec. III D.

Unlike the frictionless case, the deviations from isostaticity $`(z=4)`$ for $`\mu >0`$ cannot be attributed to corrections due to the finite stiffness of the spheres. The packings remain unambiguously hyperstatic in the hard-sphere limit.

When a static packing with a surface tilt near the angle of repose was generated by cessation of flow down an inclined plane, similar results for $`z`$ and $`\varphi `$ were obtained. For $`k_n=210^5mg/d,ϵ=0.88,\mathrm{and}\mu =0.50`$, such packings gave $`z=4.69`$ and $`\varphi =0.594`$, compared to $`z=4.90`$ and $`\varphi =0.61`$ for packings presented in this study.

### B The Radial Distribution Function

The radial distribution function (RDF), $`g(r)`$, for $`k_n=210^5mg/d`$ and $`ϵ=0.88`$, is plotted in Fig. 5 for several values of $`\mu `$. The characteristic split second peak, indicating short-range order out to second neighbors, is evident. For $`\mu =0`$, $`g(r)`$ is essentially identical to that obtained for random close packing (RCP), at volume fraction $`\varphi ^f0.64`$. As $`\mu `$ increases, the secondary peaks in $`g(r)`$ diminish, as seen in the inset.

The first peak of $`g(r)`$ is of particular interest, since near-contacts with $`r/d`$ just over 1 play an important role in the dependence of $`z`$ on the stiffness of the spheres (See Sec. III C). Figure 6 reveals a square-root singularity of the RDF near $`r/d=1`$, i.e.,

$$g(r)\left(\frac{r}{d}1\right)^\alpha ,0<\frac{rd}{d}1,$$
(6)
with $`\alpha =0.52\pm 0.03`$. This singularity has apparently not been reported elsewhere, although we have also verified its presence in the RDF of hard-sphere packings provided to us. Note that this singularity is integrable, and is distinct from particles actually in contact.

### C The Hard-Sphere Limit

We next investigate the effect of the finite stiffness of the spheres on the packings, and on the average coordination number $`z`$ in particular. Let us assume that the packings formed by these stiff elastic spheres are statistically equivalent to packings that would be obtained by first forming a truly hard-sphere packing and subsequently allowing elastic relaxation. Due to the slight compression of the spheres of finite stiffness under gravity, we expect to see an increase in coordination number during this elastic relaxation. Since the typical compressive strain of a sphere under the same loading conditions scales as $`1/k_n`$, we expect that a finite fraction of neighbors in the hard-sphere packing that were within a distance $`d[1+O(mg/k_nd)]`$ of each other form new contacts upon elastic relaxation. The number of such near-contacts in the hard-sphere packing can be computed by integrating the hard-sphere RDF $`g_{\mathrm{}}(r)`$ over the range $`r/d(1,1+mg/k_nd)`$, yielding an effective coordination number

$$z_n(k_n)=z_{\mathrm{}}+\stackrel{~}{a}_n\left(\frac{k_nd}{mg}\right)^{\alpha _n},\frac{k_nd}{mg}1,$$
(7)
where $`z_{\mathrm{}}`$ is the coordination number in the hard-sphere limit, $`\stackrel{~}{a}_n`$ is a constant, and the exponent \[cf. Eq.(6)\]

$$\alpha _n=1\alpha .$$
(8)
Thus, there is a power-law correction to the apparent coordination number, with an exponent that depends on the nature of the singularity in the first peak of $`g(r)`$. A numerical fit of the data to Eq. (7) shown in Fig. 7(a) results in $`\alpha _n=0.498\pm 0.002`$ and $`z_{\mathrm{}}=6.01\pm 0.02`$, in excellent agreement with the isostaticity hypothesis, as well as the exponent relation Eq. (8). Makse et al. also found that $`z`$ approaches 6 as the stress goes to zero in their numerical studies of compressed spheres.

Armed with this insight, we apply a similar analysis to frictional packings, with results presented in Fig. 7(b). Although the RDF of frictional packings appears to have the same square-root divergence near $`r=d`$ (see Fig. 6), the numerical fit to Eq. (7) in the presence of friction yields a different exponent $`\alpha _f1/4`$, resulting in a slower approach to the hard-sphere limit. Thus, the exponent identity Eq. (8) does not hold for frictional spheres. Moreover, in contrast to the frictionless case, the hard-sphere limits remain firmly above the isostatic value of four, and vary as a function of $`\mu `$ and $`ϵ`$.

Even though the rather unlikely scenario of a further crossover to isostaticity at extreme stiffnesses cannot be entirely ruled out, it must be pointed out that the stiffest spheres with $`k_n=210^9mg/d`$ in these packings experience strains $`\delta /d10^8`$. This should be compared to the strain of a “typical” grain, i.e. a glass sphere with a 100 micron diameter, under just its own weight on the Earth’s surface: $`\delta /d(\rho gd/E)^{2/3}`$. With the Young’s Modulus $`E610^{10}`$ Pa, $`\rho 210^3\mathrm{kg}/\mathrm{m}^3`$ and $`g10`$ m/s2, this strain is about $`10^7`$. Thus, even if isostaticity is ultimately restored, the relevance of the hard-sphere limit for real granular systems is still questionable.

### D Plasticity of Contacts

One potential explanation for the hyperstaticity of these frictional packings is the loss of degrees of freedom for the tangential forces in contacts that have become “plastic” such that $`F_t=\mu F_n`$. If a finite fraction of the contacts satisfied this condition, the isostaticity condition would need to be modified. However, as demonstrated by the distribution of the plasticity index $`\zeta \frac{F_t}{\mu F_n}`$ of the contacts in Fig. 8, almost all contacts in the static packings are below their frictional threshold $`\zeta =1`$, eliminating this possibility. A similar distribution of $`\zeta `$ was observed for a static packing created from flow arrest .

For $`\mu 1`$, the distribution of the contact force ratio $`F_t/F_n`$ indeed becomes independent of $`\mu `$, which manifests itself as a collapse of $`P(\zeta )`$ when plotted against $`F_t/F_n`$ (not shown). This result is in accord with the observation in Sec. III A that packings in this range of $`\mu `$ behave effectively the same as for systems with $`\mu =\mathrm{}`$.

### E Effect of dissipation and initial conditions

The dependence of the packing geometry on coefficient of restitution $`ϵ`$ is also interesting. Unlike $`k_n`$ and $`\mu `$, changing $`ϵ`$ would not change the configuration of a static packing after it has stopped – it only affects the relaxation dynamics by increasing the removal rate of kinetic energy. In this sense, changing $`ϵ`$ is like changing the quench rate of a supercooled liquid as it undergoes a glass transition. For very large quench rates, the system might be expected to stop immediately upon forming the minimum number of contacts necessary to achieve static mechanical equilibrium.

Similarly to the effect of $`ϵ`$, we find that the initial starting densities also affect the final packing. In Fig. 9 we plot the variation in the final packing fraction $`\varphi ^f`$, and coordination number $`z`$, as a function of the starting density $`\varphi ^i`$. We find that more dilute starting states lead to more compact final states. This behaviour may be due to the increase in potential energy the system receives when it is more dilute, converting into kinetic energy of the particles during settling, and enabling them to explore more of the phase space on their way to achieving a preferred packing. An empirical fit to the packing fraction is given by,

$$\varphi ^f=0.5778+0.0567\mathrm{exp}(4.3\varphi ^i),$$
(9)
which is similar to the empirical fit for model 2D systems proposed in Ref. . However, it should be noted that for frictional spheres, extrapolation of the coordination number to the limit $`\varphi ^i=\varphi ^f`$ does not result in isostaticity as observed in Ref. .

In light of the dependence of the final state on such parameters as $`ϵ`$ and $`\varphi ^i`$, the breakdown of the exponent identity Eq.(8) for frictional packings is perhaps not surprising. Packings obtained by the sedimentation of hard-spheres followed by elastic relaxation are probably not statistically equivalent to packings of particles that are elastic from the start, due to the strong history-dependence of the final states obtained.

In general, for hyperstatic packings the force network is not uniquely determined by the packing and the loadings on the particles. It follows that the determination of the force network for hyperstatic packings of perfectly rigid particles is an ill-posed problem. Thus the order of the limit $`k_n\mathrm{}`$ and the preparation of the packing cannot be commuted for frictional spheres. More bluntly, perfectly rigid sphere systems with friction are unlikely to capture the mechanical properties of packings of frictional, elastic spheres, even in the limit of extremely large rigidity of these latter particles.

## IV conclusions

We have studied large scale packings of spherical grains of varying hardness, friction coefficient, and coefficient of restitution, formed by sedimentation. We accounted for the systematic variation with particle stiffness and were able to infer properties of hard-sphere packings. Although frictionless hard-spheres appear to form isostatic packings regardless of construction history and restitution coefficient, frictional packings achieve a multitude of hyperstatic packings that depend on system parameters and construction history . The coordination number reduces smoothly from $`z=6`$ as the friction coefficient is increased, contrary to the hypothesis of isostaticity in such packings.

###### Acknowledgements.

Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under Contract DE-AC04-94AL85000. DL acknowledges support from US-Israel Binational Science Foundation Grant 1999235.


text:
# Binding energy in two and three-body relativistic dynamics 11footnote 1Talk presented at the Workshop ”Critical Stability of Few-Body Quantum Systems”, Les Houches, Oct.8-12, 2001

## 1 Introduction

In galilean mechanics, the binding energy of a bound state is just the energy of relative motion. It is an observable and can be expressed in terms of dynamical variables.

In the framework of special relativity, in contrast, the binding energy of a bound state is usually defined as the mass defect $`mM`$ where $`M`$ is the total mass and $`m`$ refers to the constituent masses. But it is still natural to look for some observable, depending on the relative degrees of freedom, that could precisely characterize the relative motion of the system. Of course, in case we succeed, a comparison of this quantity with $`mM`$ is desirable.

There are many different formulations of few-body relativistic dynamics  . Our approach is based upon coupled wave equations   that can be interpreted as mass-shell constraints . For scalar particles, these constraints deal with squared-mass operators that include an interaction term. In contrast, all the generators of the Poincaré algebra are free of interaction. In this report, we focus on scalar particles.

For two (resp. three) particles we start from a pair (resp. a triple) of coupled equations, for one wave function depending on two (resp. three) four-dimensional arguments. In this formulation, like in the Bethe-Salpeter (BS) equation, covariance is paid by redundant degrees of freedom.

But demanding that the total linear momentum four-vector $`P=p`$ has a sharp value $`k`$, and solving the difference(s) of the constraint equations, it is possible to factorize out one (resp. two) degrees of freedom, and one is left with a reduced wave equation, in the form of an eigenvalue problem; solving this problem determines the total mass $`M=\sqrt{kk}`$.

This procedure is straightforward for two-body systems , but in the three-body case, it requires a nontrivial transformation .

In both cases, by an elementary combination of the mass-shell constraints, one can exhibit a constant of the motion, say $`N=\mathrm{\Lambda }`$ which is homogeneous to a squared mass, depends only on relative variables (therefore is translation invariant), includes the term of mutual interaction, and is intimately related with relative motion. Our purpose is to clarify the relationship between this observable and the mass defect.

It seems natural to define slow relative motion by the condition that the eigenvalue $`\lambda `$ of $`\mathrm{\Lambda }`$ is small (in absolute value) with respect to the squares of the constituent masses, $`m_a`$ (this situation in turn implies that also $`|\lambda |M^2`$).

Notation

The dynamical variables associated with particle $`a`$ are the four-dimensional conjugate quantities $`q_a,p_a`$. Particle labels are $`a,b=1,2`$ (resp. $`a,b=1,2,3`$).

## 2 Two-Body System

For the class of models that we consider, the same interaction term $`2V`$ arises in both sqared-mass operators. This scheme seems to accomodate practically all realistic interactions inspired from quantum field theory. Its relationship with either the quasi-potential approach   or the BS equation  has been demonstrated, and fermionic generalizations have been elaborated. We introduce relative variables; they are the four-vectors

$$z=q_1q_2,y=\frac{1}{2}(p_1p_2)$$
(1)
We have the transverse parts

$$\stackrel{~}{z}=z\frac{zP}{P^2}P,\stackrel{~}{y}=y\frac{yP}{P^2}P,$$
(2)
For arbitrary masses it is convenient to set

$$\mu =\frac{1}{2}(m_1^2+m_2^2),\nu =\frac{1}{2}(m_1^2m_2^2)$$
(3)
The reduced wave equation takes on the form

$$(N+\lambda )\varphi =0$$
(4)
with

$$N=\stackrel{~}{y}^2+2V$$
(5)
$$\lambda =\frac{M^2}{4}+\frac{\nu ^2}{M^2}\mu $$
(6)
The interaction term $`V`$ has the dimension of a squared mass. It may essentially depend on $`\stackrel{~}{z}^2`$ and $`P^2`$ (possibly also on $`\stackrel{~}{y}^2,\stackrel{~}{z}\stackrel{~}{y},yP`$). On shell and using the rest frame, we simply have $`\stackrel{~}{z}^2=𝐳^2`$. Therefore (4) can be identified with a nonrelativistic equation , at least when the ”potential” $`V`$ doesnot depend on $`P^2`$. But realistic interactions may bear some dependence on $`P^2`$ (energy dependent potentials); in this case, (4) becomes a nonconventional eigenvalue problem. This complication has been discussed in the literature  .

Weak Binding

Eq. (6) can be solved for $`M^2`$. Insofar as $`\nu `$ is not too large, we obtain

$$M^2=2(\lambda +\mu )+2\sqrt{(\lambda +\mu )^2\nu ^2}$$
(7)
For slow relative motion, $`|\lambda |m_a^2`$, we can develop $`M^2`$ in powers of $`\lambda `$. According to (3) we get

$$M^2=(m_1+m_2)^2+\lambda \frac{(m_1+m_2)^2}{m_1m_2}+O(\lambda ^2)$$

hence

$$m_1+m_2M=\frac{\lambda }{2m_0}+\mathrm{}$$
(8)
It is clear that slow relative motion corresponds to weak binding.

## 3 Three-Body System

The squared mass operators are $`p_a^2+2V`$. This formulation aims at the elimination of two superfluous degrees of freedom.

For three-body systems it is difficult to find an interaction such that the mass-shell constraints are compatible among themselves, respect Poincaré invariance, reduce to three Klein-Gordon equations in the absence of interaction, and allow for eliminating two redundant degrees of freedom. All these requirements can be satisfied however, in a tractable manner, with help of a ”point transformation in momentum space”. In order to formulate this transformation, it is essential to introduce relative variables as follows .

Relative-particle indices are $`A,B=2,3`$. We define the four-vectors

$$z_A=q_1q_A,y_B=\frac{P}{3}p_B$$
(9)
The transverse part are $`\stackrel{~}{z}_A,\stackrel{~}{y}_B`$ and $`\widehat{z},\widehat{y}`$ respectively with respect tp $`P`$ and $`k`$. In the rest frame we have $`\widehat{y}_A^2=𝐲_A^2,\widehat{z}_A^2=𝐳_A^2`$, etc.

Our transformation  is characterized by

$$(p_1+p_A)(p_1p_A)=P(p_1^{}p_A^{})$$
(10)
and by the requirement that it leaves $`P`$ and $`\stackrel{~}{y}_A`$ invariant (the new relative momenta are of course $`y_A^{}=P/3p_A^{}`$ ).

This procedure generates a canonical transformation, giving rise to new configuration variables $`z_A^{}`$. The difference equations are mapped to $`y_A^{}P=c_A\mathrm{\Psi }`$ where the constants $`c_A`$ are combinations of the squared constituent masses; the dependence of $`\mathrm{\Psi }`$ on the new relative times is factorized out.

Here, for simplicity, we assume that $`m_a=m`$.

The sum of the mass-shell constraints yields the reduced equation

$$(9m^2M^2)\psi =6N\psi $$
(11)
where now

$$N=\mathrm{\Lambda }=\stackrel{~}{y}_{2}^{}{}_{}{}^{2}+\stackrel{~}{y}_3^2+\stackrel{~}{y}_2\stackrel{~}{y}_3+3V+P^2\mathrm{\Xi }$$
(12)
In the momentum representation, the reduced wave function $`\psi `$ depends only on $`\widehat{y}_2,\widehat{y}_3`$. In view of the compatibility requirement, a closed form of the interaction is available only in terms of the new variables. A typical example would be a function of $`\widehat{z}_2^{},\widehat{z}_3^{},P^2`$.

The quantity $`P^2\mathrm{\Xi }`$ has no counterpart in two-body systems. Here it stems from having added three constraints. Its exact expression in terms of the new variables amounts to solve a fourth-degree algebraic equation and would be extremely complicated (see details in ). Fortunately, it can be naturally expanded in powers of $`1/P^2`$, implying on the mass shell an expansion in powers of $`1/M^2`$. This makes our model more tractable when the constituent particles are light with respect to the total mass of the system.

When the constituent masses (although different from zero) are small with respect to the total mass $`M`$, that is $`m^2M^2`$, we can drop the last term in (12).

With this truncation, (11) is similar to a nonrelativistic equation, except for eventual complications resulting from a possible dependence of $`N`$ on $`P^2`$. Moreover, in this limit $`\widehat{z}_A^{}`$ differs very little from $`\widehat{z}_A^{}`$, which allows for a weak form of cluster separability .

But when $`m^2/M^2`$ is not small enough, the on-shell expression of $`P^2\mathrm{\Xi }`$ must be written as a Taylor expansion including several powers of $`1/M^2`$, say $`M^2\underset{¯}{\mathrm{\Xi }}={\displaystyle \frac{1}{M^2}}\mathrm{\Gamma }`$, where $`\mathrm{\Gamma }`$ is regular for $`M\mathrm{}`$, say $`\mathrm{\Gamma }=\mathrm{\Gamma }_{(0)}+O(1/M^6)`$ where

$$\mathrm{\Gamma }_{(0)}=(\widehat{y}_{2}^{}{}_{}{}^{2})^2+(\widehat{y}_{3}^{}{}_{}{}^{2})^2+(\widehat{y}_2\widehat{y}_3)^2+(\widehat{y}_{2}^{}{}_{}{}^{2}+\widehat{y}_{3}^{}{}_{}{}^{2})(\widehat{y}_2\widehat{y}_3)\widehat{y}_{2}^{}{}_{}{}^{2}\widehat{y}_{3}^{}{}_{}{}^{2}$$
(13)
Note that $`\mathrm{\Gamma }`$ is a positive operator and would survive in the absence of interaction. Irrespective of the shape of the potential, we expect that it provides a positive correction to the truncated expression of $`N`$. A rigorous statement, however, would require solving a nonconventional eigenvalue problem.

At first order in $`1/M^2`$ we have, in the rest frame

$$N\psi =(𝐲_2^2𝐲_3^2𝐲_2𝐲_\mathrm{𝟑}+3V+\frac{\mathrm{\Gamma }_{(0)}}{M^2})\psi $$

In the rest-frame $`\mathrm{\Gamma }_{(0)}`$ is bi-quadratic in $`𝐲`$.

Let us evaluate the binding energy, in the general case described by (11). Equation (11) can be viewed as a (generalized) eigenvalue equation for $`\mathrm{\Lambda }`$ with eigenvalue $`\lambda `$, if we set $`6\lambda =M^29m^2`$. The mass defect is exactly

$$mM=3m\sqrt{9m^2+6\lambda }$$
(14)
It would increase together with the eigenvalue of $`N`$.

Weak Binding

It is noteworthy that weakly bound systems are not eligible for the light-constituent approximation. Indeed they have $`m^2`$ almost one order of magnitude smaller than $`M^2`$, but this is not enough for dropping terms like $`O(1/M^2)`$.

Assuming that $`|\lambda |m^2`$ wee can develop the mass defect and find

$$3mM=\frac{\lambda }{m}+\frac{1}{6}\frac{\lambda ^2}{m^3}+\mathrm{}$$
(15)
Since we consider equal masses, then $`m=2m_0`$ where $`m_0`$ is the reduced mass of either of particles 2, 3, with respect to particle 1. Note the analogy of (15) with (8).

## 4 Conclusion

The constraint formulation of relativistic two-body dynamics admits a well- understood contact with field theory, whereas constraint three-particle dynamics is still a field of recent investigation. Nevertheless, it is possible to present in parallel ways, for both cases, the relationship of binding energy with a remarkable observable $`N`$ which naturally arises in the reduced wave equation. This situation results from the fact that, in both cases, our basic equations involve a unique interaction term and are tailored for allowing elimination of the redundant variable(s) implied by manifest covariance.

Binding energy, defined as the mass defect, has a simple relationship with the eigenvalue of $`N`$. In the case of weak binding, these quantities become proportional through the constant factor $`1/2m_0`$.

The reduced wave equation can be compared and, to some extent, identified with a nonrelativistic equation: in a straightforward manner for two particles (with arbitrary masses), but only in the light-constituent limit for three particles (with equal masses); in three-body systems, weak binding doesnot correspond to a nonrelativistic form of the wave equation.

For applications, we plan to introduce spin and to improve the contact with other approaches .